function F = lyngby_baf_test(X, R, arg1, arg2, arg3, arg4, arg5)

% lyngby_baf_test      - Barbak Ardekani's F-test
%
%       function F = lyngby_ba_ftest(X, R, 'PropertyName',
%           'PropertyValue')
%
%       Input:  X   The datamatrix
%               R   The run structure
%               Property:
%                 'Components'  { 8 } Number of Fourier components
%                               used in the test
%
%       Output: F   F-map
%
%       See also LYNGBY, LYNGBY_BAT_TEST.

% cvs : $Id: lyngby_baf_test.m,v 1.8 2000/09/27 13:58:46 fnielsen Exp $
%       $Revision: 1.8 $

    % Sizes
    [NoTime,NoVoxels] = size(X);
    
    % Default values
    Components = 8;
    
    % Parsing
    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 < 0 
	  error('The number of components should be non-negative')
	end
      else
	error('Invalid property');
      end
      n = n + 1;
    end

    M = Components;

    % Find number of scans
    [a,T] = lyngby_runinfo(R);
    
    if length(T) > 1
       error(['The Ardekani F-test cannot handle runs with ', ...
	     'different lengths.']);
    end
    msplit = floor(M/2);
    omega = 2*pi/T;
    in = (0:NoTime-1)';

    for k = 1:msplit,
      z(:,k) = cos(k*omega*in);
    end

    for k = msplit+1:M,
      q = k-msplit;
      z(:,k) = sin(q*omega*in);
    end  

    Pz = z*inv(z'*z)*z';
    Pzr = Pz/M;
    dPzr = (eye(size(Pz))-Pz)/(NoTime-M-1);

    F = zeros(1,NoVoxels);

    for v = 1:NoVoxels,
      if lyngby_mod(v,20) == 0,
	lyngby_log(sprintf('Computing F-map: %i / %i', v, ...
	    NoVoxels));
      end;

      x = X(:,v);
      if (std(x)<1e-9)
        F(1,v) = 0;
      else
        F(1,v) = (x'*Pzr*x)/(x'*dPzr*x);
      end
    end
    lyngby_log(sprintf('Done computing F-map'));

















