pca程序分析

1 本周工作

    学习球磨机数据中的PCA程序。

2 学习总结

      <1>程序设计步骤：

1、去均值

2、计算协方差矩阵及其特征值和特征向量

3、计算协方差矩阵的特征值大于阈值的个数

4、降序排列特征值

5、去掉较小的特征值

6、去掉较大的特征值(一般没有这一步）

7、合并选择的特征值

8、选择相应的特征值和特征向量

9、计算白化矩阵

10、提取主分量

<2>部分程序分析

程序说明：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;

3参考文献

      郭磊师兄的程序

      网络博客

      百度文库