www.pudn.com > zxb.rar > 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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%