function A=tfuni1R(e,w)
% A=tfuni1(e,w); 
% w: sign of the source kurtosis (-1 default)
% data e and transform A must be real
% e must be a (2xN) matrix
% S=A*e yields the estimated sources
% P.Comon, June 6 1997, revised Aug 19, 1998

if nargin<2,w=-1;end;
tol=eps*10^5; %mais on pourrait ajuster en fonction de N
N=length(e);
N2=N-1;N3=N2*(N-2)/N;
N4=N3*(N-3)/(N+1);N42=(N-2)*(N-3)/(N-1)/(N-1);;
m1=sum(e')'/N;
e=e-m1*ones(1,N);

%%%%% moments d'ordre 2
 m11=e(1,:)*e(1,:)'/N2;%cv vers 1 si standardise
 m22=e(2,:)*e(2,:)'/N2;%cv vers 1 si standardise
 m12=e(1,:)*e(2,:)'/N2;%cv vers 0 si standardise
 m21=m12;
%%%%% moments d'ordre 4
e2=e.*e;
 m1111=e2(1,:)*e2(1,:)'/N4;
 m1112=e2(1,:).*e(2,:)*e(1,:)'/N4;
 m1122=e2(1,:)*e2(2,:)'/N4;
 m1222=e2(2,:).*e(2,:)*e(1,:)'/N4;
 m2222=e2(2,:)*e2(2,:)'/N4;
%%%%% cumulants croises d'ordre 4
 g1111=m1111-N42*3*m11*m11;
 g1112=m1112-N42*3*m11*m12;
 g1122=m1122-N42*(2*m21*m12+m11*m22);
 g1222=m1222-N42*3*m12*m22;
 g2222=m2222-N42*3*m22*m22;

%%%%% variables auxiliaires
a=g1111+g2222;
b=g1222-g1112;b2=b*b;
d=g1122;d2=d*d;
f=4*g1122+g1111+g2222;

%% selection du meilleur R
R=roots([b,2*d+f-2*a,-4*b]);  % Cardan can be used to speed up rooting
Ups=a*R.*R+4*b*R+(4*d+2*f);Ups=Ups./(R.*R+4);
[Umax,iopt]=max(w*Ups);
ro=R(iopt);theta=(ro+sqrt(ro*ro+4))/2;

%%%% reconstitution de la rotation
if abs(theta)<1/N,
  A=eye(2); %fprintf('pas de rotation plane pour cette paire\n');
else,
  c=1/sqrt(1+theta*theta);
  A(1,1)=c;A(2,2)=A(1,1);A(1,2)=theta*c;A(2,1)=-conj(A(1,2));
end;

%%%%% filtrage de la sortie
% se ferait par: S=A*e;