function R = lyngby_fir_regmatrix(filterOrder, varargin)

% lyngby_fir_regmatrix - Construct regularization matrix
%
%       R = lyngby_fir_regmatrix('PropertyName', 'PropertyValue')
%
%       Input:    filterOrder Number of stick functions
%
%       Property: Reg         [ {0} | non-negative ] Regularization
%                             strength          
%                 RegMethod   [ {smooth} | ridge ] Regularization
%                             method 
%
%       Output:   R           Regularization matrix
%
%       This function produces a regularization matrix that can be
%       used in a biased estimate. It can be used in the multivariate
%       regression model (Y = X*B+U) for estimatrion of B
%            
%          B = (X'*X + R)^(-1) * X'*Y
%
%       The regularization matrix correspond to eith the 'smooth' FIR
%       filter or the ridge regularized FIR filter depending on the
%       setting of the 'RegMethod' property.
%
%       See also LYNGBY, LYNGBY_FIR_MAIN, LYNGBY_FIR_TEST_CHI2. 
%
% $Id: lyngby_fir_regmatrix.m,v 1.2 2003/02/21 14:27:05 fnielsen Exp $


    % Check input arguments
    if nargin < 1
      error('Too few input arguments')
    end


    % Default values
    regMethod = 'smooth';
    reg = 0;


    % Read properties
    n = 1;
    while n < length(varargin)
      arg = lower(varargin{n});

      if strcmp(arg, 'reg') 
	n = n + 1;
	arg = lower(varargin{n});
	if isnumeric(arg) & length(arg)==1 &  arg >= 0
	  reg = arg; 
	else
	  error(sprintf(['Argument with ''Reg'' should be '...
		'a non-negative value']));
	end

      elseif strcmp(arg, 'regmethod') 
	n = n + 1;
	arg = lower(varargin{n});
	if isstr(arg) & ( any(strmatch(arg, {'smooth' 'ridge'}, 'exact')) )
	  regMethod = arg; 
	else
	  error(sprintf(['Argument with ''RegMethod'' should be '...
		'a string either ''smooth'' or ''ridge''.']));
	end

	
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    
    
    if strcmp(regMethod, 'ridge')
      lyngby_log('Calculating Regularized Inverse');

      R = reg * eye(filterOrder);    % A unit matrix for the regularization
      
    elseif strcmp(regMethod, 'smooth');

      r = lyngby_pdf_gauss((0:filterOrder-1)', 0, (filterOrder/5)^2);
      R = reg * pinv(toeplitz(r,r) + r(1)*0.0001*eye(filterOrder));
      
    else
      error('Internal error');
    end
    





