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)