in3_out2.m

some thing more for MIMO

      % 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