自己用Matlab编写的ICA工具箱

function [Wefica, ISRef, Hef, Wsymm, ISRsymm, status]= efica(X, ini, FASTICApars, Saddlepars,EFICApars) % function [Wefica, ISRef, Hef, Wsymm, ISRsymm, status]=efica(X, ini, g, SaddleTest) %EFICA: [Wefica, ISRef, Hef, Wsymm, SIRsymm, status]=efica(X, ini, g, SaddleTest) % % version: 1.9 release: 24.11.2006 % % %improved speed % %Input data: % X ... mixed data dim x N, where dim is number of signals and N is number of % samples % ini ... starting point of the iterations for W matrix % g ... nonlinearity used in the symmetric algorithm ('tanh'-'rati'-'gaus'-'pow3') % SaddleTest ... if true (default) the test of saddle points on the demixed signals is done % % %Output data: % W_efica - demixing matrix produced by EFICA % ISRef - ISR matrix estimator of EFICA components % Hef - matrix for EFICA bias computation (in the noisy environment) % Wsymm - demixing matrix produced by symmetric Fast-ICA % ISRsymm - ISR matrix estimator of Fast-ICA components % status - one if test of saddle points was positive, zero otherwise % % Examples of usage: % % [Wefica, ISRef, Hef, Wsymm, SIRsymm]=efica(X); % % use default settings (random initialization, 'tanh' nonlinearity, % the test of saddle points will be done) % % % [Wefica, ISRef, Hef, Wsymm, SIRsymm]=efica(X, randn(dim), 'rati', true); % % random initialization, 'rati' nonlinearity used in symmetric FastICA, % the test of saddle points will be done % % % [Wefica, ISRef, Hef, Wsymm, SIRsymm]=efica(X, eye(dim), 'tanh', false); % % random initialization, 'tanh' nonlinearity used in symmetric FastICA, % without the test of saddle points % % % % % [dim N]=size(X); %Default values of parameters if nargin<5 %EFICApars % EFICApars = true; EFICAMaxIt=50; %%Maximum number of improving iterations eficaGaussianNL=3; %'ggda'; %Nonlinearity used for eficagaussian signals 'gaus'-'ggda'-'npnl' else EFICAMaxIt=EFICApars.EFICAMaxIt ; %%Maximum number of improving iterations eficaGaussianNL=EFICApars.eficaGaussianNL; %Nonlinearity used for eficagaussian signals 'gaus'-'ggda'-'npnl' end if nargin<4 %Saddlepars SaddleTest = 1; test_of_saddle_points_nonln = 4 ; %'rtsp'; MaxItAfterSaddleTest=30; %%Maximum number of iterations after a saddle point was indicated else SaddleTest = Saddlepars.SaddleTest ; test_of_saddle_points_nonln=Saddlepars.test_of_saddle_points_nonln; MaxItAfterSaddleTest=Saddlepars.MaxItAfterSaddleTest; %%Maximum number of iterations after a saddle point was indicated end if nargin<3 %FASTICApars g= 3; %'rati'; MaxIt=100; %% Maximum number of FastICA iterations else g=FASTICApars.g; MaxIt=FASTICApars.MaxIt; %% Maximum number of FastICA iterations end if nargin<2 while 1 ini=randn(dim); if cond(ini) < 10^3 break end end end epsilon=0.0001; %Stop criterion fineepsilon=1e-5; %Stop criterion for post-estimation repeat=1; rot2d=[1/sqrt(2) 1/sqrt(2);-1/sqrt(2) 1/sqrt(2)]; status=0; min_correlation=0.8; %additive noise...0.75, noise free... 0.95, turn off (unstable)...0 %removing mean X = X-mean(X,2)*ones(1,N); %fprintf('EFICA...\n'); %preprocessing C = cov(X'); CC = C^(-1/2); Z = CC*X; %FastICA W=ini; W_before_decor=W; W=W*real(inv(W'*W)^(1/2)); %Wold=zeros(dim,dim); crit=zeros(1,dim); NumIt=0; while repeat %%% Symmetric approach while (1-min(crit)>epsilon && NumIt<MaxIt)% && sum(double(changed>10))<2) Wold=W; switch g case 1 %'tanh' hypTan = tanh(Z'*W); W=Z*hypTan/N-ones(dim,1)*sum(1-hypTan.^2).*W/N; case 2 %'pow3' W=(Z*((Z'*W).^ 3))/N-3*W; case 3 %'rati' U=Z'*W; Usquared=U.^2; RR=4./(4+Usquared); Rati=U.*RR; Rati2=Rati.^2; dRati=RR-Rati2/2; nu=mean(dRati); % beta=mean(Rati2); hlp=Z*Rati/N; % mu=diag(W'*hlp)'; % J=ones(1,dim); tau=abs(nu-mu); gam=(beta-mu.^2); % ISR=(gam'*J+J'*gam+J'*tau.^2)./(tau'*J+J'*tau).^2; % ISR=ISR-diag(diag(ISR)); W=hlp-ones(dim,1)*nu.*W; case 4%'gaus' U=Z'*W; Usquared=U.^2; ex=exp(-Usquared/2); gauss=U.*ex; dGauss=(1-Usquared).*ex; W=Z*gauss/N-ones(dim,1)*sum(dGauss).*W/N; end % crit=abs(sum((W*diag(sqrt(1./sum(W.^2)))).*(W_before_decor*diag(sqrt(1./sum(W_before_decor.^2)))))); W_before_decor=W; W=W*real(inv(W'*W)^(1/2)); crit=abs(sum(W.*Wold)); NumIt=NumIt+1; end %while iteration fprintf('EFICA v1.9 NumIt: %d\n',NumIt); repeat=0; %%%The SaddleTest of the separated components if SaddleTest SaddleTest=false; %%The SaddleTest may be done only one times u=Z'*W; switch test_of_saddle_points_nonln case 1 %'tanh' table1=(mean(log(cosh(u)))-0.37456).^2; case 2 %'gaus' table1=(mean(u.*exp(-u.^2))-1/sqrt(2)).^2; case 3 %'rati' table1=(mean(2*log(4+u.^2))-3.16).^2; case 4 %'rtsp' table1=(mean(u.^2./(1+abs(u)))-0.4125).^2; case 5 %'pow3' table1=(mean((pwr(u,4)))-3).^2; end rotated=zeros(1,dim); checked=1:dim; for i=checked for j=checked(checked>i) if (~rotated(i) && ~rotated(j)) h=[u(:,i) u(:,j)]*rot2d; % rotate two columns i,j switch test_of_saddle_points_nonln case 1 %'tanh' ctrl=(mean(log(cosh(h)))-0.37456).^2; case 2 %'gaus' ctrl=(mean(exp(-h.^2/2)-1/sqrt(2))).^2; case 3 %'rati' ctrl=(mean(2*log(4+h.^2))-3.16).^2; case 4 %'rtsp' ctrl=(mean(h.^2./(1+abs(h)))-0.4125).^2; case 5 %'pow3' ctrl=(mean((pwr(h,4)))-3).^2; end if sum(ctrl)>table1(i)+table1(j) %bad extrem indicated rotated(i)=1;rotated(j)=1; %do not test the rotated signals anymore W(:,[i j])=W(:,[i j])*rot2d; repeat=1; %continue in iterating - the test of saddle points is positive NumIt=MaxIt-MaxItAfterSaddleTest; fprintf('EFICA: rotating components: %d and %d\n',i,j); status=1; end end end end end %if SaddleTest crit=zeros(1,dim); end %while repeat Wsymm=W'*CC; %estimated signals s=W'*Z; %estimate SIRs of the symmetric approach switch g case 1 %'tanh' mu=mean(s.*tanh(s),2); nu=mean(1./cosh(s).^2,2); beta=mean(tanh(s).^2,2); case 2 %'rati' ssquare=s.^2; mu=mean(ssquare./(1+ssquare/4),2); nu=mean((1-ssquare/4)./(ssquare/4+1).^2,2); beta=mean(ssquare./(1+ssquare/4).^2,2); case 3 %'gaus' aexp=exp(-s.^2/2); mu=mean(s.^2.*aexp,2); nu=mean((1-s.^2).*aexp,2); beta=mean((s.*aexp).^2,2); case 4 %'pow3' mu=mean(s.^4,2); nu=3*ones(dim,1); beta=mean(s.^6,2); end J=ones(1,dim); gam=(nu-mu)'; jm=(beta-mu.^2)'; Werr=(jm'*J+J'*jm+J'*gam.^2)./(abs(gam)'*J+J'*abs(gam)).^2; Werr=Werr-diag(diag(Werr)); ISRsymm=Werr/N; %SIRsymm=-10*l