function [v,lambda,x2,mu] = pca(x,n,mu)
% function [v,lambda,x2,mu] = pca(x,n,mu)
%
% compute n-dimensional PCA of row data x

if ~exist('mu','var')
  mu = mean(x,2);
end

x = vadd(x,-mu);

scatter = zeros(size(x,1));
for i=1:size(x,2)
  scatter = scatter + x(:,i)*x(:,i)';
end

[v,lambda] = eig(scatter);
lambda = diag(lambda);

if ~exist('n','var')
  n = size(x,1);
end

[dummy,ind] = sort(lambda);

ind = ind(end:-1:end-n+1);

lambda = lambda(ind);
v = v(:,ind);

if nargout > 2
  x2 = v'*x;
end




function v = vadd(m,b)
% function v = vadd(m,b)
% 
% Add a vector to every column of a matrix.
%  
  
  v = m;
  for i=1:size(m,2)
    v(:,i) = v(:,i) + b;
  end
  