fastica.m

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

function [S,W]=fastICA(Z,varargin)

%  
% 函数:
%     [S,W]=fastICA(Z,varargin)
% 参数说明......
%       参数名称    ||   参数取值       || 参数说明                   
%               Z                       : 白化了的信号，每行代表一个信号
%         'method'    ['defl'|'symm']   : 采用的正交方法。
%                                         'defl'为deflationary方法
%                                         'symm'为symmetric方法，默认值为'symm'                                                mm'
%          'origW'                      : 输入初始的加权矩阵，
%                                         若省略，则由系统产生初始的正交矩阵
%      'NonlinFunc'    ['pow3'|         : 采用的非线性函数。'pow3'为g(y)=y^3
%                       'tanh'|           'tanh'为g(y)=tanh(a1*y)  
%                       'gaus']           'gaus'为g(y)=y*exp(-y^2/2)                                                            
%              'a1'                     : 非线性函数 'tanh'中调用的值,1<=a1<=2，默认值1
%              'a2'                     : 非线性函数 'gaus'中调用的值，默认值1
% 'maxNumIteration'                     : 最大的迭代次数，默认值为50
%       'OverValue'                     : 决定迭代结束的数值，其越小，表明W*W'越趋于单位阵
%                                         默认值为0.0001
%          'report'    ['on'|'off']     : 是否需要打印出程序运行信息
%                                         'on'为打印，'off'为不打印。默认为'on'
%                S                      : 恢复的信号
%                W                      : 恢复信号的加权矩阵,W=[b1,b2,...bm]
%
% 调用格式：
%     如 [S,W]=fastICA(WhitenedSig,'nonlinfunc','POW3','maxNumIteration',10) 
%     输入的可选参数必须成对。先是参数名，接着是对应的参数值。注意输入是数值还是字符串。
%




% ===================================================================
% 取默认值
ICNum=size(Z,1);    % 输入为多少个信号，则恢复出多少个信号    
method='symm';
NonlinFunc='pow3';
a1=1;
a2=1;
maxNumIteration=1000;
OverValue=0.0001;
report='on';
numSamples=size(Z,2);    % 每个信号的采样点数，即信号的长度

UserInputW=0;   % 判断是否由用户输入初始的加权矩阵。为1，则由用户输入。


