mimo_main.m

这是一个关于多输入多输出系统的频域忙识别系统的matlab源码

      for jj=1:m
         Hest_phase(:,ii,jj) = abs(reshape(Hest(ii,jj,:),NF,1)).*exp(i*reshape((mod(Phase_est_unwrap(ii,jj,:)+pi,2*pi)-pi),NF,1));
      end
   end
   
   if PLOT_PHASE
      figure(3); 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(4); clf;   %%% Plot the unwraped phase
      for ii=1:m
         for jj=1:n
            subplot(m,n,(ii-1)*n+jj);
            title(sprintf('Phase Est 1 and True Unwrap',ii,jj));grid;hold on;
            plot(reshape((Phase_est_1_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(5); clf;   %%% Plot the unwraped phase
      for ii=1:m
         for jj=1:m
            subplot(m,m,(ii-1)*m+jj);
            title(sprintf('Phase Est 2 and True Unwrap',ii,jj));grid;hold on;
            plot(reshape((Phase_est_2_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
   end %%% End of if PLOT_PHASE
   
   Phase_unwrap_his(:,:,:,iloop)=Phase_est_unwrap;
   
   if SYSTEM_REAL
      Phase_est(:,:,NF/2+2:NF)=-Phase_est(:,:,NF/2:-1:2);
   end
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%%% Phase Estimation Ends here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
   
   
   
   hest_L=hest(1:L,:,:);
      
   hest_L_his(:,:,:,iloop)=hest_L;
   
   figure(6);clf
   for ii=1:m
      for jj=1:n
         subplot(m,n,(ii-1)*n+jj);hold on;
         plot(reshape(abs(H(:,ii,jj)),NF,1),'rx');
         plot(reshape(abs(Hest(ii,H_order(jj),:)),NF,1),'b-');
         plot(EIG_INDEX',0,'r*');
         legend('True Mag','Est Mag','Selected Freq');
         title(sprintf('|H(%d,%d)|',ii,jj));
         Current_axis=axis;
         Current_axis(1)=1;
         Current_axis(2)=NF;
         axis(Current_axis);
         grid
      end
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Input Reconstruction Begin here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   if RECONSTRUCT_INPUT
      
      if Minimum_Phase_System
         HREC = fft(hest_L(:,:,H_order(1:n)),NF,1);
      else
         HREC = shiftdim(Hest(:,H_order(1:n),:),2);
      end
      
      epsilon=0.0001;
      Hinv=zeros(NF,n,m);
      
      for w=1:NF
         H_dummy = shiftdim(HREC(w,:,:),1);
         H_dummy = inv(H_dummy'*H_dummy + epsilon*eye(n,n)) * H_dummy';
         
         Hinv(w,:,:) = H_dummy;
      end
      
      %hinv = ifft(Hinv,NF,1);
      hinv = real(ifft(Hinv,NF,1));
      
      for ii=1:n
         for jj=1:m
            hinv(:,ii,jj)=fftshift(hinv(:,ii,jj));
         end
      end
      
      srec=zeros(n,N);
      for ii=1:n
         for jj=1:m
            srec(ii,:) = srec(ii,:) + filter(hinv(:,ii,jj)',1,x(jj,:));
         end
      end
      
      ALL=zeros(n,n,NF);
      for w=1:NF
         ALL(:,:,w)=reshape(Hinv(w,:,:),n,m)*reshape(H(w,:,1:n),m,n);
      end
      
      h_all = real(fftshift(ifft(shiftdim(ALL,2),NF,1)));
      
      %   h_all=real(h_all);
      
      figure(7);clf
      for ii=1:n
         for jj=1:n
            subplot(2*n,n,(ii-1)*2*n+jj);hold on;grid;
            plot(real(h_all(:,ii,jj)));
            if(ii==1)
               title(sprintf('Real part of the whole system impulse response \n(The whole system should be close to identity matrix)'));
            else
               title(sprintf('Real part of the whole system impulse response'));
            end
            
            subplot(2*n,n,(ii-1)*2*n+n+jj);hold on;grid;
            plot(imag(h_all(:,ii,jj)));
            title(sprintf('Imag part of All the system impulse response'));
         end
      end
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%% Cumulant Optimization code here
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      TSTART = clock;
      
      figure(8);clf;
      for ii=1:n
         [f, c, s_eq(ii,:)] = cum3equalizer(real(srec(ii,:)), Equalizer_Length);
         for jj=1:n
            s_eq_Xcorr_s=xcorr(s_eq(ii,:),s(jj,:),MAX_CORRELATION_LENGTH);
            if MATLAB_VERSION_NUMBER(1) >= '6'
                s_eq_Xcorr_s=fliplr(s_eq_Xcorr_s);
            end
            
            max_e_Xcorr_s(ii,jj)=max(abs(s_eq_Xcorr_s));
            subplot(n,n,(ii-1)*n+jj);hold on;grid;
            plot(s_eq_Xcorr_s);
            title(sprintf('s_{eq}(%d) Xcorr s(%d)',ii,jj));
         end
         [max_e_Xcorr_s_dummy e_s_index(ii,1)]=max(max_e_Xcorr_s(ii,:));
      end
      correct(iloop)=~sum(abs(sort(e_s_index)-(1:n).'))
      
      if ~correct(iloop)
         e_s_index(2)=3-e_s_index(1);
      end
      
      hest_CUM=zeros(MAX_Minus_ORDER+MAX_Plus_ORDER+1,m,n);
      
      for ii=1:m
         for jj=1:n
            f=xcorr(s_eq(jj,:),x(ii,:),MAX_CORRELATION_LENGTH)/N;
            if MATLAB_VERSION_NUMBER(1) >= '6'
                f=fliplr(f);
            end
            
            f_xcorr_h=xcorr(f,h(:,ii,e_s_index(jj)),MAX_CORRELATION_LENGTH);
            if MATLAB_VERSION_NUMBER(1) >= '6'
                f_xcorr_h=fliplr(f_xcorr_h);
            end
            
            [f_xcorr_h_max f_xcorr_h_max_index]=max(abs(f_xcorr_h));
            
            min_f_index=MAX_CORRELATION_LENGTH-f_xcorr_h_max_index-MAX_Minus_ORDER+2;
            if min_f_index >= 1
               hest_CUM(:,ii,e_s_index(jj))=f(min_f_index+(0:MAX_Minus_ORDER+MAX_Plus_ORDER)).';
            else
               hest_CUM(1:(1-min_f_index),ii,e_s_index(jj))=0;
               hest_CUM((2-min_f_index:MAX_Minus_ORDER+MAX_Plus_ORDER+1),ii,e_s_index(jj))=...
                  f(1:MAX_Minus_ORDER+MAX_Plus_ORDER+min_f_index).';
            end
         end
      end
            
      
      hest_CUM_his(:,:,:,iloop)=hest_CUM;
      
      TIME_Optimization = etime(clock,TSTART)

      
      figure(9);clf;
      for ii=1:m
         for jj=1:n
            subplot(m,n,(ii-1)*n+jj);hold on;
            
            [h_abs_max h_abs_max_index]=max(abs(h(:,ii,jj)));
            f_abs_max_index=h_abs_max_index+MAX_Minus_ORDER;
            fh_coeff=sign(h(h_abs_max_index,ii,jj)) * sign(hest_CUM(f_abs_max_index,ii,jj)) * sqrt(sum(h(:,ii,jj).^2)/sum(hest_CUM(:,ii,jj).^2));
            plot(hest_CUM(:,ii,jj)*fh_coeff, 'b-*');
            stem(MAX_Minus_ORDER+(1:L),h(:,ii,jj),'r');
            legend('Estimated IR','True IR');
            title('Impulse Response');
            grid;
         end
      end
      
      drawnow;
      
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      %%% Cumulant Optimization code ENDS here
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      
      h_all_his(:,:,:,iloop)=h_all;
      
   end  %% End of if RECONSTRUCT_INPUT
   
   Hest_his(:,:,:,iloop)=Hest;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Input Reconstruction End here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   drawnow
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  End of Monte Carlo test here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end   % End of Monte_carlo Test iloop


hest_L_mean=mean(hest_L_his,4);
hest_L_std=std(hest_L_his,0,4);

Hest_L_his=fft(hest_L_his,NF,1);
Hest_L_abs_mean=mean(abs(Hest_L_his),4);
Hest_L_abs_std=std(abs(Hest_L_his),0,4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% In the following, we calculate the hest_ifft directly.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for iloop=1:RUN_TIMES
   for ii=1:m
      for jj=1:m
         hest_ifft_his(:,ii,jj,iloop)=real(fftshift(ifft(shiftdim(Hest_his(ii,jj,:,iloop),2))));
      end
   end
end

for iloop=1:RUN_TIMES
   for ii=1:m
      for jj=1:n
         hest_h_corr=xcorr(hest_ifft_his(:,ii,H_order(jj),iloop),h(:,ii,jj),MAX_CORRELATION_LENGTH_ifft);
         
         if MATLAB_VERSION_NUMBER(1) >= '6'
             hest_h_corr=flipud(hest_h_corr);
         end
         
         [max_xcorr_abs max_xcorr_position]=max(abs(hest_h_corr));
         %%% hest_ifft_Le_his is not in the correct column order, so H_order is needed.
         h_shift=MAX_CORRELATION_LENGTH_ifft+1-max_xcorr_position;
         hest_ifft_Le=hest_ifft_his(mod(h_shift-MAX_Minus_ORDER_ifft+(1:Le)-1,NF)+1,ii,H_order(jj),iloop)*...
            sign(hest_h_corr(max_xcorr_position));
         hest_ifft_Le=hest_ifft_Le*sqrt(sum(h(:,ii,jj).^2)/sum(hest_ifft_Le.^2));
         hest_ifft_Le_his(:,ii,H_order(jj),iloop)=hest_ifft_Le;
      end
   end
end   

hest_ifft_Le_mean = mean(hest_ifft_Le_his,4);
hest_ifft_Le_std  = std(hest_ifft_Le_his,[],4);

figure(10);clf;
   for ii=1:m
      for jj=1:n
         subplot(m,n,(ii-1)*n+jj);grid;hold on;
         title(sprintf('hest ifft(%d,%d)',ii,jj));
         h_long=zeros(Le,1);
         h_long((MAX_Minus_ORDER_ifft+(1:L)))=h(:,ii,jj);
         shadow_plot(h_long,hest_ifft_Le_mean(:,ii,H_order(jj)),hest_ifft_Le_std(:,ii,H_order(jj)));
      end
   end

for ii=1:m
   for jj=1:n
      h_long=zeros(Le,1);
      h_long((MAX_Minus_ORDER_ifft+(1:L)))=h(:,ii,jj);
      
      NMSE_ifft(ii,jj)=sum(sum((reshape(hest_ifft_Le_his(:,ii,H_order(jj),:),Le,RUN_TIMES)-...
         h_long*ones(1,RUN_TIMES)).^2)) / (h_long'*h_long) / RUN_TIMES;
   end
end


NMSE_ifft

ONMSE_ifft=sum(sum(NMSE_ifft)) / (m*n)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Calculate NMSE for the h_ifft END here
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   

if 1

hest_CUM_his_bak=hest_CUM_his;
   
if Process_hest_CUM_his
   for ii=1:m
      for jj=1:n
         [h_abs_max h_abs_max_index]=max(abs(h(:,ii,jj)));
         f_abs_max_index=h_abs_max_index+MAX_Minus_ORDER;
         for run=1:RUN_TIMES
            fh_coeff=sign(h(h_abs_max_index,ii,jj)) * sign(hest_CUM_his(f_abs_max_index,ii,jj,run)) * sqrt(sum(h(:,ii,jj).^2)/sum(hest_CUM_his(:,ii,jj,run).^2));
            %fh_coeff=sign(h(h_abs_max_index,ii,jj)) * sign(hest_CUM_his(f_abs_max_index,ii,jj,run));
            hest_CUM_his(:,ii,jj,run)=hest_CUM_his(:,ii,jj,run)*fh_coeff;
         end
      end
   end
end

hest_CUM_mean=mean(hest_CUM_his,4);
hest_CUM_std=std(hest_CUM_his,0,4);

figure(11);clf;
for ii=1:m
   for jj=1:n
      subplot(m,n,(ii-1)*n+jj);hold on;
      plot(reshape(hest_CUM_his(:,ii,jj,:),MAX_Minus_ORDER+MAX_Plus_ORDER+1,RUN_TIMES));
      stem(MAX_Minus_ORDER+(1:L),h(:,ii,jj),'r');
      title('Impulse Response h_{CUM}');
      grid;
   end
end


for ii=1:m
   for jj=1:n
      h_long=zeros(MAX_Minus_ORDER+MAX_Plus_ORDER+1,1);
      h_long((MAX_Minus_ORDER+(1:L)))=h(:,ii,jj);
      
      NMSE_CUM(ii,jj)=sum(sum((reshape(hest_CUM_his(:,ii,jj,:),MAX_Minus_ORDER+MAX_Plus_ORDER+1,RUN_TIMES)-...
         h_long*ones(1,RUN_TIMES)).^2)) / (h_long'*h_long) / RUN_TIMES;
   end
end
   
   
NMSE_CUM

ONMSE_CUM=sum(sum(NMSE_CUM)) / (m*n)


end
% clear B_* cx* Cy* Ry* Rx* V* s x noise DING* Hest_L_his W* A A1 S_x* cum* Ph*
Total_TIME = etime(clock,Program_START)