function [h] = lyngby_fir_main(x, Y, arg1, arg2, arg3, arg4, arg5, ...
    arg6, arg7, arg8, arg9, arg10, arg11, arg12);

% lyngby_fir_main      - Regularized FIR filter, main function 
%
%       function H = lyngby_fir_main(x, Y, 'PropertyName',
%          'PropertyValue');
%
%       Input:    x   Input signal to the system - the paradigm (vector)
%                 Y   Output signal of the system (vector or matrix). 
%
%       Property: FilterOrder      {7} The order of the estimated FIR
%                                  filter (integer), ie, filter length 
%                 RegMethod        [ {Ridge} | RegInverse | PCR | SVD
%                                  | Smooth ] Ridge regression
%                                  (ridge/reginverse) or principal
%                                  component regression (PCR/SVD) of
%                                  smooth FIR 
%                 PCRComp          Number of principal components for
%                                  the PCR 
%                 Regularization   {0} Regularization parameter for
%                                  ridge regularization 
%                 ConvType         [ {Skip} | Zeropad | Extend | Wrap ]
%                                  Convolution type. if 'Zeropad' use
%                                  all data points in the filter
%                                  estimation. if 'Skip' it will
%                                  disregard the FilterOrder-1 first
%                                  data points. 'Extend' will use the
%                                  first value of 'x': x(t) = x(1), t<1
%
%       Output:   H   Estimated response function
%
%       Estimation of the finite impulse response (FIR) filter, ie, a
%       linear model with 'stick' functions. This is also a "one-layer
%       feedforward linear neural network with weight decay" and an
%       "ARX(0,n) model (with ridge regression)" 
%
%       There are several regularization methods control with the
%       'RegMethod' property:
%
%       With the 'RegMethod' as 'Ridge' the regularization will
%       be ridge regression, and the 'Regularization' parameter will
%       be the ridge parameter. 
%
%       The 'PCR' regularization method will first preprocess the input
%       making basis function out of the input and use this in the
%       regression. The regularization parameter will then function
%       as a threshold parameter cutting off SVD component with a
%       singular value lower than this threshold.
%
%       Ref: Goutte et al., (2000), IEEE TMI, 19(12):1188+.
%
%       See also LYNGBY, LYNGBY_FIR_ERROR, LYNGBY_FIR_CONVOLVE,
%                LYNGBY_LZIT_MAIN, LYNGBY_POIS_MAIN.
%
% $Id: lyngby_fir_main.m,v 1.18 2003/02/21 13:39:52 fnielsen Exp $



    % Check arguments
    [DataSize, NoVoxels] = size(Y);
    if length(x) ~= DataSize,
       error(sprintf([ 'Inputs do not match in size. ' ...
	     'Y: %dx%d, x: %dx%d'], size(Y), size(x)));
    end
    
    % Default parameters
    filterOrder = 7;
    reg = 0;
    PCRComp = 3;
    regMethod = 1;
    convType = 1;
    
    % Parse Properties
    n = 1;
    while n <= nargin-2
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'filterorder')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if ~isstr(arg)
	  if length(arg) == 1
	    if arg > 0
	      filterOrder = arg;
	    else
	      error(sprintf(['''FilterOrder'' should be larger than 0. ' ...
		    'It was: %d'], arg))
	    end
	  else
	    error(['The argument with ''FilterOrder'' should be ' ...
		  'a scalar.']); 
	  end
	else
	  error('The argument with ''FilterOrder'' should be a integer.'); 
	end
      elseif strcmp(arg, 'regularization')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if ~isstr(arg)
	  if length(arg) == 1
	    if arg >= 0
	      reg = arg;
	    else
	      error(sprintf(['''Regularization'' should be >= 0. ' ...
		    'It was: %d'], arg))
	    end
	  else
	    error(['The argument with ''Regularization'' should be ' ...
		  'a scalar.']); 
	  end
	else
	  error('The argument with ''Regularization'' should be a integer.'); 
	end
      elseif strcmp(arg, 'pcrcomp')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if ~isstr(arg) & arg > 0 & arg == round(arg) & ...
	      length(arg) == 1 
	  PCRComp = arg;
	else
	  error(['The argument with ''PCRComp'' should be ' ...
		'a positive integer.']); 
	end
      elseif strcmp(arg, 'regmethod') | strcmp(arg, 'inversionmethod')
	% InversionMethod was an old name. 
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'reginverse')
	    % Old name for ridge
	    regMethod = 1;
	  elseif strcmp(arg, 'ridge')
	    regMethod = 1;
	  elseif strcmp(arg, 'pcr')
	    regMethod = 2;
	  elseif strcmp(arg, 'svd')
	    % SVD old name for PCR (principal compenent regression)
	    regMethod = 2;
	  elseif strcmp(arg, 'smooth');
	    regMethod = 3;
	  else
	    error(sprintf('Wrong arg to ''regMethod'': %s', arg))
	  end
	else
	  error('The argument with ''regMethod'' should be a string.'); 
	end
      elseif strcmp(arg, 'convtype')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if isstr(arg)
	  if strcmp(arg, 'skip')
	    convType = 1;
	  elseif strcmp(arg, 'zeropad')
	    convType = 2;
	  elseif strcmp(arg, 'extend')
	    convType = 3;
	  elseif strcmp(arg, 'wrap')
	    convType = 4;
	  else
	    error(sprintf('Wrong arg to ''ConvType'': %s', arg))
	  end
	else
	  error('The argument with ''ConvType'' should be a string.'); 
	end
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    PCRComp = min([PCRComp filterOrder]);
    
    
    h = zeros(filterOrder, NoVoxels);
        
    x = reshape(x,DataSize,1);        % Make sure that it is a column vector
        
          
    if convType == 1
      % Skip the first 
      Rows = DataSize-filterOrder+1;	% Number of equations to determine 
      Xm = zeros(Rows,filterOrder);     % The Input as a matrix
      for r = 1:Rows,
        Xm(r,:) = x(filterOrder+r-1:-1:r)';
      end
    elseif convType == 2
      % Zeropad
      Rows = DataSize;	                % Number of equations to determine 
      Xm = zeros(Rows,filterOrder);     % The Input as a matrix
      x2 = [ zeros(filterOrder-1,1); x ];
      for r = 1:Rows,
	Xm(r,:) = x2(filterOrder+r-1:-1:r)';
      end
    elseif convType == 3
      % Extend
      Rows = DataSize;	                % Number of equations to determine 
      Xm = zeros(Rows,filterOrder);     % The Input as a matrix
      x2 = [x(1)*ones(filterOrder-1,1) ; x];
      for r = 1:Rows,
        Xm(r,:) = x2(filterOrder+r-1:-1:r)';
      end
    elseif convType == 4
      % Wrap 
      Rows = DataSize;	                % Number of equations to determine 
      Xm = zeros(Rows,filterOrder);     % The Input as a matrix
      x2 = [ x(end-filterOrder:end) ; x ];
      for r = 1:Rows,
        Xm(r,:) = x2(filterOrder+r-1:-1:r)';
      end
    else
      error(sprintf('Internal error: convType %d', convType));
    end

    Xmm = Xm'*Xm;
    
    if regMethod == 1
      lyngby_log('Calculating Regularized Inverse');

      R = reg * eye(filterOrder);    % A unit matrix for the regularization
      preX = inv(Xmm + R) * Xm';     % The 'Normal equation' for ridge regression 
      
    elseif regMethod == 2          
      lyngby_log('Calculating SVD');
      [u,s,v] = svd(Xmm);
      % Compute pseudo inverse
      Mask = diag([ones(PCRComp,1) ; zeros(size(s,1)-PCRComp,1)]);
      ss = (1./(s + eps)) .* Mask;   
      preX = v*ss*u'*Xm';
      
    elseif regMethod == 3 
      lyngby_log('Calculating Smooth FIR filter');

      r = lyngby_pdf_gauss((0:filterOrder-1)', 0, (filterOrder/5)^2);
      R = reg * pinv(toeplitz(r,r) + r(1)*0.0001*eye(filterOrder));
      preX = inv(Xmm + R) * Xm';               
    else
      error(sprintf('Internal error, regMethod: %d', regMethod));
      
    end
 
    if convType == 1
      timev = filterOrder:DataSize;
    else
      timev = 1:DataSize;
    end
    
    h = preX*Y(timev,:);  % Can not use the first values






