fasticaegld.m

这是盲信号的代码 都已经通过编译了 做这方面的同仁可以参考一下 我觉得蛮惯用的

function [y,A,W]=fasticaegld(mixedsig,sub_alg,epsilon,...
			     maxNumIterations,borderBase,borderSlope);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EGLD-ICA algorithm for Independent Component Analysis
%
% This is the core program for EGLD-ICA. 
% Don't use this program directly, but 
% use the program egld_ica that call this code.
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%
% Input: 
%   mixedsig         = n*m matrix of the observed mixture 
%   sub_alg          = string for the variation of EGLD-ICA
%                      Possible values are 'egld', 'gld', and 'gldtanh'.
%   epsilon          = scalar (convergence constant)
%   maxNumIterations = integer of maximum number of iterations
%   borderBase, 
%   borderSlope      = constants for the line for switching between models:
%                      kurtosis>borderBase+borderSlope*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]=fasticaegld(mixedsig,'egld',0.0001,100,2.2,1.8)
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Jan Eriksson, Juha Karvanen, and Visa Koivunen.
%
% For details, see the files readme.txt
% and gpl.txt provided within the same package.
%
% Last modified: 8.9.2000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Extra parameter values. 
% tanhconst = the "frequency constant" for tanh function
% FName     = Output strings corresponding for different models
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Tanhconst=1;
FName=strvcat('GLD', 'GBD', 'Tanh');

[numOfIC,numSamples]=size(mixedsig);

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Method to be used.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if strcmp(lower(sub_alg),'egld'),
  % Use GLD and GBD. 
  egldtype=1;
elseif strcmp(lower(sub_alg),'gldtanh'),
  % Use GLD and tanh.
  egldtype=2;
elseif strcmp(lower(sub_alg),'gld'),
  % Use GLD only.
  egldtype=3;
else
  error('Code for the desired contrast function not found!');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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.
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  alfa3=mean(U.^3);
  alfa4Old=alfa4;
  alfa4=mean(U.^4);

  fprintf('Step no. %d; estimated scores:',round); 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Calculate nonlinearities (scores) for each signal (beginning of the loop).
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  for signalNum=1:numOfIC,
   if egldtype==1,
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
     % EGLD (GLD+GBD) used.
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     if (alfa4(signalNum)>borderBase+borderSlope*alfa3(signalNum)^2),
       FType(signalNum)=1; %GLD
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % GLD moment estimation used. See gld_momentfit.m for more details.
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       lambda=gld_momentfit(alfa3(signalNum),alfa4(signalNum));
       [score1(signalNum,:) score2(signalNum,:)]=...
	                 gld_score(U(:,signalNum)',lambda);
     else
       FType(signalNum)=2; %GBD
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       % Minimum & maximum based GBD moment estimation used. See
       % gbd_momentfit.m for more details.
       %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       samplemin=min(U);
       samplemax=max(U);
       beta=gbd_momentfit(alfa3(signalNum),alfa4(signalNum),...
                         samplemin(signalNum),samplemax(signalNum),numSamples);
       [score1(signalNum,:) score2(signalNum,:)]=...
	                 gbd_score(U(:,signalNum)',beta);
     end;
   elseif egldtype==2, 
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
     % GLD+tanh used.
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      if (alfa4(signalNum)>borderBase+borderSlope*alfa3(signalNum)^2),
        FType(signalNum)=1; %GLD
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % GLD moment estimation used. See gld_momentfit.m for more details.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        lambda=gld_momentfit(alfa3(signalNum),alfa4(signalNum));
        [score1(signalNum,:) score2(signalNum,:)]=...
	    gld_score(U(:,signalNum)',lambda);  
      else 
        FType(signalNum)=3; %tanh       
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % tanh used. See FastICA (http://www.cis.hut.fi/projects/ica/fastica/)
	% and related papers for more details.
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        score1(signalNum,:)=tanh(Tanhconst * U(:,signalNum)');
        score2(signalNum,:)=1-score1(signalNum,:).^2;
      end;
   elseif egldtype==3,
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
     % GLD only used.
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      FType(signalNum)=1; %GLD
      lambda=gld_momentfit(alfa3(signalNum),alfa4(signalNum));
      [score1(signalNum,:) score2(signalNum,:)]=...
	  gld_score(U(:,signalNum)',lambda); 
    else
     error('EGLD contrast type unknown!');
    end 
    fprintf('[%d]', signalNum); 
  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(alfa3'))),...
                 strjust(strvcat('4th moment: ', num2str(alfa4'))),...
		 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%