cCFastICA.zip_cCfastica_fastica_fastica 复数_复数fastICA_独立成份分析

function [W,shat] = cCFastICA(x,A,s) % 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 distribution % % % % % % % % m = 50000; %number of observations % % % % j = 0; % % % % % % % % % some parameters for the available distributions % % % % bino1 = 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 one W = zeros(n,n); maxcounter = 40; for k = 1:n w = 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; end QAHw(:,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 break end; % isneg should always be positive (nonnegative) if Theorem 1 is fulfilled isneg(k) = EGmuutos(k)*Eg(index); W(:,k) = w; counters(k) = counter; end; %for k shat = W'*x; return;