stbc_ofdm_est_ls.asv

是关于空时编码的nrs

%==========================================================================
% "Channel estimation via training in time domain and tracking in frequcy 
% domain for STBC-OFDM Systems"
% Sun Jie-ting, 08.11.2004, East China Uni. Science and Technology,Shanghai
% Results of SER
%==========================================================================

%clc;
function data=STBC_OFDM_EST_LS(Number_Relz,NoOFDMSymbol,NFFT,G,t_a,training_gap,f_dmax,snr_min,snr_max)

Nt = 2;
Nr = 1;
%NFFT = 128;               % FFT length
%G = 10;                   % Guard interval length

%t_a = 0.8*10^(-6);       % Sampling duration of wireless communication
%miu = 0.05;
%training_gap = 10;
%NoOFDMSymbol = 40;        %Number of OFDM symbol
OFDMSymbolbits = NFFT;      %Number of bits per OFDM symbol conveys
length_data = (NoOFDMSymbol-NoOFDMSymbol/training_gap)*NFFT;  
                            % The total data length

%-------------------------------------------------
% Parameters for Monte Carlo channel
%-------------------------------------------------



symbol_duration = NFFT * t_a;   %OFDM symbol duration
number_of_summations = 40;      % Number of summations for Monte-Carlo method



%f_dmax = 50.0;                  % Maximum Doppler frequency

load h_decimation.am -ascii;
h1_initial = h_decimation;

h2_initial = h_decimation;

N_P  = length(h_decimation);
 
%Number_Relz = 100;
step = 1;
symbol_error_rate = zeros(1,(snr_max-snr_min)/step+1);
for number_of_realization= 1: Number_Relz;

disp('number_of_realization');
disp(number_of_realization);
u1 = rand(N_P,number_of_summations); % A random variable
u2 = rand(N_P,number_of_summations); % A random variable
%---------------------------
% Data for STBC demodulation
%---------------------------
data_for_destbc = zeros(2,1,16);
data_qpsk = [1+j;-1+j;-1-j;1-j]/sqrt(2);
for k=1:4
    for m=1:4
        data_for_destbc(:,:,4*(k-1)+m) = [data_qpsk(k);data_qpsk(m)];
    end 
end
%--------------------------
% Optimal Training Sequence
%--------------------------
Training_OFDM_T1 = zeros(1,OFDMSymbolbits);
Training_OFDM_T2 = zeros(1,OFDMSymbolbits);
for i=1:OFDMSymbolbits
    if mod(i,4)==0
        Training_OFDM_T1(i) = data_qpsk(4);
    else
        Training_OFDM_T1(i) = data_qpsk(mod(i,4));
    end
    Training_OFDM_T2(i) = Training_OFDM_T1(i)*exp(j*2*pi*(i-1)*(N_P)/NFFT);
end

%----------------------------
% QPSK modulator in base band
%----------------------------------------------
BitFrameC1=(2*round(rand(Nt,OFDMSymbolbits/2,NoOFDMSymbol))-1)/sqrt(Nt);
QpskDataC1=BitFrameC1(1,:,:)+j*BitFrameC1(2,:,:);

BitFrameC2=(2*round(rand(Nt,OFDMSymbolbits/2,NoOFDMSymbol))-1)/sqrt(Nt);
QpskDataC2=BitFrameC2(1,:,:)+j*BitFrameC2(2,:,:);
%--------------------
% Alamouti space-time modulator in base band
%--------------------
Stbc1=zeros(1,OFDMSymbolbits,NoOFDMSymbol);
Stbc2=zeros(1,OFDMSymbolbits,NoOFDMSymbol);
for k=1:1:OFDMSymbolbits/2
    Stbc1(1,2*k-1,:)=QpskDataC1(1,k,:);
    Stbc1(1,2*k,:)=-conj(QpskDataC2(1,k,:));
    Stbc2(1,2*k-1,:)=QpskDataC2(1,k,:);
    Stbc2(1,2*k,:)=conj(QpskDataC1(1,k,:));
