tito_main.m

用MATLAB实现MIMO系统盲辨识

   % Computing W = svd(Ry1*Ry1')
   %===================================
   W = zeros(n,n,NF/2+1);          %W(signalIndex, signalIndex, frequency)
   
   for (w=1:NF/2+1)
      [W(:,:,w),Sdummy,W2(:,:,w)] = svd(Ry1(:,:,w));
      %[W(:,:,w),Sdummy] = eig(Cy1(:,:,w));
   end
   
   
   %
   % Estimate H
   %================================
   Hest = zeros(n,n,NF/2+1);
   for (w=1:NF/2+1)
      M = inv(V(:,:,w)) * W(:,:,w);
      Hest(:,:,w) = M * inv(diag(diag(M)));
   end
   
   %
   % Test ratio H_{21}(w), H_{12}(w)
   %==================================
   trueS = zeros(1,NF/2+1);       %%% trueS and trueRatio are the two invariants defined in
   trueRatio = zeros(1,NF/2+1);   %%% in the paper referd.
   S = zeros(1,NF/2+1);           %%% S and R are the corresponding estimates.
   R = zeros(1,NF/2+1);
   
   for (w=1:NF/2+1)
      R(w) = Hest(2,1,w) / Hest(1,2,w);
      S(w) = Hest(2,1,w)+ 1/ Hest(1,2,w);
      
      trueRatio(w) = H(w,2,1) / H(w,1,2);
      trueS(w) = H(w,2,1)+ 1/ H(w,1,2);
   end
   S_his(iloop,:)=S;
   R_his(iloop,:)=R;
   
   R=[R fliplr(conj(R(2:NF/2)))];
   S=[S fliplr(conj(S(2:NF/2)))];
   
   %.......rec of channels for R and S......
   
   num=2*length(h1)-1;
   
   %%%%%%% The following is the code for the estimation of the E(w) defined in the paper.
   
   FF=abs(S).^2+abs(sqrt(S.^2-4*R)).^2;
   DD=2*real(conj(S).*sqrt(S.^2-4*R));
   
   [P_W_t ls_err1]=rec_fromphase_part(phase(R),NF,num,Rx_cond_Index);
   P_W_A=abs(fft(P_W_t,NF));
   
   HH1_abs=abs(sqrt(P_W_A./R));
   HH2_abs=abs(sqrt(P_W_A.*R));
   
   I_HH1_abs=HH1_abs;
   I_HH2_abs=HH2_abs;
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %% The following is a iteration using new HH2_abs %%
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   for H2_loop_index=1:Iter_HH2_abs+1
      
      [H2_max index_1]=max(HH2_abs);
      index_2=index_1+1;
      if index_1 == NF
         index_2=index_1-1;
      end
      ii=(FF(1)./HH2_abs(1)^2-FF(2)./HH2_abs(2)^2)/ (abs(DD(2))/ HH2_abs(2).^2 -abs(DD(1))/ HH2_abs(1).^2);
      ii=(FF(index_1)./HH2_abs(index_1)^2-FF(index_2)./HH2_abs(index_2)^2)/ (abs(DD(index_2))/ HH2_abs(index_2).^2 -abs(DD(index_1))/ HH2_abs(index_1).^2);
      
      ii=sign(ii);   % ii is the s in the paper
      
      INDEX=find(HH2_abs>max(HH2_abs)*10^(-6/20));
      c=FF./HH2_abs.^2 +ii*abs(DD)./HH2_abs.^2;
      cc=FF(INDEX)./HH2_abs(INDEX).^2 +ii*abs(DD(INDEX))./HH2_abs(INDEX).^2;
      c=sqrt(mean(cc)/4);
      
      EE=(4*(c*HH2_abs).^2 -abs(S).^2-abs(sqrt(S.^2-4*R)).^2)./(2*real(conj(S).*sqrt(S.^2-4*R)));
      
      EE=sign(real(EE));
      
      H2_r=.5*S +.5*EE.*sqrt(S.^2-4*R);
      H1_r=H2_r./R;
      
      
      if Minimum_Phase_Channel
         h1_r=rec_frommag_complex(H1_r,Rx_cond_Index,L);
         h2_r=rec_frommag_complex(H2_r,Rx_cond_Index,L);
      else
         h1_r=rec_from_partial(H1_r,Rx_cond_Index,L);
         h2_r=rec_from_partial(H2_r,Rx_cond_Index,L);
      end
      
      h1_r=real(h1_r(1:L))*sign(real(h1_r(1))); 
      h1_r=h1_r*norm(h1)/norm(h1_r);
      h2_r=real(h2_r(1:L))*sign(real(h2_r(1))); 
      h2_r=h2_r*norm(h2)/norm(h2_r);
      
      h1_p_his(iloop,:)=h1_r;
      h2_p_his(iloop,:)=h2_r;
      
      MSE_p_h1(iloop)=20*log10(norm(h1_r-h1)/norm(h1));
      MSE_p_h2(iloop)=20*log10(norm(h2_r-h2)/norm(h2));
      
      h1_r_i=h1_r;
      h2_r_i=h2_r;
      EE_i=EE;
      
      
      for iter=1:ITER_TIMES
         
         ntimes(1,iter)=0;
         
         h1_l_norm=norm(h1_r_i);
         h2_l_norm=norm(h2_r_i);
         
         h1_r_i=real(h1_r_i(1:L));
         h2_r_i=real(h2_r_i(1:L));
         
         h1_norm=norm(h1_r_i);
         h2_norm=norm(h2_r_i);
         h1_r_i=h1_r_i*h1_l_norm/h1_norm;
         h2_r_i=h2_r_i*h2_l_norm/h2_norm;
    
         H1_r_i=fft(h1_r_i,NF);
         H2_r_i=fft(h2_r_i,NF);
    
         H2_r_c=.5*S +.5*EE_i.*sqrt(S.^2-4*R);
         H1_r_c=1./(.5*S -.5*EE_i.*sqrt(S.^2-4*R));
         for w=1:NF
            ALT1=0.5*(S(w)+sqrt(S(w)^2-4*R(w)));
            ALT2=0.5*(S(w)-sqrt(S(w)^2-4*R(w)));
            if abs(abs(ALT1)-HH2_abs(w)) < abs(abs(ALT2)-HH2_abs(w))
               H2_r_i(w)=ALT1;
            else
               H2_r_i(w)=ALT2;
            end
         end
         
         H1_r_i=H2_r_i./R;
         
         
         if Minimum_Phase_Channel
            h1_r_i=rec_frommag_complex(H1_r_i,Rx_cond_Index,L);
            h2_r_i=rec_frommag_complex(H2_r_i,Rx_cond_Index,L);
         else
            h1_r_i=rec_from_partial(H1_r_i,Rx_cond_Index,L);
            h2_r_i=rec_from_partial(H2_r_i,Rx_cond_Index,L);
         end
         
         h1_c=real(h1_r_i(1:L))*sign(real(h1_r_i(1)));h1_c=h1_c*norm(h1)/norm(h1_c);
         h2_c=real(h2_r_i(1:L))*sign(real(h2_r_i(1)));h2_c=h2_c*norm(h2)/norm(h2_c);
         
         MSE_h1(iter,iloop)=20*log10(norm(h1_c-h1)/norm(h1));
         MSE_h2(iter,iloop)=20*log10(norm(h2_c-h2)/norm(h2));
         
      end    %  End of Iteration on EE
      
      h1_r_i=h1_r_i(1:L);
      h1_r_i=h1_r_i*norm(h1)/norm(h1_r_i);
      
      h2_r_i=h2_r_i(1:L);
      h2_r_i=h2_r_i*norm(h2)/norm(h2_r_i);
      
      HH1_abs=abs(fft(real(h1_r_i),128));
      HH2_abs=abs(fft(real(h2_r_i),128));
      
      
      figure(1);clf
      subplot(221);
      plot(1:L,h1,'r-',1:L,real(h1_r_i)*sign(real(h1_r_i(1))),'b-');
      grid;title('h1');legend('true h1','Estimation');
      subplot(222);plot(1:L,h2,'r-',1:L,real(h2_r_i)*sign(real(h2_r_i(1))),'b-');
      grid;title('h2');legend('true h1','Estimation');
      subplot(223);plot(1:128,abs(H1),'r-',1:128,abs(I_HH1_abs),'g-',1:128,abs(HH1_abs),'b-');
      grid;title('H1');legend('true |H1(w)|','Initial Estimation','Final Estimation');
      subplot(224);plot(1:128,abs(H2),'r-',1:128,abs(I_HH2_abs),'g-',1:128,abs(HH2_abs),'b-');
      grid;title('H2');legend('true |H2(w)|','Initial Estimation','Final Estimation');
      drawnow;
      
   end  % End of H2_loop_index
   
   h1_his(iloop,:)=real(h1_r_i)*sign(real(h1_r_i(1)));
   h2_his(iloop,:)=real(h2_r_i)*sign(real(h2_r_i(1)));
   
   if (iloop >3 )
      h1_p_mean=mean(h1_p_his,1);
      h1_p_std=std(h1_p_his,1);
      
      h2_p_mean=mean(h2_p_his,1);
      h2_p_std=std(h2_p_his,1);
      
      h1_mean=mean(h1_his,1);
      h1_std=std(h1_his,1);
      
      h2_mean=mean(h2_his,1);
      h2_std=std(h2_his,1);
      
      figure(2);clf;
      subplot(221);shadow_plot(h1, h1_p_mean,h1_p_std);grid;title('original h1');axis_1=axis;
      subplot(223);shadow_plot(h2, h2_p_mean,h2_p_std);grid;title('original h2');axis_2=axis;
      subplot(222);shadow_plot(h1, h1_mean,h1_std);grid;title('Iterated h1');axis(axis_1);
      subplot(224);shadow_plot(h2, h2_mean,h2_std);grid;title('Iterated h2');axis(axis_2);
   end
   
end %iloop

MSE_h1_mean = mean(MSE_h1')
MSE_h1_std  = std(MSE_h1')

MSE_h2_mean = mean(MSE_h2')
MSE_h2_std  = std(MSE_h2')

S_mean=mean(S_his,1);
S_std_real=std(real(S_his),0,1);
S_std_imag=std(imag(S_his),0,1);
R_mean=mean(R_his,1);
R_std_real=std(real(R_his),0,1);
R_std_imag=std(imag(R_his),0,1);

figure(3);clf;
subplot(221);shadow_plot(real(trueS), real(S_mean),S_std_real);grid;title('Real S');
subplot(222);shadow_plot(imag(trueS), imag(S_mean),S_std_imag);grid;title('Imag S');

subplot(223);shadow_plot(real(trueRatio), real(R_mean),R_std_real);grid;title('Real R');
subplot(224);shadow_plot(imag(trueRatio), imag(R_mean),R_std_imag);grid;title('Imag R');


S_MSE=20*log10(norm(S_mean-trueS)/norm(trueS))
R_MSE=20*log10(norm(R_mean-trueRatio)/norm(trueRatio))

Sh=abs(S_his);
Rh=abs(R_his);

Sh_mean=mean(Sh);
Sh_std=std(Sh);

Rh_mean=mean(Rh);
Rh_std=std(Rh);

TOTAL_RUNNING_TIME = etime(clock,TIME_OF_PROGRAM_START)