PCA程序设计(Matlab)

把从混合信号中求出主分量（能量最大的成份）的方法称为主分量分析（PCA），而次分量（Minor Components,MCs）与主分量（Principal Components,PCs）相对，它是混合信号中能量最小的成分，被认为是不重要的或是噪声有关的信号，把确定次分量的方法称为次分量分析（MCA）.

      PCA可以用于减少特征空间维数、确定变量的线性组合、选择最有用的变量、变量辨识、识别目标或是异常值分组等。主分量子空间提供了从高维数据到低维数据在均方误差意义下的数据压缩，它能最大程度地减少方差。

    由于PCA实际计算中只涉及到输入数据概率密度分布函数（Pdf）的二阶特性（协方差矩阵），所以解出的各主分量只互相正交（不相关），但并不满足相互独立。而且信号的大部分重要特征往往包含在Pdf的高阶统计特性中，所以只有多变量观测数据是由高斯分布的源信号构成，PCA方法才有效。

    非线性PCA（NLPCA）即将高阶累积量引入标准的PCA中，是由芬兰学者Karhunen和Oja首先提出并将其应用于ICA。它的可以完成对输入信号的盲分离。高阶累积量是以隐含的方式引入计算的，采用自适应迭代方法便于工程实现。标准的PCA基于信号的协方差矩阵仅能处理高斯信号，而NLPCA可以处理非高斯信号。

.................................................................................................................................................................................... 

程序说明：y = pca(mixedsig)，程序中mixedsig为 n*T 阶混合数据矩阵，n为信号个数，T为采样点数, y为 m*T 阶主分量矩阵。

程序设计步骤：
1、去均值
2、计算协方差矩阵及其特征值和特征向量
3、计算协方差矩阵的特征值大于阈值的个数
4、降序排列特征值
5、去掉较小的特征值
6、去掉较大的特征值(一般没有这一步）
7、合并选择的特征值
8、选择相应的特征值和特征向量
9、计算白化矩阵
10、提取主分量

 程序代码

%程序说明：y = pca(mixedsig)，程序中mixedsig为 n*T 阶混合数据矩阵，n为信号个数，T为采样点数
% y为 m*T 阶主分量矩阵。

function y = pca(mixedsig)
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

%——————————————去均值————————————
meanValue = mean(mixedsig')';
mixedsig = mixedsig - meanValue * ones(1,size(meanValue,2));
[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 = cov(mixedsig',1);    %计算协方差矩阵
[E,D] = eig(covarianceMatrix);          %计算协方差矩阵的特征值和特征向量

%———计算协方差矩阵的特征值大于阈值的个数lastEig———
rankTolerance = 1e-5;
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;
end
lowerColumns = diag(D) > lowerLimitValue;

%—————去掉较大的特征值(一般没有这一步)——————
if firstEig > 1
    higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig))/2;
else
    higherLimitValue = eigenvalues(1) + 1;
end
higherColumns = 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.');
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);

+++++++转载 Dinga's Blog　http://www.dinga.cn +++++++