function C = lyngby_xcorr(X, P, arg1, arg2, arg3, arg4, arg5, arg6)

% lyngby_xcorr         - Cross correlation/covariance
%
%       function C = lyngby_xcorr(X, P, ...
%           'PropertyName', 'PropertyValue')
%
%       Input:    X           Datamatrix
%                 P           Paradigm vector
%     
%       Property: Biased      Divide by the number of time steps
%                 Components  Number of cross correlation components
%                             returned. Should be a 
%                             number. By default all components
%                             are returned.
%                 Type        [ biased | coeff | {none} ]
%
%       Output:   C           Cross correlation
%
%       This function returns the cross correlation or the
%       cross-covariance between a matrix X and a vector P.
%
%       The 'Biased' property is redundant in connection with the
%       'Type' property.
%
%       See also LYNGBY, LYNGBY_UI_MAIN.
%
% $Id: lyngby_xcorr.m,v 1.10 2002/05/08 14:25:56 fnielsen Exp $


    % Default values
    Components = size(X, 1);
    type = 'none';
    
    % Pars Properties
    n = 1;
    while n < nargin-1
      eval(sprintf('arg = lower(arg%d);', n)); 

      if strcmp(arg, 'components')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	Components = arg;
	if Components < 1
	  error('The number of components should be a value larger than 1.')
	end
      
      elseif strcmp(arg, 'biased')
	type = 'biased';
      
      elseif strcmp(arg, 'type')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'none') 
	  type = 'none';
	elseif strcmp(arg, 'biased')	  
	  type = 'biased';
	elseif strcmp(arg, 'coeff')
	  type = 'coeff';
	else
	  error(sprintf(['The argument to ''Type'' should ' ...
		'be either ''None'', ''Biased'' or' ...
		' ''coeff''. It was %s'], arg))
	end
	
      else
	error('Invalid property');
      end
      n = n + 1;
    end

    [rX, cX] = size(X);
    [rP, cP] = size(P);

    if rX ~= rP
      error('Datamatrix X and paradigm P should have same number of rows')
    end

    AllComponents = rX + rP - 1;

    if Components > AllComponents
      error(sprintf('Components = %i is too large (lyngby_xcorr)',Components));
    end

    
    nFFT = 2^nextpow2(2*rX);
    nZeroPad = nFFT - rP;

    P = [ P ; zeros(nZeroPad,1) ];
    FP = fft(P);
    
    
    C = zeros(Components, cX);
    for v = 1:cX
      if lyngby_mod(v,100) == 0
	lyngby_log(sprintf('Computing cross corr. %i / %i', v, cX));
      end;
      
      Sxp = fft([ X(:,v) ; zeros(nZeroPad,1)]) .* conj(FP);
      Cxp = real(fftshift(ifft(Sxp)));
      % 'real' used to eliminate Numerical errors 
      
      C(:,v) = Cxp((1:Components)+1+nFFT/2-ceil(Components/2),1);

    end

    if strcmp(type, 'biased')
      C = C / rX;
    elseif strcmp(type, 'coeff')
      stdP = sqrt(sum(P.^2));
      stdX = sqrt(sum(X.^2));
      C = C ./ ( ones(size(C,1),1) * (stdX * stdP + realmin) );
    end


