mfbox_pcamat.m

toolbox for spm 5 for data, model free analysis

function [E,D,N]=mfbox_pcamat(vectors,firstEig,lastEig)
% Copyright by Peter Gruber 
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

[oldDimension,samples] = size(vectors);
if oldDimension < samples % spatial ICA I guess
    [E,D] = eig((vectors*vectors')/samples);
else % temporal ICA I guess
    [U,S,V] = svd(vectors);
    D=zeros(size(S));
    Dx= diag(S);
    [Dx,p] = sort(Dx,'ascend');
    D = S(p,p)*S(p,p)';
    E = U(:,p);
end

    
rankTolerance = 1e-7;
maxLastEig = sum(diag(D)>rankTolerance);

eigenvalues = sort(diag(D),'descend');

if ((lastEig>maxLastEig) && (maxLastEig>0)), lastEig = maxLastEig; end

if (lastEig<oldDimension), lowerLimitValue = [.5,.5]*eigenvalues(lastEig:(lastEig+1));
else lowerLimitValue = eigenvalues(oldDimension)-1;
end
lowerColumns = diag(D)>=lowerLimitValue;

if (firstEig>1), higherLimitValue = [.5,.5]*eigenvalues((firstEig-1):firstEig);
else higherLimitValue = eigenvalues(1)+1;
end
higherColumns = diag(D)<=higherLimitValue;

selectedColumns = lowerColumns&higherColumns;

D = diag(D);
N = 0;
if (~all(selectedColumns)), N = mean(D(~selectedColumns)); end
E = E(:,selectedColumns);
D = diag(D(selectedColumns));

return