function [psi, a, evals] = pca(f, N, type)

% [psi, a, evals]=pca(f, N, type)
%
% performs KL analysis on snaps 
% Inputs:
% data a T x M matrix (smallest dimension first)
% N is number of eigenfunctions to calculate
% type is 'float32' [default], 'int16' or other
%
% returns:
%   psi: matrix whose columns are spatial eigenvectors
%   a:   matrix whose columns are temporal eigenvectors
%   evals:  vector of eigenvalues in descending order

% if nargin<5
%   type='float32';
% end

P=size(f,2);

cor=f*f';           % form correlation matrix (doesn't exploit symmetry)
[a,lam]=eig(cor);                % find temporal eigenfunctions/values
%a = eigenvectors, lam = diagonal matrix of eigenvalues
[evals,ind]=sort(diag(lam));
evals=flipud(evals); 
lam=diag(evals);
isig=diag(evals.^(-.5));
a=a(:,flipud(ind));
a=a(:,1:N);
psi=f'*a*isig(1:N,1:N);          % form  1st N spatial eigenfns
fic=a*diag(evals(1:N).^(0.5));
