pca.m

完整的独立分量分析(ICA)程序包,可以用来进行盲源分离

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），
%       现代信号处理实验室,[e-mail]
% 版本：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;      % 降维，取主分量后得到的新的信号矩阵