pearson_ica.m

PearsonICA程序

function [y, A, W]=pearson_ica(mixedsig, s1, v1, s2, v2,...
			       s3, v3, s4, v4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pearson-ICA algorithm for Independent Component Analysis
% 
% Input:
%   mixedsig = n*m matrix of the observed mixture
%   s1       = string corresponding for a parameter name (optional)
%   v1       = value of the paramater given in s1        (optional)
%   s2       = string corresponding for a parameter name (optional)
%   v2       = value of the paramater given in s2        (optional)
%   s3       = string corresponding for a parameter name (optional)
%   v3       = value of the paramater given in s3        (optional)
%   s4       = string corresponding for a parameter name (optional)
%   v4       = value of the paramater given in s4        (optional)
%
% Output:
%   y        = n*m matrix of the separated signals
%   A        = n*n matrix (estimated mixing matrix)
%   W        = n*n matrix (estimated demixing matrix)
%   
% pearson_ica(mixedsig) estimates the independent components from given
% multidimensional signals. Each row of the matrix mixedsig is an
% observed signal.  Pearson_ICA employs maximum likelihood approach where 
% the distributions of source signals are iteratively estimated using
% the Pearson system, see
%
% Karvanen J. and Koivunen V.:
% "Blind separation methods based on Pearson system and its extensions",
% Signal Processing, 2002, to appear
%
% Karvanen, J.,Eriksson, J., and Koivunen, V.:
% "Pearson System Based Method for Blind Separation", 
% Proceedings of Second International Workshop on
% Independent Component Analysis and Blind Signal Separation,
% Helsinki 2000, pp. 585--590
%
% The optimization is done by Hyvarinen's fixed-point algorithm,
% see http://www.cis.hut.fi/projects/ica/fastica/.
%
% Optional parameters:
%   'epsilon'           Convergence criterion
%   'maxNumIterations'  Maximum number of iterations
%   'borderBase'
%   'borderSlope'       The minimum kurtosis value for switching to 
%                       the GBD (tanh) area. Thus the algorithm
%                       uses GLD for values
%                       kurtosis>borderBase+borderSlope*skewness^2, 
%                       and GBD (tanh) for the other values. 
%                       The default slope value 1.75 comes from the
%                       relation between theoretical minimum kurtosis and
%                       the corresponding skewness. 
%                       
% Example: 
%  ICs=pearson_ica(mixedsig,'epsilon',0.00005,...
%                          'maxNumIterations',50,'borderBase',[2.5 4]);
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Juha Karvanen, Jan Eriksson, and Visa Koivunen.
%
% For details see the files readme.txt
% and gpl.txt provided within the same package.
%
% Last modified: 11.9.2001  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values for parameters. This corresponds to the standard
% Pearson-ICA algorithm such that the Pearson family is used if
% 2.6<=kurtosis<=4+skewness^2, and the contrast tanh is used otherwise.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
epsilon=0.0001;
maxNumIterations=200;
borderBase=[2.6 4];
borderSlope=[0 1];

if(rem(nargin-1,2)==1)
  error('Optional parameters should always go by pairs');
else
  for i=1:(nargin-1)/2
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Extract the name and the value of parameter number i (0<=i<=4).
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    str_param = eval (['s' int2str(i)]);
    val_param = eval (['v' int2str(i)]);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Set the value of the parameter number i.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if strcmp (str_param, 'epsilon')
      epsilon = val_param;
    elseif strcmp (str_param, 'maxNumIterations')
      maxNumIterations = val_param;
    elseif strcmp (str_param, 'borderBase')
      if size(val_param(:))~=2,
	error('You should specify two borders (up&low) for Pearson system!');
      else
        borderBase = val_param;
	end
    elseif strcmp (str_param, 'borderSlope')
      if size(val_param(:))~=2,
	error('You should specify two borders (up&low) for Pearson system!');
      else
        borderSlope = val_param;
      end
    end
  end 
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Call for the actual algorithm doing the dirty job.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[y, A, W]=fasticapearson(mixedsig,epsilon,...
		      maxNumIterations,borderBase,borderSlope);

% The end of the function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%