%************************************************************************
%Calculating MSE in MIMO-OFDM with 2x2 antenna configuration
%Alamouti Coding 
%Num_RxAnt: the number of receive antenna
%M: M-ary digital modulation
%
%
%************************************************************************

clc;
clear all;
%repeat for different receive antennas
M = 4; %M-ary QAM Modulation Scheme
SNR_max = 10;
Bn = 200; %the number of OFDM symbol will being transmitted
randn('state',0); %Remove it if you want a random start of the randn generator
randn('state',0); %Remove it if you want a random start of the randn generator
z = 1;
%Creating Fl matrix for decoding Alamouti code
for n = 1:64
    for l = 1:16
        F(n,l) = exp(-(j*2*pi*(l-1)*(n-1))/64);
    end
end
l = 0;
n = 0;
Fl = [F zeros(size(F));zeros(size(F)) F];
Cn = 1;

for k = 1:0.5:SNR_max
    snr = 10.^(k/10);
    sig = 0.5/snr;
    Smv = floor(M*rand(2,64*Bn)); %transmitted alphabet
    Sm = exp(j*2*pi/M*Smv)/sqrt(2); %transmitted symbols
    %Sm = 2x(64*Bn) matrix, the block of input is going to be transmitted with constant SNR
    
    MSE = 0;
    
    %calculating y in time domain and change it to frequency domain then apply detection part
    for p = 1:Bn
        Sof = Sm(:,(1+(p-1)*64:p*64)); %transmitted OFDM block with constant SNR
        Sofc = conj(Sof);
        %Taking IFFT of the block
        Sot = [ifft(Sof(1,:));ifft(Sof(2,:))];
        Sotc = [ifft(Sofc(1,:));ifft(Sofc(2,:))];
        %Creating channel matrix with l = 16 multipaths
        for l = 0:15
            sig2 = 0.395*exp(-0.5*l);
            h(:,l+1) = (randn(4,1)+j*rand(4,1))/sqrt(2)*sqrt(sig2);
        end

        %******************************************************************
        %LS channel estimation
        %Creating Toeplitz matrix for each of impulse response separately and calculation corresponding received signal y
        ns = sqrt(sig).*(randn(80,2)+j*randn(80,2)); %Creating noise matrix in time domain
        
        y11 = toeplitz([h(1,:) zeros(1,64)],[h(1,1) zeros(1,79)])*[Sot(1,:) Sot(1,1:16)].';
        y12 = toeplitz([h(2,:) zeros(1,64)],[h(2,1) zeros(1,79)])*[Sot(2,:) Sot(2,1:16)].';
        y1 = y11 + y12 + ns(:,1); 
        
        y21 = toeplitz([h(3,:) zeros(1,64)],[h(3,1) zeros(1,79)])*[Sot(1,:) Sot(1,1:16)].';
        y22 = toeplitz([h(4,:) zeros(1,64)],[h(4,1) zeros(1,79)])*[Sot(2,:) Sot(2,1:16)].';
        y2 = y21 + y22 + ns(:,2); 
        
        %Removing CP and transform the signals to frequency domain
        ya1 = [y1(65:79,1);y1(16:64,1)];
        Y1 = fft(ya1);
        
        ya2 = [y2(65:79,1);y2(16:64,1)];
        Y2 = fft(ya2);
        
        %Calculating CIR, first calculating diagonal matrix for input in frequency domain
        X = [diag(Sof(1,:)) diag(Sof(2,:))];
        w = X*Fl;
        
        h_1 = inv(w'*w)*w'*Y1; %h_1 = [h_11;h_12] = []32x1
        h_2 = inv(w'*w)*w'*Y2; %h_2 = [h_21;h_22] = []32x1
        
        h_11 = h_1(1:16,1);
        h_12 = h_1(17:32,1);
        
        h_21 = h_2(1:16,1);
        h_22 = h_2(17:32,1);
        
        %LS channel estimation
        h_=[h_11.';h_12.';h_21.'; h_22.'];
        
        %End of LS channel estimation
        %******************************************************************
        
        %Calculating MSE
        E = h - h_;
        mse = (trace(E'*E))/64;
        MSE =MSE + mse;
    end
    MSEF(Cn) = MSE/k;
    Cn = Cn + 1;    
end
Snr = 1:0.5:SNR_max;

%plot the bit error rate
semilogy(Snr,MSEF,'b^-','LineWidth',2.0);
xlabel('SNR [dB]');
ylabel('MSE');
title('Mean Square Error Evaluation for OFDM system with Alamouti coding, 4-QAM Modulation');
grid on;