function?[AS]=jade（Xm）
[nT]=size（X）;
?nem=m;
seuil=1/sqrt（T）/100;
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:mYl=Y（lx:）;
????for?kx=1:mYkl=Yl.*conj（Y（kx:））;
????????for?jx=1:mYjkl=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
[UD]=eig（reshape（Qm*mm*m））;
[laK]=sort（abs（diag（D）））;
M=zeros（mnem*m）;
Z=zeros（m）;
h=m*m;
for?u=1:m:nem*m
????Z（:）=U（:K（h））;
????M（:u:u+m-1）=la（h）*Z;
????h=h-1;
end;
%%%jiont?approximate?diagonalization?of?the?eigen-matrices
%%%Better