%function [U,R,E] = pcaFn(B);
%
%Principal components the normal way. Data in columns of B.
%U is a matrix containing the principal component eigenvectors in 
%     its columns.
%R is a matrix containing the principal component representations in its
%     rows. Each row is the representation for the corresponding column 
%     of B. 
%E is the vector of eigenvaules corresponding to the eigenvectors in U.

function [U,R,E] = pcaFn(B);

%Read data into columns of B;
%B = datamat';
[N,P] = size(B);

%********subtract mean
mb=mean(B');
B=B-(ones(P,1)*mb)';

%********Find eigenvectors vi of BB' (NxN)
[V,D] = eig (1/(P-1)*(B*B'));

%********Sort eigenvectors
eigvalvec = max(D);
[seigvals, index] = sort(eigvalvec); % sort goes low to high
Vsort = V(:,[fliplr(index)]);

U=Vsort;
length = sqrt (sum (U.^2));
U = U ./ (ones(N,1) * length);
R = B'*U;
E = fliplr(seigvals);