mimo_main.m

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

      for ii=1:m
         %s(ii,:)=qam(N,2,465*sum(100*clock),54.564*sum(100*clock));
      end
      
      %s = rand(m,N);
      %s = s-(mean(s'))'*ones(1,N);
      %s=sign(s);
   end
   
   s = s-(mean(s'))'*ones(1,N);
   s=diag(1./std(s,0,2))*s;       %% Input Signal
   if Less_Input_than_Output
      s(n+1:m,:)=0;
   end
   %%% The input signals are generated by now.


   %%% Generate the received signal
   x = zeros(size(s));          %%% Observed mixture , received signal
   for ii=1:m
      for jj=1:n
         x(ii,:) = x(ii,:) + filter(h(:,ii,jj),1,s(jj,:));
      end               
   end
   x = x-(mean(x'))'*ones(1,N);
   
   std_x_1=std(x(1,:));
   std_x_2=std(x(2,:));
   
   if ADD_NOISE
      for ii=1:m
         rand('state',sum(100*clock)*rand);
         noise=randn(1,N);
         noise=noise-mean(noise);
         power_coefficient = (10^(SNR/20))*std(noise)/std(x(ii,:));
         x(ii,:) = x(ii,:)+noise/power_coefficient;
      end
   end
   
   x(1,:) = std_x_1*x(1,:)/std(x(1,:));
   x(2,:) = std_x_2*x(2,:)/std(x(2,:));
   %%% The received signal are generated by now.
   
   
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   % The Process Begin here
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   %%% Estimating the cross correlation
   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
         
         if MATLAB_VERSION_NUMBER(1) >= '6'
             rx(:,ii,jj)=flipud(rx(:,ii,jj));
         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
   S_x=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);
         S_x(:,ii,jj)=fft(fftshift(rr),NF);
      end
   end
   
   if S_x_symmetry
      for w=1:NF
         S_x_dummy=reshape((S_x(w,:,:)+S_x(w,:,:))/2,m,m);
         S_x(w,:,:)=(S_x_dummy+S_x_dummy')/2;
      end
   end
   
   
   %% Estimating the cross cumulants
   TSTART = clock;
   if w2_equal_minus_w3_plus_arfa & ~Using_Bispectrum
      cx = zeros(2*C_LENGTH+1,4*C_LENGTH+1,m,m,m);
   else % w3=0
      cx = zeros(2*C_LENGTH+1,2*C_LENGTH+1,m,m,m);
   end
   
   if ~Using_Bispectrum
      cx_phase = zeros(2*C_LENGTH+1,2*C_LENGTH+1,m,m,m);
   end
   
   if ADD_CUM_WINDOW
      for kk=1:2*C_LENGTH+1
         cum_4_window_4(:,:,kk)=ones(2*C_LENGTH+1,2*C_LENGTH+1)*dd_win(kk);
      end
   end
   
   for ii=1:m
      for jj=1:m
         for kk=1:m
            for seg=0:seg_num-1
               if Using_Bispectrum
                  cx(:,:,ii,jj,kk)=cx(:,:,ii,jj,kk)+g_cc3_matrix(...
                     x(ii,seg_length*seg+1:seg_length*(seg+1)),...
                     x(jj,seg_length*seg+1:seg_length*(seg+1)),...
                     x(kk,seg_length*seg+1:seg_length*(seg+1)),C_LENGTH)/seg_num;
               else  % Using Trispectrum
                  
                  %TSTART_matrix = clock;
                  cx_3D = g_cc4_3D(...
                     x(ii,seg_length*seg+1:seg_length*(seg+1)),...
                     conj(x(jj,seg_length*seg+1:seg_length*(seg+1))),...
                     x(kk,seg_length*seg+1:seg_length*(seg+1)),...
                     conj(x(kk,seg_length*seg+1:seg_length*(seg+1))),C_LENGTH);
                  %TIME_matrix = etime(clock,TSTART_matrix)
                  
                  %TSTART_mex = clock;
                  %cx_3D = g_ccum4_3D_mex(...  %% Using the mex function for fast speed
                  %x(ii,seg_length*seg+1:seg_length*(seg+1)).',...
                  %conj(x(jj,seg_length*seg+1:seg_length*(seg+1))).',...
                  %x(kk,seg_length*seg+1:seg_length*(seg+1)).',...
                  %conj(x(kk,seg_length*seg+1:seg_length*(seg+1))).',C_LENGTH);
                  %TIME_mex = etime(clock,TSTART_mex)
                  
                  
                  cx_phase_3D = g_cc4_3D(...
                     x(ii,seg_length*seg+1:seg_length*(seg+1)),...
                     conj(x(jj,seg_length*seg+1:seg_length*(seg+1))),...
                     x(kk,seg_length*seg+1:seg_length*(seg+1)),...
                     conj(x(ii,seg_length*seg+1:seg_length*(seg+1))),C_LENGTH);
                  
                  if w2_equal_minus_w3_plus_arfa
                     if ADD_CUM_WINDOW
                        cx_3D=cx_3D.*cum_4_window_4;
                     end
                     
                     cx_3D=cx_3D.*arfa_matrix;
                     
                     cx_dummy=zeros(2*C_LENGTH+1,4*C_LENGTH+1);
                     for dd=-2*C_LENGTH:2*C_LENGTH
                        for t1=1:2*C_LENGTH+1
                           cx_dummy(t1,dd+2*C_LENGTH+1)=sum(diag(shiftdim(cx_3D(t1,:,:),1),-dd));
                        end
                     end
                     cx(:,:,ii,jj,kk) = cx(:,:,ii,jj,kk) + cx_dummy/seg_num;
                  else  % w3_equal_0
                     if ADD_CUM_WINDOW
                        cx(:,:,ii,jj,kk) = cx(:,:,ii,jj,kk) + sum(cx_3D.*cum_4_window_4,3)/seg_num;
                     else
                        cx(:,:,ii,jj,kk) = cx(:,:,ii,jj,kk) + sum(cx_3D,3)/seg_num;
                     end
                  end
                  
                  if ADD_CUM_WINDOW
                     cx_phase(:,:,ii,jj,kk) = cx_phase(:,:,ii,jj,kk) + sum(cx_phase_3D.*cum_4_window_4,3)/seg_num;
                  else
                     cx_phase(:,:,ii,jj,kk) = cx_phase(:,:,ii,jj,kk) + sum(cx_phase_3D,3)/seg_num;
                  end
                  
               end
            end
         end
      end
   end
   
   TIME_cx = etime(clock,TSTART)
   
   if CUMULANTS_ARE_REAL
      cx=real(cx);
      
      if ~Using_Bispectrum
         cx_phase=real(cx_phase);
      end
   end
   
   
   if ADD_CUM_WINDOW 
      for ii=1:m
         for jj=1:m
            for kk=1:m
               cx(:,:,ii,jj,kk) = cx(:,:,ii,jj,kk).*cum_win;
            end
         end
      end
   end
   
   
   %% Estimating the cross polyspectra matrix for magnitude estimation
   TSTART = clock;
   B_x=zeros(NF,1,m,m,m);
   B_xn=zeros(NF,Freq_number,m,m,m);
   if w2_equal_minus_w3_plus_arfa & ~Using_Bispectrum
      for ii=1:m
         for jj=1:m
            for kk=1:m
               cx_dummy=zeros(NF);
               cx_dummy((NF/2+1-C_LENGTH) : (NF/2+1+C_LENGTH) , (NF/2+1-2*C_LENGTH) : (NF/2+1+2*C_LENGTH))=cx(:,:,ii,jj,kk);
               cx_dummy=fftshift(cx_dummy);
               B_x_dummy=(fft2(cx_dummy,NF,NF));
               B_x(:,1,ii,jj,kk)=B_x_dummy(:,w2);
               B_xn(:,1:Freq_number,ii,jj,kk)=B_x_dummy(:,Ref_Frequencies);
            end
         end
      end
   else % Using Bispectrum or w3=0
      for ii=1:m
         for jj=1:m
            for kk=1:m
               %%% The following is the diagonal slices for (w, \beta-w)
               cx_dummy=zeros(NF);
               cx_dummy((NF/2+1-C_LENGTH) : (NF/2+1+C_LENGTH) , (NF/2+1-C_LENGTH) : (NF/2+1+C_LENGTH))=cx(:,:,ii,jj,kk);
               cx_dummy=fftshift(cx_dummy);
               B_x_dummy=(fft2(cx_dummy,NF,NF));
               all_w=1:NF;
               B_x(:,1,ii,jj,kk)=diag(B_x_dummy(all_w,mod(sum_w1_w2+1-all_w, NF)+1));
               
               for Freq_Index=1:Freq_number
                  B_xn(:,Freq_Index,ii,jj,kk)=diag(B_x_dummy(all_w,mod(Ref_Frequencies(Freq_Index)+1-all_w,NF)+1));
               end
            end
         end
      end
   end
   
   TIME_B_x = etime(clock,TSTART)
   
   
   %% Estimating the cross polyspectra matrix for phase retrieval
   if Using_Bispectrum
      B_x_phase=zeros(NF,m,m,m);
      for ii=1:m
         for jj=1:m
            for kk=1:m
               cx_dummy=zeros(NF);
               cx_dummy((NF/2+1-C_LENGTH) : (NF/2+1+C_LENGTH) , (NF/2+1-C_LENGTH) : (NF/2+1+C_LENGTH))=cx(:,:,ii,jj,kk);
               cx_dummy=fftshift(cx_dummy);
               B_x_dummy=(fft2(cx_dummy,NF,NF));
               
               B_x_phase(NF,ii,jj,kk)=B_x_dummy(2,mod(k_arfa-1, NF)+1);
               B_x_dummy(2:NF-1,:)=B_x_dummy(NF:-1:3,:);
               B_x_phase(1:NF-k_arfa,ii,jj,kk)=diag(B_x_dummy(1:NF-k_arfa, mod(k_arfa,NF)+1:NF));
               if k_arfa > 1
                  B_x_phase(NF-k_arfa+1:NF-1,ii,jj,kk)=diag(B_x_dummy(NF-k_arfa+1:NF-1, 1:k_arfa-1));
               end % if k_arfa > 1
            end % kk
         end % jj
      end % ii
   else % Using Trispectrum
      B_x_phase=zeros(NF,m,m,m);
      for ii=1:m
         for jj=1:m
            for kk=1:m
               cx_dummy=zeros(NF);
               cx_dummy((NF/2+1-C_LENGTH) : (NF/2+1+C_LENGTH) , (NF/2+1-C_LENGTH) : (NF/2+1+C_LENGTH))=cx_phase(:,:,ii,jj,kk);
               cx_dummy=fftshift(cx_dummy);
               B_x_dummy=(fft2(cx_dummy,NF,NF));
               
               B_x_phase(NF,ii,jj,kk)=B_x_dummy(2,mod(k_arfa-1, NF)+1);
               B_x_dummy(2:NF-1,:)=B_x_dummy(NF:-1:3,:);
               B_x_phase(1:NF-k_arfa,ii,jj,kk)=diag(B_x_dummy(1:NF-k_arfa, mod(k_arfa,NF)+1:NF));
               if k_arfa > 1
                  B_x_phase(NF-k_arfa+1:NF-1,ii,jj,kk)=diag(B_x_dummy(NF-k_arfa+1:NF-1, 1:k_arfa-1));
               end
            end
         end
      end
   end
   
   Rx=reshape(shiftdim(S_x(:,:,:),1),m,m,NF);    %% The cross power spectrum matrix
   Rx1=reshape(shiftdim(B_x(:,1,Rx1_index,:,:),2),m,m,NF);  %% The cross polyspectra matrix
   
   
   for Freq_Index=1:Freq_number  %% Rxn is a "m x m x NF x Freq_number" matrix, the last m is the index ii
      Rxn(:,:,:,Freq_Index)=reshape(shiftdim(B_xn(:,Freq_Index,Rx1_index,:,:),3),m,m,NF);  %% The cross polyspectra matrix
   end
   
   jj=0;
   for ii=1:m
      if ii ~= Rx1_index
         jj=jj+1;
         for Freq_Index=1:Freq_number
            Rxn(:,:,:,Freq_Index+jj*Freq_number)=reshape(shiftdim(B_xn(:,Freq_Index,ii,:,:),3),m,m,NF);
         end
      end
   end
   
   
   if Rx_diagonal_real
      for w=1:NF
         Rx(:,:,w) = Rx(:,:,w)-i*imag(diag(diag(Rx(:,:,w))));
      end
   end
   
   %%% Calculate the condition number of the cross-power spectrum matrix
   for w=1:NF
      Rx_eig(:,w)=sort(abs(real(svd(Rx(:,:,w)))));
      Rx_cond(w,1)=real(Rx_eig(m,w)/Rx_eig(1+m-n,w));
   end

   [Rx_cond_sort Rx_cond_sort_index]=sort(Rx_cond);
   
   Rx_cond_max=Rx_cond(Rx_cond_sort_index(floor(NF*Rx_cond_selection_percentage/100)));
   
   Rx_cond_Index=[]; %% Select freqs by the conditiona number of Rx(w)
   for w=1:NF
      if Rx_cond(w,1) < Rx_cond_max
         Rx_cond_Index=[Rx_cond_Index w];
      end
   end
   
   Rx_INDEX=zeros(1,NF);
   Rx_INDEX(Rx_cond_Index)=1;
   
   %%% Calculate V(w) and V^{-1}(w)
   if Less_Input_than_Output
      for w=1:NF
         [Urx Drx Vrx]=svd(Rx(:,:,w));
         Dv=diag(Drx);
         Dv(1+n:m)=0.000001;
         Dv_1=sqrt(Dv);
         Dv(1:n)=1./sqrt(Dv(1:n));
         V(:,:,w) = Urx*diag(Dv)*Vrx';
         V_1(:,:,w) = Urx*diag(Dv_1)*Vrx';
      end
   else
      for w=1:NF
         V_1(:,:,w) = sqrtm(Rx(:,:,w));
         if Rx_cond(w,1) < Rx_cond_max |1
            V(:,:,w) = inv(V_1(:,:,w));
         else
            [Urx Drx Vrx]=svd(Rx(:,:,w));
            Dv=diag(Drx);
            Dv(1:m-1)=1./sqrt(Dv(1:m-1));
            Dv(m)=0.000001;
            V(:,:,w) = Urx*diag(Dv)*Vrx';
         end
      end
   end
   
   if V_diagonal_real
      for w=1:NF
         V(:,:,w) = V(:,:,w)-i*imag(diag(diag(V(:,:,w))));
      end
   end
   
   
   %%% Interpolate V and V^{-1} in frequencies with high condition numbers
   if interpolate_V
      for ii=1:m
         for jj=1:m
            V_inter_real(:,ii,jj)=interpolate_func(reshape(real(V(ii,jj,:)),NF,1),Rx_cond_Index);
            V_inter_imag(:,ii,jj)=interpolate_func(reshape(imag(V(ii,jj,:)),NF,1),Rx_cond_Index);
         end
      end
      
      V_inter=shiftdim(V_inter_real+i*V_inter_imag,1);
      
      V=V_inter;
   end %% End of if interpolate_V