user_alg6.m

The ICA/BSS algorithms are pure mathematical formulas, powerful, but rather mechanical procedures: T

function W=jade_fobi(X)
% JADE with random matrices 
%**************************************************
% Blind identification by joint diagonalization   
% of cumulant matrices generated by random matrices                       
%                       			  *
% REFERENCES:
% J.-P. Cardoso, "High Order Contrasts for independent Component Analysis", Neural Computation, 11, 157-192 
% (1999), 157-192.
% P. Georgiev, A, Cichocki, "Robust Blind Source Separation and Dispersing Algorithms", 
% ICASSP, Orlando, Floeida, 2002. 
% This code is addaptation of SOBI algorithm for cumulant matrices of fourth order
% by Pando Georgiev - January 8, 2003, modified by A. Cichocki March 5
% 2003.
% The code assume for simplicity that the number of sources is eaqual to the number of sensors.
p=0;
%X=x;
[m,N]=size(X);
L=20;
n=m;
%%% whitening
X=X-(mean(X')'*ones(1,N)); %Removing  mean value
 Rxx=(X*X')/N;
%  [Ux,Dx,Vx]=svd(Rxx);
%  Dx=diag(Dx);
% % n=11;
%  if n<m, % under assumption of additive white noise and
%         %when the number of sources are known or can a priori estimated
%   Dx=Dx-real((mean(Dx(n+1:m))));
%   Q= diag(real(sqrt(1./Dx(1:n))))*Ux(:,1:n)';
% else    % under assumption of no additive noise and when the 
%         % number of sources is unknown
%    n=max(find(Dx>1e-99)); %Detection the number of sources
%    Q= diag(real(sqrt(1./Dx(1:n))))*Ux(:,1:n)';
% end;
Q=inv(sqrtm(Rxx));
X=Q*X;
%Ag=Q*H;
[m,N]=size(X);
Xp=X(:,p+1:N);size(Xp);
X0=X(:,1:N-p);size(X0);
%---------------------------------
%Generating cumulant matrices
x4 = X0 * Xp' / (N-p);
x5 = Xp * X0' / (N-p);
x2 = X * X' / N;
%Joint diagonalization
pm=L*m;
k=1;
i=0;
for u=1:m:pm; i=i+1;
B=randn(m);
x1 = X0 * diag(sparse(sum(Xp .*(B*Xp), 1))) * X0' / (N-p);
x3	=   sum(sum((B*Xp) .*Xp, 1), 2) / (N-p);
C(:,:,i) = x1 - x2 * x3 -  x4 * (B+B') * x5;    
M(:,u:u+m-1)=(C(:,:,1+(u-1)/m));
end;
%%%joint diagonalization
seuil=1/sqrt(N)/100; encore=1; 
m;
U=eye(m);
yes=1;
while encore, encore=0;
 for p=1:m-1,
  for q=p+1:m,
   %%% Givens rotations
   g=[   M(p,p:m:pm)-M(q,q:m:pm)  ;
         M(p,q:m:pm)+M(q,p:m:pm)  ;
      i*(M(q,p:m:pm)-M(p,q:m:pm)) ];
	  [Ucp,D] = eig(real(g*g'));  [la,K]=sort(diag(D));
   angles=Ucp(:,K(3)); angles=sign(angles(1))*angles;
   c=sqrt(0.5+angles(1)/2);
   sr=0.5*(angles(2)-j*angles(3))/c;  sc=conj(sr);
   oui = abs(sr)>seuil ;
   encore=encore | oui ;
   if oui , %%%update of the A and V matrices 
       yes=yes+1;
    colp=M(:,p:m:pm); colq=M(:,q:m:pm);
    M(:,p:m:pm)=c*colp+sr*colq; M(:,q:m:pm)=c*colq-sc*colp;
    rowp=M(p,:); rowq=M(q,:);
    M(p,:)=c*rowp+sc*rowq; M(q,:)=c*rowq-sr*rowp;
    Up=U(:,p);
    U(:,p)=c*U(:,p)+sr*U(:,q);U(:,q)=c*U(:,q)-sc*Up;
   end%% if
  end%% q loop
 end%% p loop
end%% while
  W=U'*Q;