www.pudn.com > ICA.rar > PCA.m, change:2003-11-16,size:5919b


function [NewSig,E,D] = PCA(OldSig,varargin) 
 
% PCA ---- Performing PCA algorithm to the input signals 
% 调用格式: 
%     [NewSig,E,D] = PCA(OldSig,'PCNum',4,'minIndex',1,'maxIndex',4,'report','on') 
%     [NewSig,E,D] = PCA(OldSig,'PCNum',10,'report','off') 
% 
% 参数说明: 
%   参数名称   ||   参数取值   ||  参数说明  
%     'PCNum'      (整数)        要求出的主分量数,如果用户输入该参数,则minIndex和maxIndex无效 
%  'minIndex'      (整数)        选择的最小的特征值的序号(在EIG分解后的特征值矩阵中的排列序号) 
%  'maxIndex'      (整数)        选择的最大的特征值的序号(在EIG分解后的特征值矩阵中的排列序号)   
%    'report'      [ 'on' |       'on'表示显示运行的相关信息 
%                    'off']       'off'表示不显示信息 
%     OldSig                       输入的信号矩阵。每行代表一个信号;每列代表各个信号的一次观测 
%     NewSig                       输出的经过PCA处理后的信号矩阵 
%          E                       PCA处理后的特征向量矩阵 
%          D                       PCA处理后的特征值矩阵 
% 
% 注意: 
%     可选参数必须成对输入。参数名在前,对应的参数值紧跟在其后。注意数值和字符串的区别。 
% 
% 参考: 
%     SELCOL  FFPICA 
% 
% 作者:张智林(Zhang Zhi-Lin), 
%       现代信号处理实验室,zhangzl@vip.163.com 
% 版本:1.1 
% 日期:2003年11月16日 
 
 
 
OldDim=size(OldSig,1); 
%====================================================================== 
% 默认值设置 
PCNum=OldDim; 
minIndex=1; 
maxIndex=OldDim; 
report='on'; 
 
enterPCNum=0;      % 标志用户是否输入PCNum的值。为1,则表示用户已输入PCNum的值, 
                   % 此时用户输入的minIndex和maxIndex均无效,程序根据用户输入的 
                   % PCNum值来进行PCA计算 
 
                    
%====================================================================== 
% 获取可选参数 
if(mod(length(varargin),2)==1) 
    error('Optional parameters should always go by pairs.\n'); 
else 
    for i=1:2:(length(varargin)-1) 
        switch lower(varargin{i}) 
            case 'pcnum' 
                PCNum=varargin{i+1}; 
                enterPCNum=1;  % 用户输入了PCNum的值,此时用户输入的minIndex和maxIndex值无效 
            case 'minindex' 
                minIndex=varargin{i+1}; 
            case 'maxindex' 
                maxIndex=varargin{i+1}; 
            case 'report' 
                report=lower(varargin{i+1}); 
            otherwise 
                error(['Unrecognized parameter: ''' varargin{i} '''']); 
        end 
    end 
end 
 
 
%==================================================================== 
% 检测输入参数的正确性 
if ~( strcmp(report,'on') | strcmp(report,'off') ) 
    error(sprintf('Illegal value [ %s ] for parameter: ''report''\n', report)); 
end 
 
if ~isreal(OldSig) 
    error('Input matrix has an imaginary part.\n'); 
end 
 
if ~isnumeric(minIndex)  
    error(sprintf('Illegal value for parameter: ''minIndex'', it must be number\n', minIndex)); 
elseif minIndex > OldDim | minIndex < 1 
    error(sprintf('Illegal value for parameter: ''minIndex'', input correct value\n', minIndex)); 
end 
 
if ~isnumeric(maxIndex)  
    error(sprintf('Illegal value for parameter: ''maxIndex'', it must be number\n', maxIndex)); 
elseif maxIndex > OldDim | maxIndex < minIndex  
    error(sprintf('Illegal value for parameter: ''maxIndex'', input correct value\n', maxIndex));     
end 
 
if ~isnumeric(PCNum) | PCNum > OldDim 
    error(sprintf('Illegal value for parameter: ''PCNum'', it must be number or correct value\n', PCNum)); 
end 
 
 
%------------------------------------------------------------------------- 
% 程序根据用户输入的PCNum值进行计算 
if enterPCNum ==1         
    maxIndex = OldDim; 
    minIndex = OldDim-PCNum+1; 
end 
 
 
%========================================================================= 
% 计算PCA 
if strcmp(report,'on') fprintf('Calculating PCA ...\n'); end 
 
covarianceMatrix=cov(OldSig',1);          % 计算协相关阵 
nonZeroEig=rank(covarianceMatrix,1e-9);   % 计算协相关阵的秩,也即独立的列数或行数 
[E,D]=eig(covarianceMatrix);              % 计算特征向量矩阵和特征值矩阵,注意特征值由小到大排列 
nonZeroIndex = OldDim - nonZeroEig + 1;   % 非零的最小特征值的序号 
 
 
%-------------------------------------------------------------------------- 
% 检测是否降维 
if minIndex < nonZeroIndex 
    minIndex = nonZeroIndex;        % 使最小特征值的序号 等于 最小的非零特征值的序号 
    if strcmp(report,'on')  
        fprintf('Dimension reduced from %d to %d because of the singularity of covaricance matrix\n', ...  
                 OldDim,maxIndex - minIndex +1);  
    end 
else   
    if strcmp(report,'on')  
        if OldDim == (maxIndex - minIndex + 1) 
            fprintf('Dimension has not reduced.\n'); 
        else 
            fprintf('Dimension reduced from %d to %d\n',OldDim,maxIndex - minIndex +1); 
        end 
    end 
end 
  
 
%-------------------------------------------------------------------------- 
% 选取特征值 
selectedColumns = zeros(OldDim,1); 
 
for i=1:OldDim    % 列向量selectedColumns中的元素1表示选取对应序号的特征值和特征向量 
    if (minIndex <= i) & (i <= maxIndex) 
        selectedColumns(i)=1; 
    end 
end 
 
sumTotal = sum(diag(D));                            % 原输入信号的特征值总和 
sumRemove = sum( diag(D).*(~selectedColumns) );     % 舍去的特征值之和 
 
if strcmp(report,'on') 
    fprintf('Smallest remaining(non-zero) eigenvalue: %d\n',D(minIndex,minIndex)); 
    fprintf('Largest remaining(non-zero) eigenvalue: %d\n',D(maxIndex,maxIndex)); 
    fprintf('Sum of removed eigenvalue: %d\n',sumRemove ); 
end 
 
% 真正计算选取的特征值矩阵和特征向量矩阵 
E=selcol( E, selectedColumns ); 
D=selcol( selcol(D,selectedColumns)', selectedColumns); 
 
sumUsed = sum(diag(D)); 
 
if strcmp(report,'on') 
    fprintf('%g%% of (non-zero)eigenvalues retained.\n',sumUsed/sumTotal*100); 
end 
 
NewSig = E'*OldSig;      % 降维,取主分量后得到的新的信号矩阵