end

%------------------------------------------------------------------------
% Transmitted signal of trasmitt antena 1 to receive antenna (channel:h1)
%------------------------------------------------------------------------
%snr_min =12;
%snr_max =12;


for snr = snr_min:step:snr_max; 
    disp('snr = ');
    disp(snr);


h1_frame = [];
rs1MP_frame = [];
initial_time=0;                 % Initial time

for i=1:NoOFDMSymbol;

   if mod(i,training_gap) == 0  | i==1
%       disp('Insert Training Sequence');
       OFDM_signal_tem = OFDM_Modulator(Training_OFDM_T1,NFFT,G);
   else 
       OFDM_signal_tem = OFDM_Modulator(Stbc1(1,:,i),NFFT,G);
   end
   % OFDM signal from the first transmitt antenna is created
   
  [h1, t] = MCM_channel_model(u1, initial_time, number_of_summations, symbol_duration, ...,
       f_dmax, h1_initial);
  h1_frame = [h1_frame;h1];
  initial_time = t;
   
   rs1MP = conv(OFDM_signal_tem, h1);
%   H1_i = fft([h1 zeros(1,NFFT-N_P)])
%   pause
%   r3 = Stbc1(1,:,i).*H1_i
%   pause
   % The received signal over multhipath channel is created
%   rs1 = awgn(rs1,snr,'measured','dB');
   rs1MP_frame = [rs1MP_frame; rs1MP];
   clear OFDM_signal_tem;
end;


%------------------------------------------------------------------------
% Transmitted signal of antenna 2 to receive antenna 1 (channel: h2)
%-----------------------------------------------------------------------

h2_frame = [];
rs2MP_frame = [];

initial_time=0;                 % Initial time
for i=1:NoOFDMSymbol;

   if mod(i,training_gap) == 0 | i==1
%       disp('Insert Training Sequence');
       OFDM_signal_tem = OFDM_Modulator(Training_OFDM_T2,NFFT,G);
   else 
       OFDM_signal_tem = OFDM_Modulator(Stbc2(1,:,i),NFFT,G);
   end
   

   % OFDM signal from the second antenna is created
   [h2, t] = MCM_channel_model(u2, initial_time, number_of_summations, symbol_duration, ...,
       f_dmax, h2_initial);
   h2_frame = [h2_frame;h2];
   initial_time = t;
%   H2_i = fft([h2 zeros(1,NFFT-N_P)])
%   pause
%   r4 = Stbc2(1,:,i).*H2_i
%   pause   
   rs2MP = conv(OFDM_signal_tem, h2);
   % The received signal over multhipath channel is created
%   rs2 = awgn(rs2,snr,'measured','dB');
   rs2MP_frame = [rs2MP_frame; rs2MP];
   clear OFDM_signal_tem;
end;


%-----------------
% Recever: OFDM demodulator, channel estimation
%-----------------


for i=1:NoOFDMSymbol;
   rsBeforeAwgn_i = rs1MP_frame(i,:) + rs2MP_frame(i,:);
   rs_i = awgn(rsBeforeAwgn_i,snr,'measured','dB');
   Demodulated_rs_i = OFDM_Demodulator(rs_i,NFFT,NFFT,G);

    if mod(i,training_gap)==0 | i==1
%       disp('Training Estimation');
       p1 = zeros(1,N_P);
       p2 = zeros(1,N_P);
       for l = 1:N_P
           for k = 1:NFFT
               p1(l) = p1(l) + Demodulated_rs_i(k)*conj(Training_OFDM_T1(k))*exp(j*2*pi*(k-1)*(l-1)/NFFT);
               p2(l) = p2(l) + Demodulated_rs_i(k)*conj(Training_OFDM_T2(k))*exp(j*2*pi*(k-1)*(l-1)/NFFT);               
           end
       end
       h1_estimation = 1/NFFT*p1;
       h2_estimation = 1/NFFT*p2;
       H1_i_estimation = fft([h1_estimation zeros(1,NFFT-N_P)]);
       H2_i_estimation = fft([h2_estimation zeros(1,NFFT-N_P)]);

       H1_i = fft([h1_frame(i,:) zeros(1,NFFT-N_P)]);
       H2_i = fft([h2_frame(i,:) zeros(1,NFFT-N_P)]);
    else

