% mainS3C2, P.Comon, JUL 2, 1998
% Identification d'un melange de 3 sources recues sur 2 capteurs.
% On utilise d'abord le tenseur non ciculaire, pour trouver la
% decomposition, a une variete pres de dim 1. Puis on utilise
% le tenseur circulaire pour lever cette indetermination de dim 
% 1 (complexe) qui reste. cette indetermination est fixee a l'aide
% de pzz2pX.m carr2PR2.m et mod2PC2PR.m. La premiere decomposition
% utilise Quantics/binarydec3.m
%
%
% Modified by Binning Chen on March 11, 2001. To use C^{31} instead of C^{40}
%


clear
global VECTEURSNOYAU
global POLYNOMECIRCULAIRE

%path(path,'/home/comon/Matlab/Fonctions')
format compact;ii=sqrt(-1);jj=-1/2+ii*sqrt(3)/2;

phi11=pi/7;phi22=-pi/3;A11=exp(ii*phi11)*0.9;A22=exp(ii*phi22);
theta3=pi/3;phi13=pi/4;phi23=0;A13=cos(theta3)*exp(ii*phi13);
A23=sin(theta3)*exp(ii*phi23);
%A=[A11 0 A13;0 A22 A23];
randn(100);
A=[A11 0.8*randn A13;-0.97815*randn A22 A23*rand];

randn(100);

B=randn(2,3)+randn(2,3)*i

% cum sources, supposes inconnus dans l'identification
Kurt=[1 1 1];   % cum symetriques des sources
%KurtC=[-1 -1 -1];  % cum circulaires reels des sources
KurtC=[1 1 1];  % cum circulaires reels des sources

%%% CUMULANTS SYMETRIQUES NON CIRCULAIRES
% Cumulant des observations (Symetrique Complexe)
for i=1:2, 
   for j=1:2, 
      for k=1:2, 
         for l=1:2,
            %T(i,j,k,l) = (A(i,:).*A(j,:).*A(k,:).*A(l,:))*Kurt.';
            T(i,j,k,l) = (A(i,:).*A(j,:).*A(k,:).*conj(B(l,:)))*Kurt.';
         end;
      end;
   end;
end;
T1111=T(1,1,1,1);
T1112=T(1,1,1,2);
T2111=T(2,1,1,1);
T1122=T(1,1,2,2);
T2211=T(2,2,1,1);
T1222=T(1,2,2,2);
T2221=T(2,2,2,1);
T2222=T(2,2,2,2);

% polynome normalise associe a T
%P=[T1111 T1112 T1122 T1222 T2222];
P=[T1111 T2111 T2211 T2221 T2222];
d=4;% ADD by Binning Chen.% fd=facto(d);c=ones(1,d+1);
% ADD by Binning Chen.% for i=1:d-1,c(i+1)=fd/facto(i)/facto(d-i);end;
% le polynome associe a T est c.*P, soit P sous forme normalisee

% construction de la matrice de Sylvester
r=3;
%M=hankel(P(1:d-r+1),P(d-r+1:d+1));
M=[T1111 T2111 T2211 T2221; T1112 T1122 T1222 T2222];

% obtention des deux vecteurs du noyau
V=null(M);v1=V(:,1);v2=V(:,2);
VECTEURSNOYAU=[v1 v2];

%%% CUMULANTS CIRCULAIRES
% Cumulant des observations (Complexe Circulaire)
for i=1:2, 
   for j=1:2, 
      for k=1:2, 
         for l=1:2,
            TC(i,j,k,l) = (A(i,:).*A(j,:).*conj(A(k,:)).*conj(A(l,:)))*KurtC.';
         end;
      end;
   end;
end;
TC1111=TC(1,1,1,1);TC1112=TC(1,1,1,2);TC1122=TC(1,1,2,2);
TC1221=TC(1,2,2,1);TC1222=TC(1,2,2,2);TC2222=TC(2,2,2,2);

% polynome reel associe a TC (en dim double: 4 variables)
PC=pzz2pX(TC1111,TC1112,TC1122,TC1221,TC1222,TC2222);  
POLYNOMECIRCULAIRE=PC;
% PC est de taille 35=nombre de monomes de degre 4 en 4 variables


%%% RECHERCHE DE LA COMBINAISON OPTIMALE
smin=[];theta0=0;phi0=0;

TSTART_optimization = clock;

Optim_Options=optimset('MaxFunEvals', 3000);
%X=fmins('objectif',[theta0,phi0]);
X=fminunc('objectif',[theta0,phi0],Optim_Options);

TIME_optimization = etime(clock,TSTART_optimization)


%%%% PRESENTATION DU RESULTAT
u=v1*cos(X(1))+v2*sin(X(1))*exp(ii*X(2));
if abs(u(1))<2*eps, 
   q=roots(u(2:n));L=[1 q(:).';0 ones(1,2)];
else 
   q=roots(u); L=[q(:).';ones(1,3)]; 
end;

A_roots=sort(A(1,:)./A(2,:))
L_roots=sort(L(1,:))