MIMO_OFDM.m

%------------------------------------------
% EE359 final project Fall 2002
% Channel estimation for a MIMO-OFDM system
% By Shahriyar Matloub               
%------------------------------------------

clear all;
%close all;
i=sqrt（-1）;
Rayleigh=1;
AWGN=0;                             % for AWGN channel 
MMSE=0;                             % estimation technique
Nsc=64;                             % Number of subcarriers
Ng=16;                              % Cyclic prefix length
SNR_dB=[0 5 10 15 20 25 30 35 40];  % Signal to noise ratio
Mt=2;                               % Number of Tx antennas
Mr=2;                               % Number of Rx antennas
pilots=[1:Nsc/Ng:Nsc];              % pilot subcarriers 
DS=5;                              % Delay spread of channel
iteration_max=200;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Channel impulse response %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if （Rayleigh）
    N=50;
    fm=100;
    B=20e3;
    fd=（rand（1N）-0.5）*2*fm;
    theta=randn（1N）*2*pi;
    c=randn（1N）;
    c=c/sum（c.^2）;
    t=0:fm/B:10000*fm/B;
    Tc=zeros（size（t））;
    
    Ts=zeros（size（t））;
    for k=1:N
       Tc=c（k）*cos（2*pi*fd（k）*t+theta（k））+Tc;
       Ts=c（k）*sin（2*pi*fd（k）*t+theta（k））+Ts;
    end
    r=ones（Mt*Mr1）*（Tc.^2+Ts.^2）.^0.5;
    index=floor（rand（Mt*MrDS）*5000+1）;
end

MEE1=zeros（1length（SNR_dB））;
MEE2=zeros（1length（SNR_dB））;

for snrl=1:length（SNR_dB）
    snrl
    estimation_error1=zeros（Mt*MrNsc）;
    estimation_error2=zeros（Mt*MrNsc）;
    R1=besselj（02*pi*fm*（Nsc+Ng）/B）;
    sigma2=10^（-SNR_dB（snrl）/10）;
    aa=（1-R1^2）/（1-R1^2+sigma2）;
    bb=sigma2*R1/（1-R1^2+sigma2）;

    for iteration=1:iteration_max
        %iteration    
        if AWGN==1
            h=ones（Mt*Mr1）;
        else
            phi=rand*2*pi;
            h=r（index+iteration）*exp（j*phi）;
            %h=rand（Mt*MrDS）;
            h=h.*（ones（Mt*Mr1）*（exp（-0.5）.^[1:DS]））;
            h=h./（sqrt（sum（abs（h）.^22））*ones（1DS））;
        end



        CL=size（h2）;                                               % channel length
        data_time=zeros（MtNsc+Ng）;
        data_qam=zeros（MtNsc）;
        data_out=zeros（MrNsc）;
        output=zeros（MrNsc）;

        for tx=1:Mt
            data_b=0*round（rand（4Nsc））;                                  % data
            data_qam（tx:）=j*（2*（mod（data_b（1:）+data_b（2:）2）+2*data_b（1:））-3）+...
            2*（mod（data_b（3:）+data_b（4:）2）+2*data_b（3:））-3;
            for loop=1:Mt 
                data_qam（txpilots+loop-1）=（1+j）*（loop==tx）;              % pilots
            end
            data_time_temp=ifft（data_qam（tx:））;
            data_time（tx:）=[data_time_temp（end-Ng+1:end） data_time_temp];
        end
    
        for rx=1:Mr
            for tx=1:Mt
                output_temp=conv（data_time（tx:）h（（rx-1）*Mt+tx:））;
                output（rx:）=output_temp（Ng+1:Ng+Nsc）+output（rx:）;
            end
            np=（sum（abs（output（rx:））.^2）/length（output（rx:）））*sigma2;
            noise=（randn（size（