jade盲源分离算法

function [AS]=jade（Xm）
% Source separation of complex signals with JADE.
% Jade performs ‘Source Separation‘ in the following sense:
%   X is an n x T data matrix assumed modelled as X = A S + N where
% 
% o A is an unknown n x m matrix with full rank.
% o S is a m x T data matrix （source signals） with the properties
%    	a） for each t the components of S（:t） are statistically
%    	   independent
% 	b） for each p the S（p:） is the realization of a zero-mean
% 	   ‘source signal‘.
% 	c） At most one of these processes has a vanishing 4th-order
% 	   cumulant.
% o  N is a n x T matrix. It is a realization of a spatially white
%    Gaussian noise i.e. Cov（X） = sigma*eye（n） with unknown variance
%    sigma.  This is probably better than no modeling at all...
%
% Jade performs source separation via a 
% Joint Approximate Diagonalization of Eigen-matrices.  
%
% THIS VERSION ASSUMES ZERO-MEAN SIGNALS
%
% Input :
%   * X: Each column of X is a sample from the n sensors
%   * m: m is an optional argument for the number of sources.
%     If ommited JADE assumes as many sources as sensors.
%
% Output :
%    * A is an n x m estimate of the mixing matrix
%    * S is an m x T naive （ie pinv（A）*X）  estimate of the source signals
%
%
% Version 1.5.  Copyright: JF Cardoso.  
%
% See notes references and revision history at the bottom of this file



[nT]	= size（X）;

%%  source detection not implemented yet !
if nargin==1 m=n ; end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A few parameters that could be adjusted
nem	= m;		% number of eigen-matrices to be diagonalized
seuil	= 1/sqrt（T）/100;% a statistical threshold for stopping joint diag


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% whitening
%
if m 	[UD] 	= eig（（X*X‘）/T）; 
	[puissk]=sort（diag（D））;
 	ibl 	= sqrt（puiss（n-m+1:n）-mean（puiss（1:n-m）））;
 	bl 	= ones（m1） ./ ibl ;
 	W	= diag（bl）*U（1:nk（n-m+1:n））‘;
 	IW 	= U（1:nk（n-m+1:n））*diag（ibl）;
else    %assumes no noise
 	IW 	= sqrtm（（X*X‘）/T）;
 	W	= inv（IW）;
end;
Y	= W*X;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Cumulant estimation


R	= （Y*Y‘ ）/T ;
C	= （Y*Y.‘）/T ;

Yl	= zeros（1T）;
Ykl	= zeros（1T）;
Yjkl	= zeros（1T）;

Q	= zeros（m*m*m*m1） ;
index	= 1;

for lx = 1:m ; Yl 	= Y（lx:）;
for kx = 1:m ; Ykl 	= Yl.*conj（Y（kx:））;
for jx = 1:m ; Yjkl	= Ykl.*conj（Y（jx:））;
for ix = 1:m ; 
	Q（index） = ...
	（Yjkl * Y（ix:）.‘）/T -  R（ixjx）*R（lxkx） -  R（ixkx）*R（lxjx） -  C（ixlx）*conj（C（jxkx））  ;
	index	= index + 1 ;
end ;
end ;
end ;
end

%% If you prefer to use more memory and less CPU you may prefer this
%% code （due to J. Galy of ENSICA） for the estimation the cumulants
%ones_m = ones（m1） ; 
%T1 	= kron（ones_mY）; 
%T2 	= kron（Yones_m）;  
%TT 	= （T1.* conj（T2）） ;
%TS 	= （T1 * T2.‘）/T ;
%R 	= （Y*Y‘）/T  ;
%Q	= （TT*TT‘）/T - kron（Rones（m））.*kron（ones