% tracS3C2.m, P.Comon, July 2, 1998
% Idem que benchS3C2.m mais avec boucle sur plusieurs valeurs de T, et trace

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
TT=[200 300 500 1000 2000 3000 5000];
NBtirages=15;
for jT=1:length(TT),T=TT(jT);   %%%%%%%% BOUCLE SUR T
for tirage = 1:NBtirages,
fprintf('T = %g \n',T);
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;
X=fmins('objectif',[theta0,phi0]);

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;

Distance(jT,tirage)=gap(A,L);
end;  %% fin tirages
end;    %%%% FIN BOUCLE SUR T



MND=mean(Distance.');
figure(1);plot(TT,MND,'r',TT,MND,'or');
ylabel('Mean Gap');xlabel('Sample length T');
AX=axis;AX(3)=0.2;axis(AX);
print MeanGap.eps -deps

ecarts=Distance-MND(:)*ones(1,NBtirages);
MSTD=sqrt(mean((ecarts.^2).'));
figure(2);plot(TT,MSTD,'r',TT,MSTD,'or');
ylabel('Std Dev of Mean Gap');xlabel('Sample length T');
AX=axis;AX(3)=0.05;AX(4)=0.35;axis(AX);
print StdMeanGap.eps -deps