www.pudn.com > zxb.rar > pca.m

function y=pca(mixedsig) 
%x is nxT matrix 
%n is the number of signals 
%T is the number of samples 
%2005.12.1 
 
if nargin == 0, 
  error ('You must supply the mixed data as input argument.'); 
end 
 
if length (size (mixedsig)) > 2, 
  error ('Input data can not have more than two dimensions.'); 
end 
 
if any (any (isnan (mixedsig))), 
  error ('Input data contains NaN''s.'); 
end 
 
% Remove the mean and check the data 
 
meanValue = mean (mixedsig')'; 
mixedsig = mixedsig - meanValue * ones (1,size (mixedsig, 2)); 
 
 
[Dim, NumOfSampl] = size(mixedsig); 
oldDimension = Dim; 
fprintf('Number of signals: %d\n', Dim); 
fprintf('Number of samples: %d\n', NumOfSampl); 
fprintf('Calculate PCA ....\n'); 
firstEig          = 1; 
lastEig           = Dim; 
%[E, D]=pcamat(mixedsig, firstEig, lastEig, interactivePCA, verbose); 
% Calculate the covariance matrix. 
covarianceMatrix = cov(mixedsig', 1); 
% Calculate the eigenvalues and eigenvectors of covariance matrix. 
[E, D] = eig (covarianceMatrix); 
%calculate the number of eigenvalue that bigger than tolerance 
rankTolerance = 1e-7; 
maxLastEig = sum (diag (D) > rankTolerance); 
lastEig = maxLastEig; 
% Sort the eigenvalues - decending. 
eigenvalues = flipud(sort(diag(D))); 
 
% Drop the smaller eigenvalues 
if lastEig < oldDimension 
  lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1)) / 2; 
else 
  lowerLimitValue = eigenvalues(oldDimension) - 1; 
end 
 
lowerColumns = diag(D) > lowerLimitValue; 
 
% Drop the larger eigenvalues 
if firstEig > 1 
  higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig)) / 2; 
else 
  higherLimitValue = eigenvalues(1) + 1; 
end 
higherColumns = diag(D) < higherLimitValue; 
 
% Combine the results from above 
selectedColumns = lowerColumns & higherColumns; 
 
fprintf ('Selected [ %d ] dimensions.\n', sum (selectedColumns)); 
fprintf ('Smallest remaining (non-zero) eigenvalue [ %g ]\n', eigenvalues(lastEig)); 
fprintf ('Largest remaining (non-zero) eigenvalue [ %g ]\n', eigenvalues(firstEig)); 
fprintf ('Sum of removed eigenvalues [ %g ]\n', sum(diag(D) .* (~selectedColumns))); 
 
E = selcol(E, selectedColumns); 
D = selcol(selcol(D, selectedColumns)', selectedColumns); 
 
 
whiteningMatrix = inv (sqrt (D)) * E'; 
dewhiteningMatrix = E * sqrt (D); 
y =  whiteningMatrix * mixedsig; 
 
function newMatrix = selcol(oldMatrix, maskVector); 
 
% newMatrix = selcol(oldMatrix, maskVector); 
% 
% Selects the columns of the matrix that marked by one in the given vector. 
% The maskVector is a column vector. 
 
% 15.3.1998 
 
if size(maskVector, 1) ~= size(oldMatrix, 2), 
  error ('The mask vector and matrix are of uncompatible size.'); 
end 
 
numTaken = 0; 
 
for i = 1 : size (maskVector, 1), 
  if maskVector(i, 1) == 1, 
    takingMask(1, numTaken + 1) = i; 
    numTaken = numTaken + 1; 
  end 
end 
 
newMatrix = oldMatrix(:, takingMask);