%************************************************************************
%Alamouti Coding for OFDM with 2x2 antenna configuration
%
%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];

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
    
    %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
        %******************************************************************
        
        %Creating Toeplitz matrix for each of impulse response separately and calculation corresponding received signal y
        ns = sqrt(sig).*(randn(80,4)+j*randn(80,4)); %Creating noise matrix in time domain
        
        y011 = toeplitz([h(1,:) zeros(1,64)],[h(1,1) zeros(1,79)])*[Sot(1,:) Sot(1,1:16)].';
        y012 = toeplitz([h(2,:) zeros(1,64)],[h(2,1) zeros(1,79)])*[Sot(2,:) Sot(2,1:16)].';
        y0 = y011 + y012 + ns(:,1); %at time t1 for received antenna 1
        
        y111 = toeplitz([h(1,:) zeros(1,64)],[h(1,1) zeros(1,79)])*-[Sotc(2,:) Sotc(2,1:16)].';
        y112 = toeplitz([h(2,:) zeros(1,64)],[h(2,1) zeros(1,79)])*[Sotc(1,:) Sotc(1,1:16)].';
        y1 = y111 + y112 + ns(:,2); %at time t2 for received antenna 1
        
        y021 = toeplitz([h(3,:) zeros(1,64)],[h(3,1) zeros(1,79)])*[Sot(1,:) Sot(1,1:16)].';
        y022 = toeplitz([h(4,:) zeros(1,64)],[h(4,1) zeros(1,79)])*[Sot(2,:) Sot(2,1:16)].';
        y2 = y021 + y022 + ns(:,3); %at time t1 for received antenna 2
        
        y121 = toeplitz([h(3,:) zeros(1,64)],[h(3,1) zeros(1,79)])*-[Sotc(2,:) Sotc(2,1:16)].';
        y122 = toeplitz([h(4,:) zeros(1,64)],[h(4,1) zeros(1,79)])*[Sotc(1,:) Sotc(1,1:16)].';
        y3 = y121 + y122 + ns(:,4); %at time t2 for received antenna 2
        
        %Removing CP and transform the signals to frequency domain
        ya0 = [y0(65:79,1);y0(16:64,1)];
        YA0 = fft(ya0);
        
        ya1 = [y1(65:79,1);y1(16:64,1)];
        YA1 = fft(ya1);
        
        ya2 = [y2(65:79,1);y2(16:64,1)];
        YA2 = fft(ya2);
        
        ya3 = [y3(65:79,1);y3(16:64,1)];
        YA3 = fft(ya3);
        
        Y = [YA0 conj(YA1) YA2 conj(YA3)];
        
        %Creating subcarrier matrix in frequency domain
        ht = [(h(1,:)).' (h(2,:)).';(h(3,:)).' (h(4,:)).']; %Perfect CSI
        H = Fl*ht;
        
        
        %h_=[h_11.';h_12.';h_21.';h_22.']  LS channel estimation
        ht_ = [h_11 h_12;h_21 h_22];
        H_ = Fl*ht_;
        %Alamouti Decoding
        for Kc = 1:64
            YC = Y(Kc,:);
            Yc = YC.';
            
            %signals with perfect CSI
            Hc = [H(Kc,1) H(Kc,2);H(Kc,2)' -H(Kc,1)';H(Kc+64,1) H(Kc+64,2);H(Kc+64,2)' -H(Kc+64,1)'];            
            X_perf = (Hc)'*Yc;             
            %Selecting 2 symbols from OFDM block, send it and retrive it until the end of the OFDM block as well as for Bn OFDM block
            angp = angle(X_perf);
            S_perf = mod(round(angp/(2*pi/M)),M); %received alphabet, its another
            Sp_det_mat(:,(64*(p-1)+Kc)) = S_perf; %Output one OFDM transmission block
            
            %signals with LS channel estimation
            Hc_ = [H_(Kc,1) H_(Kc,2);H_(Kc,2)' -H_(Kc,1)';H_(Kc+64,1) H_(Kc+64,2);H_(Kc+64,2)' -H_(Kc+64,1)'];
            X_est = (Hc_)'*Yc; 
            ang_est = angle(X_est); 
            S_est = mod(round(ang_est/(2*pi/M)),M); 
            S_det_mat_est(:,(64*(p-1)+Kc)) = S_est; 
        end
    end
        
        %Calculating bit error rate
        ERRORp = 0;
        ERROR_est = 0;
        for p = 1:2
            for u = 1:64*Bn
                if Smv(p,u) ~= Sp_det_mat(p,u)
                    ERRORp = ERRORp + 1;
                end
                if Smv(p,u) ~= S_det_mat_est(p,u)
                    ERROR_est = ERROR_est + 1;
                end
            end
        end
        berp(z) = ERRORp/(2*64*Bn);
        ber_est(z) = ERROR_est/(2*64*Bn);
        z = z + 1;
end
Snr = 1:0.5:SNR_max;
ber2p = berp;
ber2_est = ber_est;

%plot the bit error rate
semilogy(Snr,ber2p,'g*-',Snr,ber2_est,'b^-','LineWidth',2.0);
xlabel('SNR [dB]');
ylabel('BER');
title('Bit Error Rate Evaluation for OFDM system with Alamouti coding, 4-QAM Modulation');
legend('Perfect CSI','LS channel estimation');
grid on;