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 = 0.5;
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);