% PCA.m

% It is assumed that data is already normalized to mean zero and variance one!!!!
% the input proj may be empty, ie., []

function [PCAdata,proj] = PCA(data,proj,odim)

% get size of pattern matrix. # patterns and # features (dimensionality)
[npat,dim] = size(data)
if (dim < odim)
   error('Desired output dimensionality > input dimensionality');
end

if isempty(proj)		% it means the data is training data
   % compute covariance matrix and perform eigenanalysis
   c = cov(data);
   [vecs,vals] = eig(c); % note vals is a dim x dim diagonal mtx
   % calculate total variance in this data = sum of eigenvalues
   s2 = trace(vals);
   % sort eigenvalues. index is the array of permutation indices
   [sv,index]=sort(diag(vals));
   % form projection matrix.  First coord has highest variance, etc.
   proj = fliplr(vecs(:,index(dim-odim+1:dim)));    % yecccch
   % compute % of variance retained;
   s3 = sum(sv(dim-odim+1:dim));
   fprintf(1,'%f%% of variance retained\n',100.0*s3/s2);
end
% in case the input projection matrix has more dimentions than desired
proj=proj(1:dim,1:odim);

% project matrix
PCAdata = data * proj;

return;