%***********************************************************************%
% Nom du fichier : tfuni1.m															%
%																								%
% Ce fichier est tetan_nv mais il a été renommé									%
%***********************************************************************%

function [u,theta]=tfuni1(z,w);
% Search for a unitary direct transform, u,
% for separating 2 sources in presence of unknown noise.
% Contrast used=sum of marginal cumulants (without squares);
% Needs sources to have cumulants of the same sign, w.
% u: plane complex roatation separating the sources
% theta: complex angle of u
% Observation z and matrix u can be either real or complex,
% but if everything is real, use instead tfuni1R or tfuni4.
% Version P.Comon, July 23, 1996, 
% revised V.Capedevielle March 1997, 
% revised V.Capedevielle Sept 19, 1997
% revised P.Comon Aug 18, 1998
% Related bibliography: 
% Comon-Chevalier-Capdevielle, SPAWC'97.
% Comon-Moreau, ICASSP'97.
% Comon-Grellier, ICA'99, Aussois, France, Jan 11-15, 1999.

ii=sqrt(-1);
N=length(z);

tol=eps*10^4/N;

z1=z(1,:);
z1c=conj(z1);
z2=z(2,:);
z2c=conj(z2);

%Calcul des moments d'ordre 2
%-----------------------------------------
m11=mean(z1.*z1c);
m11n=mean(z1.*z1);
m22=mean(z2.*z2c);
m22n=mean(z2.*z2);
m12=mean(z1c.*z2);
m12n=mean(z1.*z2);

%Calcul des moments d'ordre 4 (le contraste ne fait intervenir que les cumulants
%ayant autant de puissances complexes que de puissances complexes conjuguees)
%--------------------------------------------------------------------------- ---------------------------------------

m1111=mean(z1.*z1c.*z1c.*z1);
m1112=mean(z1.*z1c.*z1c.*z2);
m2112=mean(z2.*z1c.*z1c.*z2);
m1122=mean(z1.*z1c.*z2c.*z2);
m2122=mean(z2.*z1c.*z2c.*z2);
m2222=mean(z2.*z2c.*z2c.*z2);

%Calcul des cumulants a l'ordre 4
%--------------------------------------------
c1111=m1111-2*m11^2-abs(m11n)^2;

c1112=m1112-2*m11*m12-m12n*conj(m11n);

c2112=m2112-2*m12^2-m22n*conj(m11n);

c1122=m1122-m11*m22-m12*conj(m12)-m12n*conj(m12n);

c2122=m2122-2*m12*m22-conj(m12n)*m22n;

c2222=m2222-2*m22^2-abs(m22n)^2;

c1122=real(c1122);c1111=real(c1111);c2222=real(c2222);


%disp('les valeurs des cumulants');
%disp(c1111);disp(c1112);disp(c2112);disp(c1122);disp(c2122);disp(c2222);


%variables auxiliaires
%---------------------
a=c1111+c2222;
b=c2122-c1112;
beta=angle(b);b=abs(b);b2=b*b;
delta=angle(c2112);d=abs(c2112);d2=d*d;
f=4*c1122+c1111+c2222;
mu=delta-2*beta;cm2=cos(mu)^2;
tempo=d*d*sin(mu)*cos(mu);


if b<tol,
r=exp(-ii*delta/2);r=[r ii*r 0];
Ups=[f/2+d f/2-d a];
[Umax,iopt]=max(w*Ups);theta=r(iopt);
u=(1/sqrt(2))*[1 theta; -theta' 1];
break;
end

%recherche de la phase phi du sinus de la rotation
%-------------------------------------------------------------------
%%recherche de la tangente de (phi+beta)/2
if abs(cm2-1)<tol,
 Phi=[];R=[];Ups=[];Ups0=[];Theta=[];
 Phi=[-delta/2;-delta/2+pi/2];
 Theta=exp(ii*Phi);
 Ups=[f/2+d;f/2-d];
 cpsi=[-1;1];
 sigma=(f-2*a+2*d*cos(mu))./(b*cpsi);
 R=(-sigma+sqrt(sigma.*sigma+16))/2;
 Phi=[pi-beta;-beta];
 r=(-R+sqrt(R.*R+4))/2;
 phi=Phi+pi;Theta0=r.*exp(ii*phi);
 Ups0=(a*R.*R+4*b*R.*cpsi+4*d*cos(mu)+2*f)./(R.*R+4);
 Theta=[Theta;Theta0;0];
 Ups=[Ups;Ups0;a];
 [Umax,iopt]=max(w*Ups);theta=Theta(iopt);
 u=[1 theta;-theta' 1]/sqrt(1+theta*theta');
 break
end

h(1)=(-1+cm2)*d2;
h(2)=(-f+2*a)*sin(mu)*d+6*tempo;
h(3)=(6-14*cm2)*d2+(-8*a+4*f)*cos(mu)*d+4*b2;
h(4)=(-10*a+5*f)*sin(mu)*d-14*tempo;
h(5)=0;
h(6)=h(4);
h(7)=-h(3);
h(8)=h(2);
h(9)=-h(1);

T=roots(h);T(find((imag(T))>0))=[];T=real(T);
%phase du sinus des rotations admissibles
Phi=2*atan(T)-beta;

%Recherche du module R de t-1/conj(t)
%cas ou phi+beta=0 est solution
Ups0=[];J=find(abs(sin(Phi+beta))<0.01);
if length(J)>0,R0=[];Phi0=[];
 for q=1:length(J),
  j=J(q);
  phij=Phi(j);
  cpb=cos(phij+beta);%peut valoir 1 ou -1
  Rj=roots([b*cpb,2*d*cos(2*phij+delta)+f-2*a,-4*b*cpb]);
  R0=[R0;Rj];Phi0=[Phi0;phij;phij];
 end
 JJ=find((imag(R0))>0);Phi0(JJ)=[];R0(JJ)=[];R0=real(R0);
 Phi(J)=[];%Il ne reste plus dans phi que les autres cas
end

if length(R0)>0,
 R02=R0.*R0;
 Ups0=a*R02+4*b*R0+(4*d*cos(mu)+2*f);Ups0=Ups0./(R02+4);
end

%Cas ou phi+beta est non nul
R=-2*d/b*sin(2*Phi+delta)./sin(Phi+beta);R2=R.*R;
Ups1=a*R2+4*b*cos(Phi+beta).*R+4*d*cos(2*Phi+delta)+2*f;
Ups1=Ups1./(R2+4);R=R(:);


%Selection du meilleur angle (donnant maxi contraste)

Ups=[Ups0;Ups1];
Phi1=Phi;Phi=[Phi0;Phi1(:)];
ro=[R0;R];
[Umax,iopt]=max(w*Ups);

%reconstitution de la rotation optimale
ro=ro(iopt);r=(ro+sqrt(ro*ro+4))/2;
theta=r*exp(ii*Phi(iopt));
if Umax<w*a,theta=0;r=0;Umax=w*a;end;


%%%%%%%%%%%%%%%%%%% LIGNE RAJOUTEE %%%%%%%%%%%%%%
%-CHOIX DE LA SOLUTION THETA DANS LE DISQUE UNITE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if r>1,theta=-1/conj(theta);r=abs(theta);end
cc=1/sqrt(1+r*r);
u(1,1)=cc;
u(2,2)=u(1,1);
u(1,2)=theta*cc;
u(2,1)=-conj(u(1,2));