%************************************************************************
%Alamouti Coding with 2xNr antenna configuration
%
%Num_RxAnt: the number of receive antenna
%M: M-ary digital modulation
%
%
%************************************************************************

clc;
clear all;
Num_RxAnt = 6;
M = 4; %QPSK Modulation Scheme
SNR_max = 16;
N = 2000; %the pair number of 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;
for k = 1:SNR_max
    A = floor(M*rand(2,N)); %transmitted alphabet
    Str = exp(j*2*pi/M*A)/sqrt(2); %transmitted symbols
    %Str=2xN matrix, the block of input is going to be transmitted with constant SNR
    for i = 1:N
        S = [Str(1,i);Str(2,i)]; %making STC matrix
        snr = 10.^(k/10);
        sig = 0.5/snr;
        Ns = sqrt(sig).*(rand(2.*Num_RxAnt,1) + j*(rand(2.*Num_RxAnt,1)));
        H = [];
        Hsp_code = [];
        for w = 1:Num_RxAnt
            h = (rand(1,2) + j*(rand(1,2)))/sqrt(2);
            H = [H;h(1) h(2)]; %channel matrix
            Hsp_code = [Hsp_code;h(1) h(2);h(2)' -h(1)'];
        end

        %******************************************************************
        %Channel estimation
        Z = hadamard(2); %Orthogonal Training data
        noise = sqrt(sig).*(randn(Num_RxAnt,2) + j*(randn(Num_RxAnt,2)));
        S_train = H*Z + noise;
        H_est = S_train*Z'*(inv(Z*Z'));
        w = 0;
        %H_est: estimated channel matrix using LS algorithm and optimal orthogonal training data
        Hls = [];
        for w = 1:Num_RxAnt
            Hls = [Hls;H_est(w,1) H_est(w,2);H_est(w,2)' -H_est(w,1)']; %channel matrix
        end

        %end of channel estimation
        %******************************************************************
        
        %Perfect CSI signals
        r = Hsp_code*S + Ns; %receive signal on the receiver
        S_p = Hsp_code'*r;
        angp = angle(S_p); %receive angles
        B = mod(round(angp/(2*pi/M)),M); %receive alphabet, its another
        %overal output for the all of the transmission symbol
        Sesp(:,i) = [B(1);B(2)];
        %Str = transmitted symbols
        %Sesp = Perfect CSI received symbols
        
        %LS signals
        r = Hsp_code*S + Ns; %receive signal on the receiver
        S_ls = Hls'*r;
        ang_ls = angle(S_ls); %receive angles
        Bls = mod(round(ang_ls/(2*pi/M)),M); %receive alphabet, its another
        %overal output for the all of the transmission symbol
        Ses_ls(:,i) = [Bls(1);Bls(2)];
        %Str = transmitted symbols
        %Sesls = LS received symbols
        end
        
        
        %calculating bit error rate
        BER = 0;
        ERRORp = 0;
        ERRORls = 0;
        for p = 1:2
            for i = 1:N
                if A(p,i) ~= Sesp(p,i)
                    ERRORp = ERRORp + 1;
                end
                if A(p,i) ~= Ses_ls(p,i)
                    ERRORls = ERRORls + 1;
                end
            end
        end
        berp(z) = ERRORp/(2*N);
        berls(z) = ERRORls/(2*N);
        z = z + 1;
end
Snr = 1:SNR_max;
berp1 = berp;
berls1 = berls;

%plot the bit error rate
semilogy(Snr,berp1,'b-',Snr,berls1,'g+-','LineWidth',2.0);
xlabel('SNR [dB]');
ylabel('BER');
title('Bit Error Rate Evaluation for MIMO Alamouti coding with 2x6 antenna configuration, QPSK Modulation');
legend('Perfect CSI','LS channel estimation');
grid on;