% benchS3C2, P.Comon, July 2, 1998
% Idem que mainS3C2.m mais avec des cumulants estimes sur T echantillons

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


%%% GENERATION DES 3 SOURCES ET DES 2 OBSERVATIONS
T=500;
S=[vcircN(T,4);vcircN(T,4);vcircN(T,4)];
Y=A*S;

%%% CUMULANTS SYMETRIQUES COMPLEXES NON CIRCULAIRES
T11 = Y(1,:)*Y(1,:).'/T;
T12 = Y(1,:)*Y(2,:).'/T;T21=T12;
T22 = Y(2,:)*Y(2,:).'/T;
% Cumulant des observations (Symetrique Complexe)
T1111 = (Y(1,:).*Y(1,:))*(Y(1,:).*Y(1,:)).'/T;
T1112 = (Y(1,:).*Y(1,:))*(Y(1,:).*Y(2,:)).'/T;
T1122 = (Y(1,:).*Y(1,:))*(Y(2,:).*Y(2,:)).'/T;
T1222 = (Y(1,:).*Y(2,:))*(Y(2,:).*Y(2,:)).'/T;
T2222 = (Y(2,:).*Y(2,:))*(Y(2,:).*Y(2,:)).'/T;
T1111=T1111-3*T11*T11;
T1112=T1112-3*T11*T12;
T1122=T1112-2*T12*T12-T11*T22;
T1222=T1222-3*T12*T22;
T2222=T2222-3*T22*T22;

% 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 COMPLEXES CIRCULAIRES
TC11 = Y(1,:)*Y(1,:)'/T;
TC12 = Y(1,:)*Y(2,:)'/T;TC21=conj(TC12);
TC22 = Y(2,:)*Y(2,:)'/T;
% Cumulant des observations (Complexe Circulaire)
TC1111=(Y(1,:).*Y(1,:))*(Y(1,:).*Y(1,:))'/T;
TC1112=(Y(1,:).*Y(1,:))*(Y(1,:).*Y(2,:))'/T;
TC1122=(Y(1,:).*Y(1,:))*(Y(2,:).*Y(2,:))'/T;
TC1221=(Y(1,:).*Y(2,:))*(Y(2,:).*Y(1,:))'/T;
TC1222=(Y(1,:).*Y(2,:))*(Y(2,:).*Y(2,:))'/T;
TC2222=(Y(2,:).*Y(2,:))*(Y(2,:).*Y(2,:))'/T;
TC1111=TC1111-2*TC11*TC11-T11*conj(T11);
TC1112=TC1112-2*TC11*TC12-T11*conj(T12);
TC1122=TC1122-2*TC12*TC12-T11*conj(T22);
TC1221=TC1221-TC12*TC21-TC11*TC22-T12*conj(T21);
TC1222=TC1222-2*TC12*TC22-T12*conj(T22);
TC2222=TC2222-2*TC22*TC22-T22*conj(T22);

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