% 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

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];

% cum sources, supposes inconnus dans l'identification
Kurt=[1 1 1];   % cum symetriques 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.';
end;end;end;end;
T1111=T(1,1,1,1);T1112=T(1,1,1,2);T1122=T(1,1,2,2);
T1222=T(1,2,2,2);T2222=T(2,2,2,2);

% polynome normalise associe a T
P=[T1111 T1112 T1122 T1222 T2222];
d=4;fd=facto(d);c=ones(1,d+1);
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));

% 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;

%Theta=[0:0.1:1]*pi;
%Phi=[0:0.1:2]*pi;
%%Theta=[1.53:0.0005:1.535];
%%Phi=[1.7335:0.0005:1.7345];
%fprintf('\n it = ');
%for it=1:length(Theta),theta=Theta(it);fprintf('%g  ',it)
%  for ip=1:length(Phi);phi=Phi(ip);
%  smin(it,ip)=objectif([theta,phi]);
%end;end;fprintf('\n ');
%contour(Phi,Theta,smin);xlabel('Phi');ylabel('Theta');
%[sminI,I]=min(smin);[smin2,J]=min(sminI);
%theta0=Theta(I(J));phi0=Phi(J);
%fprintf('Smini = %g, pour theta = %g et phi = %g\n',smin2,theta0,phi0)

X=fmins('objectif',[theta0,phi0]);


%%%% 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;

OBTENU=normal(L)
VRAI=normal(A)
DISTANCE=gap(A,L)
fprintf('\n V.S.mini = %g, pour theta = %g et phi = %g\n',objectif(X),X(1),X(2))