% Cumulant des observations (Complexe Circulaire)
      
      TC = T_22_true;
      TC1111=TC(1,1,1,1);TC2111=TC(2,1,1,1);TC2112=TC(2,1,1,2);
      TC2211=TC(2,2,1,1);TC2221=TC(2,2,2,1);TC2222=TC(2,2,2,2);
      
      % polynome reel associe a TC (en dim double: 4 variables)
      %PC=pzz2pX(TC1111,TC1112,TC1122,TC1221,TC1222,TC2222); %% Here we change 1 to 2, and vice versa. 
      PC=pzz2pX(TC2222,TC2221,TC2211,TC2112,TC2111,TC1111);  %% That is also change 2 to 1.
      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]);
      Optim_OPTIONS=optimset('MaxFunEvals', 2000);
      X=fminunc('objectif',[theta0,phi0],Optim_OPTIONS);
      
      
      %%%% PRESENTATION DU RESULTAT
      u=v1*cos(X(1))+v2*sin(X(1))*exp(i*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;
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%% Comon's code END here 2
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      
      
      
      
      %%% The system trasnfer function H(w) is
      %%% /                                  \
      %%% | L_roots(1) L_roots(2) L_roots(3) |
      %%% |     1         1           1      |
      %%% \                                  /
      
      %%% Note: The following Hroots and Lroots are actually the minus roots.
      [H2roots H2roots_index(w,:)]=sort(H(2,:)./H(1,:));
      [L2roots H2roots_index(w,:)]=sort(L(1,:));
      %H2roots,L2roots;
      poly_H2roots=poly(H2roots);
      poly_L2roots=poly(L2roots);
      poly_L2roots=poly_L2roots(4:-1:1)/poly_L2roots(4);
      H2_roots(w,:)=H2roots;
      L2_roots(w,:)=L2roots;
      

      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%% De Lathauwer's codes BEGIN here
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      Hr_true=H_roots(w,:)
      Lr_t=L_roots(w,:)
      [qest1roots, qest2roots, qest3roots, qest4roots] = lieven4(Tw_31(:,:,:,:,w), TC);
      
      
      Q1_roots(w,:)=qest1roots.';
      Q2_roots(w,1:length(qest2roots))=qest2roots.';
      Q3_roots(w,:)=qest3roots.';
      Q4_roots(w,:)=qest4roots.';

      
      %disp('Press Any Key when ready ...');
      %pause;
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%% De Lathauwer's codes END here
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      
            
   end %%% End of loop w=1:NF
   
   for ii=1:n
      L_roots_MSE(:,ii)=20*log10(abs(L_roots(:,ii)-H_roots(:,ii))./abs(H_roots(:,ii)));
      Q1_roots_MSE(:,ii)=20*log10(abs(Q1_roots(:,ii)-H_roots(:,ii))./abs(H_roots(:,ii)));
      Q2_roots_MSE(:,ii)=20*log10(abs(Q2_roots(:,ii)-H_roots(:,ii))./abs(H_roots(:,ii)));
      Q3_roots_MSE(:,ii)=20*log10(abs(Q3_roots(:,ii)-H_roots(:,ii))./abs(H_roots(:,ii)));
      Q4_roots_MSE(:,ii)=20*log10(abs(Q4_roots(:,ii)-H_roots(:,ii))./abs(H_roots(:,ii)));
   end
   
   for ii=1:n
      L2_roots_MSE(:,ii)=20*log10(abs(L2_roots(:,ii)-H2_roots(:,ii))./abs(H2_roots(:,ii)));
   end
   
   figure(1);clf;
   for ii=1:n
      subplot(n,2,(ii-1)*2+1);hold on;grid;
      plot(0:NF-1,real(H_roots(:,ii)),'r-');
      plot(0:NF-1,real(L_roots(:,ii)),'b-');
      if ii==1
         title('Real H_{11}(\omega)/H_{12}(\omega)');
      else 
         if ii==2
            title('Real H_{21}(\omega)/H_{22}(\omega)');
         else
            title('Real H_{31}(\omega)/H_{32}(\omega)');
         end
      end
      axis([1 NF -4 3]);
      
      subplot(n,2,(ii-1)*2+2);hold on;grid;
      plot(0:NF-1,imag(H_roots(:,ii)),'r-');
      plot(0:NF-1,imag(L_roots(:,ii)),'b-');
      if ii==1
         title('Imag H_{11}(\omega)/H_{12}(\omega)');
      else 
         if ii==2
            title('Imag H_{21}(\omega)/H_{22}(\omega)');
         else
            title('Imag H_{31}(\omega)/H_{32}(\omega)');
         end
      end
      axis([1 NF -4 3]);
   end
   
   figure(2);clf;
   for ii=1:n
      subplot(n,1,ii);hold on;%grid;
      plot((0:NF-1)*2*pi/NF,abs(H_roots(:,ii)),'r:');
      plot((0:NF-1)*2*pi/NF,abs(L_roots(:,ii)),'b-');
      if ii==1
         title('|H_{11}(\omega)/H_{12}(\omega)|');
      else 
         if ii==2
            title('|H_{21}(\omega)/H_{22}(\omega)|');
         else
            title('|H_{31}(\omega)/H_{32}(\omega)|');
         end
      end
      axis([0 2*pi 0 5]);
   end
      
   figure(3);clf;
   for ii=1:n
      subplot(n,2,(ii-1)*2+1);hold on;grid;
      plot(0:NF-1,real(1./H_roots(:,n+1-ii)),'r-');
      plot(0:NF-1,real(L2_roots(:,ii)),'b-');
      if ii==1
         title('Real H2_{11}(\omega)/H2_{12}(\omega)');
      else 
         if ii==2
            title('Real H2_{21}(\omega)/H2_{22}(\omega)');
         else
            title('Real H2_{31}(\omega)/H2_{32}(\omega)');
         end
      end
      axis([1 NF -4 3]);
      
      subplot(n,2,(ii-1)*2+2);hold on;grid;
      plot(0:NF-1,imag(1./H_roots(:,n+1-ii)),'r-');
      plot(0:NF-1,imag(L2_roots(:,ii)),'b-');
      if ii==1
         title('Imag H2_{11}(\omega)/H2_{12}(\omega)');
      else 
         if ii==2
            title('Imag H2_{21}(\omega)/H2_{22}(\omega)');
         else
            title('Imag H2_{31}(\omega)/H2_{32}(\omega)');
         end
      end
      axis([1 NF -4 3]);
   end
   
   figure(4);clf;
   for ii=1:n
      subplot(n,1,ii);hold on;grid;
      plot(0:NF-1,abs(1./H_roots(:,n+1-ii)),'r-');
      plot(0:NF-1,abs(L2_roots(:,ii)),'b-');
      if ii==1
         title('|H2_{11}(\omega)/H2_{12}(\omega)|');
      else 
         if ii==2
            title('|H2_{21}(\omega)/H2_{22}(\omega)|');
         else
            title('|H2_{31}(\omega)/H2_{32}(\omega)|');
         end
      end
      axis([1 NF 0 5]);
   end
      
   figure(5);clf;
   subplot(211);
   plot(L_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of L Roots in dB');
   subplot(212);
   plot(L2_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of L2 Roots in dB');
   
   figure(51);clf;
   subplot(321);
   plot(L_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of L Roots in dB');
   subplot(322);
   plot(L2_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of L2 Roots in dB');
   subplot(323);
   plot(Q1_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of Q1 Roots in dB');
   subplot(324);
   plot(Q2_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of Q2 Roots in dB');
   subplot(325);
   plot(Q3_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of Q3 Roots in dB');
   subplot(326);
   plot(Q4_roots_MSE);grid
   axis([1 NF -40 20]);
   title('MSE of Q4 Roots in dB');
   
   
   if USING_cum_true
      Error_index=find(abs(sum(H_roots-L_roots,2))>0.01).'
      Error_index_2=find(abs(sum(1./H_roots(:,3:-1:1)-L2_roots,2))>0.01).'
      Error_index_Q1=find(abs(sum(H_roots-Q1_roots,2))>0.01).'
      Error_index_Q2=find(abs(sum(H_roots-Q2_roots,2))>0.01).'
      Error_index_Q3=find(abs(sum(H_roots-Q3_roots,2))>0.01).'
      Error_index_Q4=find(abs(sum(H_roots-Q4_roots,2))>0.01).'
   else
      Error_index=find(abs(sum(H_roots-L_roots,2))>0.01).';
      Error_index_2=find(abs(sum(1./H_roots(:,3:-1:1)-L2_roots,2))>0.01).';
      Error_index_Q1=find(abs(sum(H_roots-Q1_roots,2))>0.01).';
      Error_index_Q2=find(abs(sum(H_roots-Q2_roots,2))>0.01).';
      Error_index_Q3=find(abs(sum(H_roots-Q3_roots,2))>0.01).';
      Error_index_Q4=find(abs(sum(H_roots-Q4_roots,2))>0.01).';
   end
   
   for w=Error_index
      w;
      H_roots(w,:);L_roots(w,:);
   end
   
   Error_number=length(Error_index)
   Error_number_2=length(Error_index_2)
   Error_number_Q1=length(Error_index_Q1)
   Error_number_Q2=length(Error_index_Q2)
   Error_number_Q3=length(Error_index_Q3)
   Error_number_Q4=length(Error_index_Q4)
   
   if DEBUG
      L_roots(Error_index,:)=H_roots(Error_index,:);
   end
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%% The roots of each column of H(w) has been got by now.
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%% The Following part is the frequency dependent scaling compensation
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%   Calculating the cross-power spectrum
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   rx = zeros(2*R_LENGTH+1,m,m);
   for ii=1:m
      for jj=1:m
         for seg=0:seg_num-1
            rx(:,ii,jj) = rx(:,ii,jj) + reshape(xcorr(...
               x(ii,seg_length*seg+1:seg_length*(seg+1)),...
               x(jj,seg_length*seg+1:seg_length*(seg+1)),R_LENGTH,'unbiased'),2*R_LENGTH+1,1);
         end
      end
   end
   rx=rx/seg_num;   
   
   if ADD_WINDOW
      for ii=1:m
         for jj=1:m
            rx(:,ii,jj)=rx(:,ii,jj).*r_window;
         end
      end
   end
   
   
   %% Estimating the cross Power spectrum
   Pw=zeros(NF,m,m);
   for ii=1:m
      for jj=1:m
         rr=zeros(NF,1);
         rr(NF/2-R_LENGTH+1:NF/2+R_LENGTH+1)=reshape(rx(:,ii,jj),2*R_LENGTH+1,1);
         Pw(:,ii,jj)=fft(fftshift(rr),NF);
      end
   end
   
   if Pw_symmetry
      for w=1:NF
         Pw_dummy=reshape((Pw(w,:,:)+Pw(w,:,:))/2,m,m);
         Pw(w,:,:)=(Pw_dummy+Pw_dummy')/2;
      end
   end
   
   Pw=shiftdim(Pw, 1);
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %%%   Calculating the cross-power spectrum ENDS here
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   if Using_DeLathauwer_ALS_method
      
      switch ALS_Method_index
      case {1}
         L_roots=Q1_roots;
         disp('ALS Method 1')
      case {2}
         L_roots=Q2_roots;
         disp('ALS Method 2')
      case {3}
         L_roots=Q3_roots;
         disp('ALS Method 3')
      case {4}
         L_roots=Q4_roots;
         disp('ALS Method 4')
      otherwise
         disp('Unknown ALS method.')
      end