% Complex FastICA% Ella Bingham 1999% Neural Networks Research Centre, Helsinki University of Technology% This is simple Matlab code for computing FastICA on complex valued signals.% The algorithm is reported in:% Ella Bingham and Aapo Hyv鋜inen, "A fast fixed-point algorithm for % independent component analysis of complex valued signals", International % Journal of Neural Systems, Vol. 10, No. 1 (February, 2000) 1-8.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%eps = 1;	%epsilon in G% Generate complex signals s = r .* (cos(f) + i*sin(f)) where f is uniformly% distributed and r is distributed according to a chosen distributionm = 50000; %number of observationsj = 0;% some parameters for the available distributionsbino1 = max(2, ceil(20*rand)); bino2 = rand;exp1 = ceil(10*rand);gam1 = ceil(10*rand); gam2 = gam1 + ceil(10*rand);f1 = ceil(10*rand); f2 = ceil(100*rand);poiss1 = ceil(10*rand);nbin1 = ceil(10*rand); nbin2 = rand;hyge1 = ceil(900*rand); hyge2 = ceil(20*rand); hyge3 = round(hyge1/max(2,ceil(5*rand)));chi1 = ceil(20*rand);beta1 = ceil(10*rand); beta2 = beta1 + ceil(10*rand);unif1 = ceil(2*rand); unif2 = unif1 + ceil(2*rand);gam3 = ceil(20*rand); gam4 = gam3 + ceil(20*rand);f3 = ceil(10*rand); f4 = ceil(50*rand);exp2 = ceil(20*rand);rayl1 = 10;unid1 = ceil(100*rand);norm1 = ceil(10*rand); norm2 = ceil(10*rand);logn1 = ceil(10*rand); logn2 = ceil(10*rand);geo1 = rand;weib1 = ceil(10*rand); weib2 = weib1 + ceil(10*rand);j = j + 1; r(j,:) = binornd(bino1,bino2,1,m); j = j + 1; r(j,:) = gamrnd(gam1,gam2,1,m); j = j + 1; r(j,:) = poissrnd(poiss1,1,m); j = j + 1; r(j,:) = hygernd(hyge1,hyge2,hyge3,1,m); j = j + 1; r(j,:) = betarnd(beta1,beta2,1,m); %j = j + 1; r(j,:) = exprnd(exp1, 1, m); %%j = j + 1; r(j,:) = unidrnd(unid1,1,m); %j = j + 1; r(j,:) = normrnd(norm1,norm2,1,m); %j = j + 1; r(j,:) = geornd(geo1,1,m); [n,m] = size(r);for j = 1:n	f(j,:) = unifrnd(-2*pi, 2*pi, 1, m);end;s = r .* (cos(f) + i*sin(f));[n,m] = size(s);% Whitening of s:s = inv(diag(std(s'))) * s;% Mixing using complex mixing matrix A: each coefficient a_{jk} is complex.A = rand(n,n) + i*rand(n,n);%A = orth(A);xold = A * s;% Whitening of x:[Ex, Dx] = eig(cov(xold'));Q = sqrt(inv(Dx)) * Ex';x = Q * xold;%x = x - mean(x,2)*ones(1,m);% Condition in Theorem 1 should be < 0 when maximising and > 0 when % minimising E{G(|w^Hx|^2)}. for j = 1:n;	g(j,:) = 1./(eps + abs(s(j,:)).^2) .* (1 - abs(s(j,:)).^2 .* ...		(1/2 * 1./(eps + abs(s(j,:)).^2) + 1));	Eg(j) = mean(g(j,:));end;% Fixed point algorithm: components one by oneW = zeros(n,n);maxcounter = 40;for k = 1:nw = rand(n,1) + i*rand(n,1);clear EG;counter = 0;wold = zeros(n,1);absAHw(:,k) = ones(n,1);while min(sum(abs(abs(wold) - abs(w))), maxcounter - counter) > 0.001;	wold = w;	g = 1./(eps + abs(w'*x).^2);	dg = -1/2 * 1./(eps + abs(w'*x).^2).^2;	w = mean(x .* (ones(n,1)*conj(w'*x)) .* (ones(n,1)*g), 2) - ...		mean(g + abs(w'*x).^2 .* dg) * w;	w = w / norm(w);	% Decorrelation:	w = w - W*W'*w;	w = w / norm(w);	counter = counter + 1;	G = log(eps + abs(w'*x).^2);	EG(counter) = mean(G,2);	if k < n		figure(3), subplot(floor(n/2),2,k), plot(EG)	end;endQAHw(:,k) = (Q*A)'*w;absQAHw = abs(QAHw);% which component was found  = index[maximum, index] = max(absQAHw(:,k));% Is EG increasing or decreasing; did we find a maximum or a minimum?EGmuutos(k) = EG(counter) - EG(counter-1);if EGmuutos(k) == NaN	breakend;% isneg should always be positive (nonnegative) if Theorem 1 is fulfilledisneg(k) = EGmuutos(k)*Eg(index);W(:,k) = w;counters(k) = counter;end; %for kshat = W'*x;