fasticapearson.m

I developed an algorithm for using local ICA in denoising multidimensional data. It uses delay embed

function [y,A,W]=fasticapearson(mixedsig,epsilon, ...
                                maxNumIterations,borderBase,borderSlope, ...
				numOfIC,firstEig,lastEig,s_verbose);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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: 11.9.2000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extra parameter values. 
% tanhconst = the "frequency constant" for tanh function
% betaBase, betaSlope
%           = If , for the Pearson distribution family, 
%             alpha4<betaBase+betaSlope*alpha3^2,
%             then the distribution is estimated through
%             the generalized beta distribution (using
%             the minimum and maximum values).
% FName     = Output strings corresponding to different models
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tanhconst=1;
betaBase=2.2;
betaSlope=1.2;
FName=strvcat('Pearson', 'GBD', 'Tanh');

[Dim,numSamples]=size(mixedsig);

[mixedsig, mixedmean] = remmean(mixedsig);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Whitening data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values for 'pcamat' parameters
interactivePCA    = 'off';

[E, D]=pcamat(mixedsig, firstEig, lastEig, interactivePCA, s_verbose);
[X, whiteningMatrix, dewhiteningMatrix] = whitenv (mixedsig, E, D, s_verbose);

[Dim,numSamples]=size(X);
numOfIC=min(numOfIC,Dim);

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

switch lower(s_verbose)
 case 'on'
  b_verbose = 1;
 case 'off'
  b_verbose = 0;
 otherwise
  error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Start (main loop).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for round = 1 : maxNumIterations + 1,
  if round == maxNumIterations + 1,
    if (b_verbose>0) warning(sprintf('No convergence after %d steps', maxNumIterations)); end
    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);

  if (b_verbose>0) warning(sprintf('Step no. %d; estimated scores:',round)); end 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 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;
    elseif 0
%    elseif (alpha4(i)<betaBase+betaSlope*alpha3(i)^2 && (alpha4(i)>1+alpha3(i)^2) && (alpha4(i)<3+2*alpha3(i)^2))
        FType(i)=2; % GBD
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % GBD sample min && max estimation used. 
	% For details, see gbd_momentfit.m  and
	%   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
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        samplemin=min(U);
        samplemax=max(U);
        beta=gbd_momentfit(alpha3(i),alpha4(i),...
	                 samplemin(i),samplemax(i),numSamples);
        [score1(i,:) score2(i,:)]=gbd_score(U(:,i)',beta);
    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   
    if (b_verbose>0) warning(sprintf('[%d]', i)); end 
  end
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate nonlinearities (scores) for each signal (end of the loop).
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  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.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if (all(isfinite(B(:))) == 0)
      B = randn(size(B));
      B = B * real((B' * B)^(-1/2));
  end
  B = B * real((B' * B)^(-1/2));
  while (all(isfinite(B(:))) == 0)
      B = randn(size(B));
      B = B * real((B' * B)^(-1/2));
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Show the progress...
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  minAbsCos = min(abs(diag(B'*BOld)));
  if (b_verbose>0)
   warning(sprintf('Step no. %d, change in value of estimate: %.6f \n', round, 1 - minAbsCos)); 
   warning(EstimationValues); 
  end

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % 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)),
    if (b_verbose>0) warning(sprintf('Convergence after %d steps', round)); end 
    break;
  end  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the results.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
W = B'*whiteningMatrix; 
A = dewhiteningMatrix * B;
y = W*mixedsig;

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

return