equalizer_mmse.m

采用MMSE准则设计最佳均衡器的源程序

% Equalizer design based on the MMSE cost function

clear all;
close all;

% Modulation order
M =4;
bps =log2(M);

% N path tap delay line channel (assuming ISI from N-1 previous symbols)
N=5;

% Your range of SNR
EbNodB=[0:1:25];
[a,datapoints]=size(EbNodB);
EbNo=10.^(EbNodB/10);
EsNo=bps*EbNo;

%The noise variance (this is now a complex random variable)
var=1./(EsNo);

%Bits per frame
NumberBits=2000;
NumberSymbols=NumberBits/bps;

% Number of frame errors per SNR (since these are Monte Carlo trials, the  more errors you log, the more accurate your simulated results)
err_vector=[75*ones(1,datapoints)];
max_errors = err_vector;
minBER = 10^-5;

% initialize storage
frame_errors = zeros( 1,datapoints);
bit_errors = zeros(1,  datapoints);
trials = zeros(1, datapoints );
P_e = zeros(1, datapoints);
P_b = zeros(1,  datapoints);

fname = strcat( 'MMSE_EQ','Paths',sprintf('%d', N),'.mat' );
% % % %   see if this is a continuation (this allows you to start your % simulation from where it was halted during the previous run)
%      if (fopen(fname, 'r') > 0)
%             load (fname);
%         end
% reseed the random number generator
 randn( 'state', cputime); 

 
 mapper=de2bi(bitxor([0:1:M-1],floor([0:1:M-1]/2)), log2(M),'left-msb'); 
 symbols=sqrt(0.5)*[-1-sqrt(-1),-1+sqrt(-1),1+sqrt(-1),1-sqrt(-1)];  
 [rows, cols]=size(mapper);
 for snrpoint = 1:datapoints
    fprintf( '\nEbNo = %f dB\n', EbNodB(snrpoint) );
        
    while  frame_errors( snrpoint) < max_errors( snrpoint ) 
        trials(snrpoint) = trials(snrpoint) + 1 ;
       
        % Generate random data bits
        cin=round(rand(1, NumberBits)); 
         k=1;
          for i=1:bps:NumberBits
              ind=cin(i:i+bps-1);
             for j =1:rows
               if ((ind(1,1:bps)==mapper(j,1:bps)))
                 sym =j;
                 break;        
               end
             end
              tak(k) =symbols(sym);  
              k=k+1;
          end
          
          pad =N-1;
          H=zeros(2*(pad+1), (2*(pad+1)-1));
          i=0;
          h1= randn(1,N)+sqrt(-1)*randn(1,N);
          h2= randn(1,N)+sqrt(-1)*randn(1,N); 
          for r=1:N
             H(r+i,r:(r-1)+N)=(h1)/sqrt(sum(abs(h1).^2));
             H(r+i+1,r:(r-1)+N)=(h2)/sqrt(sum(abs(h2).^2));
             i=i+1;
           end   
    
           d_hat=zeros(1,NumberSymbols);    
          for k=pad+1:NumberSymbols-pad
              
              x_temp =[tak(k+pad:-1:k-pad)];  
              n=sqrt(var(snrpoint)/2)*(randn(2*N,1)+sqrt(-1)*randn(2*N,1));        
              y=H*x_temp.'+n;              
              % Equalizer taps              
              filter_taps= inv(H*H'+ var(snrpoint)*eye(2*N))*H(:,N);
              %Equalizer outputs   
              d_hat(k)=filter_taps'*y;              
          end
          cout_r=(real(d_hat)>0);
          cout_i=(imag(d_hat)>0);
          
          cout(1:bps:NumberBits)=cout_r;
          cout(2:bps:NumberBits)=cout_i;
          
         errors=0;
         [err,vec]=find(cin((pad*bps)+1:(NumberSymbols-pad)*bps)~=cout((pad*bps)+1:(NumberSymbols-pad)*bps));
         errors=sum(err); 
  
        frame_errors(snrpoint) = frame_errors( snrpoint ) + (errors >0);
        bit_errors(snrpoint ) =   bit_errors(snrpoint ) + errors;       
        P_e(snrpoint) = frame_errors(snrpoint)/trials( snrpoint );% Frame error rate
        P_b(snrpoint) = bit_errors(snrpoint)/( trials(snrpoint)*(NumberBits-pad*bps) );% Bit error rate        
        save (fname, 'frame_errors','bit_errors','trials', 'P_e', 'P_b','EbNodB');
    end% trials    
        fprintf('%f  %f   \n',EbNodB(1,snrpoint),P_b(1,snrpoint) );
       if (P_b(snrpoint) <1e-5)
           break;
       end

end%snrpoint

 
 
semilogy(EbNodB, P_b);
grid on;