pca.rar_PCA 协方差_PCA 协方差_PCA数据降维_协方差矩阵_矩阵降维

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