% ===================================================================
% 获取可选参数
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 'method'
                method=varargin{i+1}; method=lower(method);
            case 'origw'  % origW
                W=varargin{i+1}; 
                UserInputW=1;   % 由用户输入初始的加权矩阵
            case 'nonlinfunc'   % NonlinFunc
                NonlinFunc=varargin{i+1}; NonlinFunc=lower(NonlinFunc);
            case 'a1'
                a1=varargin{i+1};  
            case 'a2'
                a2=varargin{i+1};  
            case 'maxnumiteration'   % maxNumIteration
                maxNumIteration=varargin{i+1}; 
            case 'overvalue'         % OverValue
                OverValue=varargin{i+1};  
            case 'report'
                report=varargin{i+1};  report=lower(report);
            otherwise
                % 输入参数名称有错误
                error(['Unrecognized parameter: ''' varargin{i} '''']);
        end
    end
end

if UserInputW == 0
    % 由系统产生初始的正交加权矩阵W
    W=orth(rand(ICNum)-0.5);
end

% ===================================================================
% 检测输入参数的正确性
if ~( strcmp(report,'on') | strcmp(report,'off') )
    error(sprintf('Illegal value [ %s ] for parameter: ''report''\n', report));
end

if ~isreal(Z)
    error('Input matrix has an imaginary part.\n');
end

if ~( strcmp(method,'symm') | strcmp(method,'defl'))
    error(sprintf('Illegal value [ %s ] for parameter: ''method''\n', method));
end

if ~( strcmp(NonlinFunc,'pow3') | strcmp(NonlinFunc,'tanh')| ...
        strcmp(NonlinFunc,'gaus') | strcmp(NonlinFunc,'skew') )
    error(sprintf('Illegal value [ %s ] for parameter: ''NonlinFunc''\n', NonlinFunc));
end

if ~isnumeric(a1)
    error(sprintf('Illegal value for parameter: ''a1'',it must be number\n', a1));
end

if ~isnumeric(a2) 
    error(sprintf('Illegal value for parameter: ''a2'',it must be number\n', a2));
end

if (~isnumeric(maxNumIteration) ) |  maxNumIteration < 1 
    error(sprintf('Illegal value for parameter: ''maxNumIteration''\n', maxNumIteration));
end

if ~isnumeric(OverValue) 
    error(sprintf('Illegal value for parameter: ''OverValue'', it must be number.\n', OverValue));
end

if ( size(W,1)~=size(W,2) ) | ( ~isreal(W) ) | size(W,1)~=ICNum 
    error(sprintf('Illegal value for parameter: ''origW''\n'));
end

% ===================================================================
% 打印各种信息
if strcmp(report,'on') 
    fprintf('\n=======================================\n');
    fprintf('   Information about parameters...\n');
    fprintf('=======================================\n');
    fprintf('Number of input signals: %d\n',size(Z,1));
    fprintf('Number of estimated signals: %d\n',ICNum);
    fprintf('Length of signals:%d\n',numSamples);    
    fprintf('method = %s\n',method);
    fprintf('initial W = \n');disp(W);
    fprintf('NonlinFunc = %s\n',NonlinFunc);
    fprintf('a1 = %d\n',a1);
    fprintf('a2 = %d\n',a2);
    fprintf('maxNumIteration = %d\n',maxNumIteration);
    fprintf('OverValue = %d\n',OverValue);
    fprintf('report = %s\n',report);
end

 
if strcmp(report,'on') 
    fprintf('\n=======================================\n');
    fprintf('    Starting ICA calculation ...\n'); 
    fprintf('=======================================\n');
end

% ===================================================================
% symmetric 正交方法
if strcmp(method,'symm')
    if( strcmp(report,'on') )   fprintf('Using symmetric method \n'); end
    
    WOld=zeros(size(W));  % WOld用来存贮上一次的W的值，以便和现在的W值相比较
    
    %--------------------------------------------
    % fast fixed-point ICA 算法的迭代
    for loop=1:maxNumIteration + 1
        
        % 若迭代次数到了maxNumIteration次仍没有收敛，则打印信息并返回当前的W，S值。
        if loop == maxNumIteration + 1  
            fprintf('Cannot convergence even after %d steps \n',maxNumIteration);
            S=W'*Z;
            return;
        end
        
        % Symmetric 正交化
        W=W*real(inv(W'*W)^(1/2)); 
        
        % 判断是否符合递归结束的条件。在这里，也考虑到了方向正好相反的向量
      %  minAbsCos = min(abs(diag(W'*WOld)));  
        minAbsCos = min(abs(diag(W'*WOld)));  
        
        if(strcmp(report,'on'))
            % 打印迭代过程等信息
            fprintf('-------------------\n');
            fprintf('Step No.%d, change in value of estimate: %.6f \n',loop,1-minAbsCos);           
        end
        
        if( 1-minAbsCos < OverValue)
            if( strcmp(report,'on') )  
                fprintf('\n=======================================\n');
                fprintf('  Convergence after %d steps\n',loop);
                fprintf('=======================================\n');
            end
            S=W'*Z;
            return;                %跳出整个循环－－－－－－－－
        end
        
        WOld=W;
        
        % 采用不同的非线性函数计算
        switch NonlinFunc
            % pow3
            case 'pow3'
                W=(Z*((Z'*W).^3))/numSamples-3*W;   % 在整个观测序列上求平均
            case 'tanh'
                hypTan=tanh(a1*Z'*W);
                W=Z*hypTan/numSamples -ones(size(W,1),1)*sum(1-hypTan.^2).*W/numSamples*a1;       
            case 'gaus'
                Y=Z'*W;
                Y2=Y.^2;
                ex=exp(-a2*Y2/2);
                gauss=Y.*ex;
                degauss=(1-a2*Y2).*ex;
                W=Z*gauss/numSamples - ones(size(W,1),1)*sum(degauss).*W/numSamples;
        end
        
    end   % 迭代结束
end   % symmetric method 结束


% ===================================================================
% deflation 正交方法,加权矩阵B的每个列向量一个接一个的求出来
if strcmp(method,'defl')
    if( strcmp(report,'on') )   fprintf('Using deflation method \n'); end
    
    B=zeros(ICNum);   % 最后要求出的加权矩阵
    
    % 处理初始的加权向量
    if UserInputW==0  % 由系统产生初始的加权向量
        w=rand(ICNum,1)-0.5;
    else    % 由用户输入初始的加权矩阵
        w=W(:,1);
    end
    w=w-B*B'*w;     % 正交化
    w=w/norm(w);    % 归一化
    wOld=zeros(size(w));   % wOld用来存贮上一次的w值
    
    %--------------------------------------------
    for i=1:ICNum                % 依次求出ICNum个独立分量
        
        % ------------------------------
        % fast fixed-point ICA 算法的迭代
        for loop=1:maxNumIteration+1
            % 已经找到的基向量生成了一个向量空间。现在把当前的向量投影在这个向量空间的
            % 正交空间里。从而实现正交化。注意到之所以能够利用矩阵B正确进行投影，是因
            % 为矩阵B中的零向量对投影不产生任何影响。
            w=w-B*B'*w;     % 正交化
            w=w/norm(w);    % 归一化
            
            if loop == maxNumIteration+1     % 当迭代超过maxNumIteration时，则退出返回
                fprintf('IC No.%d cannot converge in %d iteration.\n',i,maxNumIteration);
                S=B'*Z;
                return;
            end
            
            % 判断是否收敛。当收敛时，w和wOld方向相同，或者正好相反。
            % 所以，判断条件考虑了这两种情况：
            if norm(w-wOld) < OverValue | norm(w+wOld) < OverValue
                if( strcmp(report,'on') )  
                    fprintf('IC No.%d  ---------- Computed after %d iteration\n',i,loop-1);
                end
                B(:,i)=w;  % 将当前的w值保存在B中
                break;
            end
            
            wOld=w;
            
            % 采用不同的非线性函数计算
            switch NonlinFunc
                
                case 'pow3'
                    w = ( Z*( (Z'*w).^3 ) )/numSamples - 3*w;
                    
                case 'tanh'
                    hypTan = tanh(a1*Z'*w);
                    w = (Z*hypTan - a1*sum(1-hypTan.^2)'*w)/numSamples;
                    
                case 'gaus'
                    u=Z'*w;
                    u2=u.^2;
                    ex=exp(-a2*u2/2);
                    gauss=u.*ex;
                    degauss=(1-a2*u2).*ex;
                    w=(Z*gauss - sum(degauss)'*w)/numSamples;
            end
         
        end  % 对B的某个列向量w循环maxNumIteration 次 
       
    end  % 求ICNum个独立分量结束
    
    W=B;   % 加权矩阵的输出变量是W
    S=W'*Z; % 恢复的信号
    
end % deflation method 结束