%       H1_i = fft([h1_frame(i,:) zeros(1,NFFT-N_P)]);
%       H2_i = fft([h2_frame(i,:) zeros(1,NFFT-N_P)]);
        H1_i = H1_i_estimation;
        H2_i = H2_i_estimation;
 %--------------------------
 % Demodulated signal 
 %--------------------------
        estimated1_OFDM = [];
        estimated2_OFDM = [];        
        for l=1:OFDMSymbolbits/2;
           R1=Demodulated_rs_i(2*l-1);
           R2=Demodulated_rs_i(2*l);
           R2conj = conj(R2);
           R_n2l=[R1;R2conj];
           H_n2l = [H1_i(2*l-1) H2_i(2*l-1);conj(H2_i(2*l)) -conj(H1_i(2*l))];
           minNorm = 10000;
           num_min = 0;
           for m=1:16
               dec_norm = norm(R_n2l-H_n2l*data_for_destbc(:,:,m))^2;
               if dec_norm < minNorm;
                   minNorm = dec_norm;
                   num_min = m;
               end
            end
            deStbc = data_for_destbc(:,:,num_min)
            estimated1_OFDM = [estimated1_OFDM deStbc(1) -conj(deStbc(2))];
            estimated2_OFDM = [estimated2_OFDM deStbc(2) conj(deStbc(1))];            
            TrSTBC = [Stbc1(1,2*l-1,i);Stbc2(1,2*l-1,i)]
            if deStbc==TrSTBC
%                disp('yes');
            else
                symbol_error_rate(snr-snr_min+1) = symbol_error_rate(snr-snr_min+1) + 1;
%                disp('no');
            end
        end;
        P1 = zeros(1,N_P);
        P2 = zeros(1,N_P);            
        for l = 1:N_P
            for k = 1:NFFT
                P1(l) = P1(l) + Demodulated_rs_i(k)*conj(estimated1_OFDM(k))*exp(j*2*pi*(k-1)*(l-1)/NFFT);
                P2(l) = P2(l) + Demodulated_rs_i(k)*conj(estimated2_OFDM(k))*exp(j*2*pi*(k-1)*(l-1)/NFFT);               
             end
         end            
         Q12 = zeros(N_P,N_P);
         Q21 = zeros(N_P,N_P);
         for l = 1:N_P
             for k = 1:N_P
                 Q12(l,k) = sum(estimated1_OFDM.*conj(estimated2_OFDM).*exp(j*2*pi*[0:1:NFFT-1]*(k-l)/NFFT));
                 Q21(l,k) = sum(estimated2_OFDM.*conj(estimated1_OFDM).*exp(j*2*pi*[0:1:NFFT-1]*(k-l)/NFFT));
             end
         end            
         h1_estimation = 1/NFFT*(P1-h2_estimation*Q21);
         h2_estimation = 1/NFFT*(P2-h1_estimation*Q12);
         H1_i_estimation = fft([h1_estimation zeros(1,NFFT-N_P)]);
         H2_i_estimation = fft([h2_estimation zeros(1,NFFT-N_P)]);          
    end;
end;    %NoOFDMSymbol
end;    %snr
end;    %realization
ser = symbol_error_rate/(length_data*Number_Relz);
%semilogy(snr_min:step:snr_max, ser,'b');

%ylabel('SER');
%xlabel('SNR');
data = [snr_min:step:snr_max; ser]
save STBC_OFDM_EST_LS_data.am data -ascii;
clear all;