PCA源代码 基于MATLAB的


01.function y = pca(mixedsig)

02.

03.%程序说明:y = pca(mixedsig),程序中mixedsig为 n*T 阶混合数据矩阵,n为信号个数,T为采样点数

04.% y为 m*T 阶主分量矩阵。

05.% n是维数,T是样本数。

06.

07.if nargin == 0

08.    error('You must supply the mixed data as input argument.');

09.end

10.if length(size(mixedsig))>2

11.    error('Input data can not have more than two dimensions. ');

12.end

13.if any(any(isnan(mixedsig)))

14.    error('Input data contains NaN''s.');

15.end

16.

17.%——————————————去均值————————————

18.meanValue = mean(mixedsig')';

19.[m,n] = size(mixedsig);

20.%mixedsig = mixedsig - meanValue*ones(1,size(meanValue)); %当数据本身维数很大时容易出现Out of memory

21.for s =  1:m

22.    for t = 1:n

23.        mixedsig(s,t) = mixedsig(s,t) - meanValue(s);

24.    end

25.end

26.[Dim,NumofSampl] = size(mixedsig);

27.oldDimension = Dim;

28.fprintf('Number of signals: %d\n',Dim);

29.fprintf('Number of samples: %d\n',NumofSampl);

30.fprintf('Calculate PCA...');

31.firstEig = 1;

32.lastEig = Dim;

33.covarianceMatrix = corrcoef(mixedsig');    %计算协方差矩阵

34.[E,D] = eig(covarianceMatrix);          %计算协方差矩阵的特征值和特征向量

35.

36.%———计算协方差矩阵的特征值大于阈值的个数lastEig———

37.%rankTolerance = 1;

38.%maxLastEig = sum(diag(D) >= rankTolerance);

39.%lastEig = maxLastEig;

40.lastEig = 10;

41.

42.%——————————降序排列特征值——————————

43.eigenvalues = flipud(sort(diag(D)));

44.

45.%—————————去掉较小的特征值——————————

46.if lastEig < oldDimension

47.    lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1))/2;

48.else

49.    lowerLimitValue = eigenvalues(oldDimension) - 1;

50.end

51.lowerColumns = diag(D) > lowerLimitValue;

52.

53.%—————去掉较大的特征值(一般没有这一步)——————

54.if firstEig > 1

55.    higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig))/2;

56.else

57.    higherLimitValue = eigenvalues(1) + 1;

58.end

59.higherColumns = diag(D) < higherLimitValue;

60.

61.%—————————合并选择的特征值——————————

62.selectedColumns =lowerColumns & higherColumns;

63.

64.%—————————输出处理的结果信息—————————

65.fprintf('Selected [%d] dimensions.\n',sum(selectedColumns));

66.fprintf('Smallest remaining (non-zero) eigenvalue[ %g ]\n',eigenvalues(lastEig));

67.fprintf('Largest remaining (non-zero) eigenvalue[ %g ]\n',eigenvalues(firstEig));

68.fprintf('Sum of removed eigenvalue[ %g ]\n',sum(diag(D) .* (~selectedColumns)));

69.

70.%———————选择相应的特征值和特征向量———————

71.E = selcol(E,selectedColumns);

72.D = selcol(selcol(D,selectedColumns)',selectedColumns);

73.

74.%——————————计算白化矩阵———————————

75.whiteningMatrix = inv(sqrt(D)) * E';

76.dewhiteningMatrix = E * sqrt(D);

77.

78.%——————————提取主分量————————————

79.y = whiteningMatrix * mixedsig;

80.

81.%——————————行选择子程序———————————

82.function newMatrix = selcol(oldMatrix,maskVector)

83.if size(maskVector,1)~= size(oldMatrix,2)

84.    error('The mask vector and matrix are of uncompatible size.');

85.end

86.numTaken = 0;

87.for i = 1:size(maskVector,1)

88.    if maskVector(i,1) == 1 

89.        takingMask(1,numTaken + 1) = i;

90.        numTaken = numTaken + 1;

91.    end

92.end

93.newMatrix = oldMatrix(:,takingMask);
复制代码

展开阅读全文

没有更多推荐了,返回首页