in3_out2.m

some thing more for MIMO

   end 
    
    
   for w=1:NF 
       
      w; 
       
      H=HH(:,:,w); 
      Pw_true(:,:,w)=H*H'; 
       
       
      H_lambda(1,:)=abs(L_roots(w,:)).^2; 
      H_lambda(2,:)=L_roots(w,:); 
      H_lambda(3,:)=conj(L_roots(w,:)); 
      H_lambda(4,:)=[1 1 1]; 
       
      if Using_Pw_true 
         P_column=[Pw_true(1,:,w) Pw_true(2,:,w)].'; 
      else 
         P_column=[Pw(1,:,w) Pw(2,:,w)].'; 
      end 
       
      lambda=sqrt(inv(H_lambda'*H_lambda)*H_lambda'*P_column); 
       
      lambda=abs(lambda); 
       
      Hestp(:,:,w)=[L_roots(w,1) L_roots(w,2) L_roots(w,3); 1 1 1]*diag(lambda); 
       
      H_abs=abs(H); 
      Hestp_abs=abs(Hestp(:,:,w)); 
   end 
    
   if LIMIT_HESTP_MAGNITUDE 
      Hestp=min(5,abs(Hestp)).*exp(i*angle(Hestp)); 
   end 
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%% The frequency dependent scaling compensation ENDS here 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%% The Following part is the frequency dependent permutation 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
   for w=1:NF 
       
      w; 
       
      H=HH(:,:,w); 
       
      Cl_w=H*diag(abs(HH(l_index,:,Ref_Freq)).^2)*H'; 
       
      H_lambda(1,:)=abs(Hestp(1,:,w)).^2; 
      H_lambda(2,:)=Hestp(1,:,w).*conj(Hestp(2,:,w)); 
      H_lambda(3,:)=Hestp(2,:,w).*conj(Hestp(1,:,w)); 
      H_lambda(4,:)=abs(Hestp(2,:,w)).^2; 
       
      Cl_column=[Cl_w(1,:) Cl_w(2,:)].'; 
       
      lambda=inv(H_lambda'*H_lambda)*H_lambda'*Cl_column; 
       
      [lambda_sort lambda_index]=sort(lambda); 
       
      Hest(:,:,w)=Hestp(:,lambda_index,w); 
       
   end 
    
   if MODIFY_Hest_at_DC_Freq 
      Hest(:,:,1)=real(Hest(:,:,2)); 
      if mod(NF,2)==0 
         Hest(:,:,NF/2+1)=real(Hest(:,:,NF/2)); 
      end 
   end 
    
   figure(8);clf 
   for ii=1:m 
      for jj=1:n 
         subplot(m,n,(ii-1)*n+jj);hold on;%grid; 
         plot((0:NF-1)/NF*2*pi,reshape(abs(Hest(ii,H_order(jj),:)),NF,1),'b-'); 
         plot((0:NF-1)/NF*2*pi,reshape(abs(HH(ii,jj,:)),NF,1),'r:'); 
         if ii==1 & jj==1 
            title('|H_{11}(\omega)|'); 
            axis([0 2*pi 0 2.5]); 
         end 
         if ii==1 & jj==2 
            title('|H_{11}(\omega)|'); 
            axis([0 2*pi 0 5]); 
         end 
         if ii==1 & jj==3 
            title('|H_{11}(\omega)|'); 
            axis([0 2*pi 0 5]); 
         end 
         if ii==2 & jj==1 
            title('|H_{11}(\omega)|'); 
            axis([0 2*pi 0 5]); 
         end 
         if ii==2 & jj==2 
            title('|H_{11}(\omega)|'); 
            axis([0 2*pi 0 5]); 
         end 
         if ii==2 & jj==3 
            title('|H_{11}(\omega)|'); 
            axis([0 2*pi 0 3]); 
         end 
         xlabel('Frequency (\omega)'); 
      end 
   end 
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%% The frequency dependent permutation ENDS here 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%% The Following part phase ambiguity recovery.  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
   if 0 
      Hest_bak=Hest; 
      Hest(:,:,Error_index)=HH(:,H_order_i,Error_index); 
   end 
    
   Phase_hat=angle(Hest);       %%% the estimated phase with phase ambiguity 
    
   for w=1:NF 
       
      C_phase_w=conj(HH(:,:,w))*... 
         diag(HH(l_1, :, mod(1-k_arfa-w3, NF)+1).*conj(HH(l_2, :, mod(1-w3, NF)+1)))*... 
         HH(:,:,mod(k_arfa+w-1, NF)+1).'; 
       
       
      H_lambda(1,:)=conj(Hest(1,:,w)).*Hest(1,:,mod(w+k_arfa-1,NF)+1); 
      H_lambda(2,:)=conj(Hest(1,:,w)).*Hest(2,:,mod(w+k_arfa-1,NF)+1); 
      H_lambda(3,:)=conj(Hest(2,:,w)).*Hest(1,:,mod(w+k_arfa-1,NF)+1); 
      H_lambda(4,:)=conj(Hest(2,:,w)).*Hest(2,:,mod(w+k_arfa-1,NF)+1); 
       
      C_phase_column=[C_phase_w(1,:) C_phase_w(2,:)].'; 
       
      lambda_phase(:,w)=inv(H_lambda'*H_lambda)*H_lambda'*C_phase_column; 
       
       
   end 
    
   Psi=angle(lambda_phase); 
    
   Phase_1_sum=sum(Psi,2); 
   Phase_1_est=-Phase_1_sum/NF; 
    
   Phi=zeros(n,NF); 
    
   for ii=1:n 
      Phi(ii,2:NF) = (A1*reshape(Psi(ii,2:NF),NF-1,1)+Phase_1_est(ii)*sum(A1,2)).'; 
   end 
    
   for w=1:NF 
      for ii=1:n 
         Phase_est(:,ii,w)=Phase_hat(:,ii,w)+Phi(ii,w)*ones(m,1); 
      end 
   end 
    
   for ii=1:m 
      for jj=1:n 
         Phase_est_unwrap(ii,jj,:)=unwrap(Phase_est(ii,jj,:)); 
         Phase_true_unwrap(ii,jj,:)=unwrap(Phase_true(ii,jj,:)); 
      end 
   end 
    
   figure(10); clf;   %%% Plot the unwraped phase 
   for ii=1:m 
      for jj=1:n 
         subplot(m,n,(ii-1)*n+jj); 
         title(sprintf('Phase Est and True Unwrap',ii,jj));grid;hold on; 
         plot(reshape((Phase_est_unwrap(ii,H_order(jj),:)),NF,1)/pi,'b-'); 
         plot(reshape((Phase_true_unwrap(ii,jj,:)),NF,1)/pi,'r-'); 
         ylabel('Phase in PI'); 
         Current_axis=axis; 
         Current_axis(1)=1; 
         Current_axis(2)=NF; 
         axis(Current_axis); 
      end 
   end    
    
   figure(11); clf;   %%% Plot the unwraped phase difference (Linear) 
   for ii=1:m 
      for jj=1:n 
         subplot(m,n,(ii-1)*n+jj); 
         title(sprintf('Phase Est and True Unwrap difference',ii,jj));grid;hold on; 
         plot(reshape((Phase_est_unwrap(ii,H_order(jj),:)),NF,1)/pi-... 
              reshape((Phase_true_unwrap(ii,jj,:)),NF,1)/pi,'b-'); 
         ylabel('Phase in PI'); 
         Current_axis=axis; 
         Current_axis(1)=1; 
         Current_axis(2)=NF; 
         axis(Current_axis); 
      end 
   end    
    
   Delay_Times=[2 4 4] 
   Phase_shift=[1 1 0] 
    
   for ii=1:n 
      for w=1:NF 
         Phase_est_mod(:,H_order(ii),w)=mod(Phase_est(:,H_order(ii),w)-Delay_Times(ii)*pi*(w-1)/NF+... 
            Phase_shift(ii)*pi+pi, 2*pi)-pi; 
      end 
   end 
    
   figure(12); clf;   %%% Plot the phase 
   for ii=1:m 
      for jj=1:n 
         subplot(m,n,(ii-1)*n+jj); 
         title('Phase Estimation'); 
         hold on; 
         plot((0:NF-1)*2*pi/NF, reshape((Phase_est_mod(ii,H_order(jj),:)),NF,1),'b-'); 
         plot((0:NF-1)*2*pi/NF, reshape((Phase_true(ii,jj,:)),NF,1),'r:'); 
         xlabel('Frequency (\omega)'); 
         axis([0 2*pi -pi pi]); 
      end 
   end 
    
   Hest_all=abs(Hest).*exp(i*Phase_est_mod); 
    
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%% The Phase estimation ENDS here 
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%% The Following is time-momain impluse response retrieval using IFFT.  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
   if DEBUG 
      Hest_all=abs(HH(:,H_order_i,:)).*exp(i*Phase_est_mod); 
   end 
    
    
   for ii=1:m 
      for jj=1:n 
         hest_ifft(:,ii,jj)=real(fftshift(ifft(shiftdim(Hest_all(ii,jj,:),2)))); 
      end 
   end 
    
   hest_ifft_short=zeros(2*CL+1,m,n); 
   h_search_distance=max(2,CL-0); 
    
   for ii=1:m 
      for jj=1:n 
         [h_true_max h_true_max_index]=max(abs(h(:,ii,jj))); 
         h_long=zeros(2*CL+1,1); 
         h_long((CL+1:2*CL)+1-h_true_max_index)=h(:,ii,jj); 
         h_long_matrix=repmat(h_long,1,2*h_search_distance+1); 
          
         hest_ifft_dummy=hest_ifft(:,ii,H_order(jj)); 
         if Limit_Delay_time 
            Delay_limit=round(NF/5); 
            hest_ifft_dummy(1:Delay_limit)=0; 
            hest_ifft_dummy(NF-Delay_limit:NF)=0; 
         end 
         [h_max h_max_index]=max(abs(hest_ifft_dummy)); 
 
         h_candidate=zeros(2*CL+1,2*h_search_distance+1); 
         for h_index=h_max_index-h_search_distance:h_max_index+h_search_distance 
            h_candidate(:,h_index-h_max_index+h_search_distance+1) = hest_ifft_dummy(h_index-CL:h_index+CL); 
         end 
          
         [h_error_p h_p_index]=min(sum(abs(h_candidate - h_long_matrix).^2,1)); 
         [h_error_n h_n_index]=min(sum(abs(h_candidate + h_long_matrix).^2,1)); 
          
         if h_error_p <= h_error_n 
            h_ifft_short=hest_ifft_dummy(h_max_index+h_p_index-h_search_distance-1+(-CL:CL)); 
         else 
            h_ifft_short=-hest_ifft_dummy(h_max_index+h_n_index-h_search_distance-1+(-CL:CL)); 
         end 
          
         hest_ifft_short(:,ii,H_order(jj))=h_ifft_short; 
      end 
   end 
    
   figure(14);clf; 
   for ii=1:m 
      for jj=1:n 
         subplot(m,n,(ii-1)*n+jj);hold on;%grid; 
         [h_true_max h_true_max_index]=max(abs(h(:,ii,jj))); 
         h_long=zeros(2*CL+1,1); 
         h_long((CL+1:2*CL)+1-h_true_max_index)=h(:,ii,jj); 
         plot(0:2*CL, h_long,'r:o'); 
         plot(0:2*CL, hest_ifft_short(:,ii,H_order(jj)),'b-*'); 
         %legend('True', 'Estimation'); 
         title(sprintf('h_{%d%d}(n)',ii,jj)); 
         AX=axis; 
         AX(1)=0; 
         AX(2)=2*CL; 
         axis(AX); 
      end 
   end 
    
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
   %%% The time-momain impluse response retrieval using IFFT ENDS here.  
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
   L_roots_his(:,:,run)=L_roots; 
   L2_roots_his(:,:,run)=L2_roots; 
    
   Q1_roots_his(:,:,run)=Q1_roots; 
   Q2_roots_his(:,:,run)=Q2_roots; 
   Q3_roots_his(:,:,run)=Q3_roots; 
   Q4_roots_his(:,:,run)=Q4_roots; 
 
   hest_ifft_short_his(:,:,:,run)=hest_ifft_short; 
    
   Hest_all_his(:,:,:,run)=Hest_all; 
    
end %%% End of loop run 
 
hest_ifft_short_mean=mean(hest_ifft_short_his,4); 
hest_ifft_short_std=std(hest_ifft_short_his,0,4); 
 
L_roots_mean=mean(L_roots_his,3); 
L_roots_std=std(L_roots_his,0,3); 
 
L2_roots_mean=mean(L2_roots_his,3); 
L2_roots_std=std(L2_roots_his,0,3); 
 
Q1_roots_mean=mean(Q1_roots_his,3); 
Q1_roots_std=std(Q1_roots_his,0,3); 
Q2_roots_mean=mean(Q2_roots_his,3); 
Q2_roots_std=std(Q2_roots_his,0,3); 
Q3_roots_mean=mean(Q3_roots_his,3); 
Q3_roots_std=std(Q3_roots_his,0,3); 
Q4_roots_mean=mean(Q4_roots_his,3); 
Q4_roots_std=std(Q4_roots_his,0,3); 
 
SUM_ABS_cum_true_all 
SUM_ABS_cum_all_error 
 
Total_time = etime(clock, Program_Start)