function q=tfuni1V(A1,A2,Kum)
% q=tfuni1V(A1,A2,Kum)
% A1 and A2 are the rows of the mixing matrix for the pair considered.
% Seek for a unitary transform, q, for separating 2 sources
% in presence of unknown noise.
% Complex version of tfuniV.m, using a contrast without squares
% Source cumulants, Kum(k), should have the same sign.
% P.Comon, 25 Nov 1996
% Related bibliography: Comon-Moreau, ICASSP'97,
% Comon-Chevalier-Capdevielle SPAWC'97
% Comon Eurasip Signal Processing'94

ii=sqrt(-1);Kum=Kum(:);
w=sign(sum(Kum)); %% choix automatique: optionnel

% cumulants des observations
A1c=conj(A1);A2c=conj(A2);
A11=A1.*A1c;A12=A1.*A2c;A21=conj(A12);A22=A2.*A2c;
 g1111=(A11.*A11)*Kum; g1111=real(g1111);
 g1112=(A11.*A21)*Kum;
 g1122=(A11.*A22)*Kum; g1122=real(g1122);
 g1222=(A12.*A22)*Kum;
 g2222=(A22.*A22)*Kum; g2222=real(g2222);
 g2112=((A2.*A1c).^2)*Kum;
 g2122=conj(g1222);
g=[g1111 g1112 g1122 g1222 g2222 g2112];


%%%%% variables auxiliaires
a=g1111+g2222;
b=g2122-g1112;beta=angle(b);b=abs(b);b2=b*b;
delta=angle(g2112);d=abs(g2112);d2=d*d;
f=4*g1122+g1111+g2222;
mu=delta-2*beta;cm2=cos(mu)^2;
tempo=d*d*sin(mu)*cos(mu);
fprintf('Contraste en entree: %g\n',upsilon1(g,0)); %pour test

%%%%% CAS ou b=0 
if b<eps, fprintf('***cas b=0\n');
  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);
  q=[1 theta;-theta' 1]/sqrt(2);
  fprintf('Contraste en sortie: %g\n',Umax);%pour test
  break;
end;

%%%%% CAS sin(mu)=0
if abs(cm2-1)<eps, fprintf('***cas sin(mu)=0\n');
 Phi=[];R=[];Ups=[];Ups0=[];Theta=[];
 % alors T=[0;1;-1] ie mu=k*pi
 % si cos(psi)=0, ie T=+-1
 Phi=[-delta/2;-delta/2+pi/2];
 Theta=exp(ii*Phi);
 Ups=[f/2+d;f/2-d];
 % si sin(psi)=0, ie T=0
 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);
 % choix du meilleur theta
 Theta=[Theta;Theta0;0];
 Ups=[Ups;Ups0;a];
 [Umax,iopt]=max(w*Ups);theta=Theta(iopt);
 q=[1 theta;-theta' 1]/sqrt(1+theta*theta');
 fprintf('Contraste en sortie: %g\n',Umax);%pour test
 %fprintf('Contraste en sortie: %g\n',upsilon1(g,theta));%pour test
 break;
end;

%%%%% CAS generique
fprintf('***cas generique\n');
%%%%% RECHERCHE DE LA PHASE phi du sinus de la rotation
%% recherche de la tangente de (phi+beta)/2
 h(1)=(-1+cm2)*d2;                              %% t^8
 h(2)=(-f+2*a)*sin(mu)*d+6*tempo;		%% t^7
 h(3)=(6-14*cm2)*d2+(-8*a+4*f)*cos(mu)*d+4*b2;	%% t^6
 h(4)=(-10*a+5*f)*sin(mu)*d-14*tempo;           %% t^5
 h(5)=0;				        %% t^4
 h(6)=h(4);;		                        %% t^3
 h(7)=-h(3);                            	%% t^2
 h(8)=h(2);		                        %% t
 h(9)=-h(1);			         	%% 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 valloir 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;
  cc=1/sqrt(1+r*r);
  A(1,1)=cc;A(2,2)=A(1,1);A(1,2)=theta*cc;A(2,1)=-conj(A(1,2));
q=A;
%fprintf('Contraste en sortie: %g\n',upsilon1(g,theta));%pour test
fprintf('Contraste en sortie: %g\n',Umax);%pour test