PCA实现代码

        整了几天，终于把pca的matlab代码完成了。你妹哦，网上的代码太坑人了！！

        

function [y whiteningMatrix] = pcaRecogFace(mixedsig,rankTolerance)%程序说明：y = pca(mixedsig)，程序中mixedsig为 n*T 阶混合数据矩阵，n为信号个数，T为采样点数% y为 m*T 阶主分量矩阵。% n是维数，T是样本数。%%rankTolerance表示最小特征值。%——————————————去均值————————————meanValue = mean(mixedsig')';[m,n] = size(mixedsig);for s =  1:m    for t = 1:n        mixedsig(s,t) = mixedsig(s,t) - meanValue(s);    endend[Dim,NumofSampl] = size(mixedsig);oldDimension = Dim;fprintf('Number of signals: %d\n',Dim);fprintf('Number of samples: %d\n',NumofSampl);fprintf('Calculate PCA...');firstEig = 1;lastEig = Dim;covarianceMatrix = corrcoef(mixedsig');    %计算协方差矩阵[E,D] = eig(covarianceMatrix);          %计算协方差矩阵的特征值和特征向量%———计算协方差矩阵的特征值大于阈值的个数lastEig———maxLastEig = sum( diag(D) > rankTolerance );lastEig = maxLastEig;%——————————降序排列特征值——————————eigenvalues = flipud(sort(diag(D)));%—————————去掉较小的特征值——————————if lastEig < oldDimension    lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1))/2;else    lowerLimitValue = eigenvalues(oldDimension) - 1;endlowerColumns = diag(D) > lowerLimitValue;%—————去掉较大的特征值(一般没有这一步)——————if firstEig > 1    higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig))/2;else    higherLimitValue = eigenvalues(1) + 1;endhigherColumns = diag(D) < higherLimitValue;%—————————合并选择的特征值——————————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 eigenvalue[ %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)if size(maskVector,1)~= size(oldMatrix,2)    error('The mask vector and matrix are of uncompatible size.');endnumTaken = 0;for i = 1:size(maskVector,1)    if maskVector(i,1) == 1          takingMask(1,numTaken + 1) = i;        numTaken = numTaken + 1;    endendnewMatrix = oldMatrix(:,takingMask);

           将上述代码在贝母中，按第一主成分、第二主成分和第三主成分分别形成图像，效果如下：

                                  

          可见，通过不同主成分投影，结果是使得数据更有差异。在代码中为何采用白化矩阵进行投影还有待后续研究。