function [Xtrunc, Vtrunc, Y] = PCAfun (data, PC)
% inputs load data and establish number of components to use

 [m n] = size(data); 	% obtain dimensions of data matrix
 mu = mean(data,2); 	% compute row mean
 X = data-repmat(mu,1,n); 	% subtract row mean to obtain X
 Z = 1/sqrt(n-1)*X'; 	% create matrix, Z
 covar = Z'*Z; 	% covariance matrix of Z
 numeigs = min(size(covar)); % number of eigenvalues to plot
 
% Extract eigenvalues and eigenvectors
 [V, Lambda] = eigs(covar, numeigs,'la'); % compute eigenvectors and eigenvalues
 evals = diag(Lambda); % store eigenvalues in vector form
 fig2 = figure('Color',[1 1 1]); 
 bar(evals);	% bar plot of eigenvalues (variances)
 
 % compute percentage of variance included
 totvals = sum(evals);	%total sum of eigenvalues
 runtot = cumsum(evals);	%cumulative sum of eigenvalues - largest to smallest
 percentinc = runtot(PC)/totvals 	% percentage included

% Extract principal components
 Vtrunc = V(:,1:PC); % extract first PC vectors in new basis
 Y = Vtrunc' * X; % project data onto new basis
 Xtrunc = Vtrunc*Y; 	% convert back to original basis
 Xtrunc = Xtrunc+repmat(mu,1,n); 	% add the row means back
 
