fasticapearson.m

PearsonICA程序

function [y,A,W]=fasticapearson(mixedsig,epsilon,...
			     maxNumIterations,borderBase,borderSlope);
           
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pearson-ICA algorithm for Independent Component Analysis
%
% This is the core program for Pearson-ICA. 
% Don't use this program directly, but 
% use the program pearson_ica that call this code.
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
%
% Input: 
%   mixedsig         = n*m matrix of the observed mixture 
%   epsilon          = scalar (convergence constant)
%   maxNumIterations = integer of maximum number of iterations
%   borderBase, 
%   borderSlope      = constants for the line for switching between
%                      Pearson and tanh are:
%                      Pearson is used iff
%                      kurtosis>borderBase(1)+borderSlope(1)*skewness^2 or
%                      kurtosis<borderBase(2)+borderSlope(2)*skewness^2. 
% Output:
%   y                = n*m matrix of the separated sources (rows)
%   A                = estimated n*n mixing matrix
%   W                = estimated n*n separation matrix
%
% The function estimates the independent components from given 
% multidimensional signals. Each row of the matrix mixedsig is 
% an observed signal. For optimization Hyvarinen's fixed-point algorithm,
% see http://www.cis.hut.fi/projects/ica/fastica/, is used. 
% Some of the code in this function is based on the FastICA
% package distributed under the same lisence.           
%
% example: [y, A, W]=fasticapearson(mixedsig,0.0001,100,[2.6 4],[0 1])
%
% 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: 9.10.2001
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extra parameter values. 
% tanhconst = the "frequency constant" for tanh function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tanhconst=1;
FName=strvcat('Pearson', 'GBD', 'Tanh');

[numOfIC,numSamples]=size(mixedsig);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
alpha3=zeros(1,numOfIC);
alpha4=3*ones(1,numOfIC);
FTypeOld=zeros(1,numOfIC);
EstimationValues=[];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Random starting point (rotation matrix).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B = orth(rand(numOfIC) - .5);
BOld = B;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Whitening data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zeromeansig=mixedsig-mean(mixedsig')'*ones(1,size(mixedsig,2));
whiteningMatrix=inv(real(sqrtm(cov(zeromeansig'))));
X=whiteningMatrix*zeromeansig;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start (main loop).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for round = 1 : maxNumIterations + 1,
  if round == maxNumIterations + 1,
    fprintf('No convergence after %d steps\n', maxNumIterations);
    break;
  end
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate the independent component estimates (u_i's),
  % u_i = b_i' x = x' b_i. 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  U = X' * B;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Estimate 3rd and 4th moments.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  alpha3=mean(U.^3);
  alpha4Old=alpha4;
  alpha4=mean(U.^4);

  fprintf('Step no. %d; estimated scores:',round); 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate nonlinearities (scores) for each signal (beginning of the loop).
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for i=1:size(U,2),
    if (alpha4(i)<borderBase(1)+borderSlope(1)*alpha3(i)^2) ...
	  | (alpha4(i)>borderBase(2)+borderSlope(2)*alpha3(i)^2),
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      % tanh used. See FastICA (http://www.cis.hut.fi/projects/ica/fastica/)
      % and related papers for more details.
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      FType(i)=3; %tanh
      score1(i,:)=tanh(Tanhconst * U(:,i)');
      score2(i,:)=1-score1(i,:).^2;
    else
        FType(i)=1; %Pearson       
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Pearson moment estimation used. 
	     % See pearson_momentfit.m for more details.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        [pearsona pearsonb]=pearson_momentfit(alpha3(i),alpha4(i));
        [score1(i,:) score2(i,:)]=pearson_score(U(:,i)',pearsonb,pearsona,20);
       end   
  %fprintf('[%d]', i); 
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate nonlinearities (scores) for each signal (end of the loop).
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fprintf('\n'); 

  EstimationValues=strcat(strjust(strvcat('Signal no.',...
					num2str([1:numOfIC]'))),...
                 strjust(strvcat('3rd moment: ', num2str(alpha3'))),...
                 strjust(strvcat('4th moment: ', num2str(alpha4'))),...
		 strjust(strvcat('Model used: ',FName(FType,:))));
  
  B=X*score1'/numSamples-...
    ones(size(B,1),1)*sum(ones(size(U)).*score2').*B/numSamples ;
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Symmetric orthogonalization.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  B = B*real(sqrtm(inv(B'*B))); 

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Show the progress...
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  minAbsCos = min(abs(diag(B'*BOld)));
  fprintf('Step no. %d, change in value of estimate: %.6f \n',...
	     round, 1 - minAbsCos);
  disp(EstimationValues);
  disp(blanks(2)');
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Test for termination condition. Note that we consider opposite
  % directions here as well. Program is not terminated if any of
  % the model distributions is changed in the last iteration.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ((1 - minAbsCos < epsilon)  & (max(abs(FType-FTypeOld))==0)),
    fprintf('Convergence after %d steps\n', round); 
    break;
  end  

  BOld = B;
  FTypeOld=FType;  
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End (main loop).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = B'*whiteningMatrix; 
A=inv(W);
y=W*mixedsig;

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