PCA程式設計

把從混合信号中求出主分量（能量最大的成份）的方法稱為主分量分析（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);