% By Alger 2007
% 2Tx+1Rx SFBC+OFDM system in Rayleigh fading channel using Jakes' model
% Simplified Version

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%      Initialing       %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [BER_ZF, BER_DF, BER_SML, BER_JML, BER_ZF_new, BER_DF_new, BER_SML_new, BER_JML_new] = SFBCs(SNRindB, fm, d, M)
%clear;SNRindB = 0;fm = 0;d = 1;M = 12;

% OFDM parameters
fs = 800000;                                      % sampling rate/Bandwidth
Ncarr = 128;                                      % number of subcarriers
GIlen = 160;                                      % length of OFDM symbol plus cyclic prefix
CPlen = GIlen - Ncarr;                            % length of cyclic prefix

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%    Channel source     %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% generate infomation bits
X_info = randint(1,Ncarr);                        % binary bits
X_send = 2*X_info - 1;                            % BPSK modulation

% separate the sent bits to 2 branches
X_send1(1:2:Ncarr) =       X_send(1:2:Ncarr) ;    % transmitter1
X_send1(2:2:Ncarr) = -conj(X_send(2:2:Ncarr));    % [X1, -conj(X2)]

X_send2(1:2:Ncarr) =       X_send(2:2:Ncarr) ;    % transmitter2
X_send2(2:2:Ncarr) =  conj(X_send(1:2:Ncarr));    % [X2,  conj(X1)]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%         IFFT          %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ifft with length Slen
z_ifft1 = sqrt(Ncarr)*ifft(X_send1,Ncarr);        % sqrt(Ncarr) should be mutiplied to maintain the Power
z_ifft2 = sqrt(Ncarr)*ifft(X_send2,Ncarr);        % sqrt(Ncarr) should be mutiplied to maintain the Power

% adding cyclic prefix
Z_ofdm1 = [z_ifft1(Ncarr-CPlen+1:Ncarr) z_ifft1]; % final OFDM sampled signals with CP
Z_ofdm2 = [z_ifft2(Ncarr-CPlen+1:Ncarr) z_ifft2]; % final OFDM sampled signals with CP

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%        channel        %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% generate random channel matrix "h" with rayleigh distribution and Doppler shift by Jakes' model
h1 = zeros(M , Ncarr);                            % channel matrix
h2 = zeros(M , Ncarr);                            % channel matrix
r =  zeros(1 , Ncarr);                            % received signal

% multi paths channel following exponential delay model
TapGain = sqrt( (1-exp(-1/d)) / (1-exp(-M/d)) / 2 );
for i = 1:M
    h1(i,:) = TapGain * jakes(Ncarr,fm,fs);
    h2(i,:) = TapGain * jakes(Ncarr,fm,fs);
    r = r + h1(i,:).*Z_ofdm1( (CPlen-i+2):(GIlen-i+1) ) + h2(i,:).*Z_ofdm2( (CPlen-i+2):(GIlen-i+1) );
    TapGain = TapGain * exp(-1/d/2); 
end

R = awgn(r,SNRindB);                              % AWGN
%R = r;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%         IFFT          %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% fft
z_fft = fft(R,Ncarr)/sqrt(Ncarr);                 % 1/sqrt(Ncarr) should be multiplied

% separate the signals on odd and even subcarriers
X_odd  = z_fft(1:2:Ncarr);                        % R1 =  h11*     X1  + h21*     X2
X_even = z_fft(2:2:Ncarr);                        % R2 = -h12*conj(X2) - h22*conj(X1)

% Calculate the frequecy domain response
% | C A |
% | B C |
h_est1 = sum(h1,2)/Ncarr;
C1 = fft([h_est1;  zeros(Ncarr-M,1)]);
h_est1 = h1 * exp( sqrt(-1)*2*pi*(0:(Ncarr-1)).'/Ncarr) /Ncarr;
A1 = fft([h_est1;  zeros(Ncarr-M,1)]);
h_est1 = h1 * exp(-sqrt(-1)*2*pi*(0:(Ncarr-1)).'/Ncarr) /Ncarr;
B1 = fft([h_est1;  zeros(Ncarr-M,1)]);

h_est2 = sum(h2,2)/Ncarr;
C2 = fft([h_est2;  zeros(Ncarr-M,1)]);
h_est2 = h2 * exp( sqrt(-1)*2*pi*(0:(Ncarr-1)).'/Ncarr) /Ncarr;
A2 = fft([h_est2;  zeros(Ncarr-M,1)]);
h_est2 = h2 * exp(-sqrt(-1)*2*pi*(0:(Ncarr-1)).'/Ncarr) /Ncarr;
B2 = fft([h_est2;  zeros(Ncarr-M,1)]);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Conventional detection %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

H11_1 = C1(1:2:Ncarr); H21_1 = C2(1:2:Ncarr);      % R1 =  h11*     X1 + h21*     X2
H12_1 = C1(2:2:Ncarr); H22_1 = C2(2:2:Ncarr);      % R2 = -h12*conj(X2)+ h22*conj(X1)

X_receive_11 =  DF_decode(H11_1, H12_1, H21_1, H22_1, X_odd, X_even, Ncarr);
X_receive_12 =  ZF_decode(H11_1, H12_1, H21_1, H22_1, X_odd, X_even, Ncarr);
X_receive_13 = SML_decode(H11_1, H12_1, H21_1, H22_1, X_odd, X_even, Ncarr);
X_receive_14 = JML_decode(H11_1, H12_1, H21_1, H22_1, X_odd, X_even, Ncarr);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%     New detection     %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

H11_2 = C1(1:2:Ncarr) + A2(2:2:Ncarr); H21_2 = C2(1:2:Ncarr) - A1(2:2:Ncarr);
H12_2 = C1(2:2:Ncarr) - B2(1:2:Ncarr); H22_2 = C2(2:2:Ncarr) + B1(1:2:Ncarr);

X_receive_21 =  DF_decode(H11_2, H12_2, H21_2, H22_2, X_odd, X_even, Ncarr);
X_receive_22 =  ZF_decode(H11_2, H12_2, H21_2, H22_2, X_odd, X_even, Ncarr);
X_receive_23 = SML_decode(H11_2, H12_2, H21_2, H22_2, X_odd, X_even, Ncarr);
X_receive_24 = JML_decode(H11_2, H12_2, H21_2, H22_2, X_odd, X_even, Ncarr);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%      Demodulate       %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% demodulate the receive symbols
X_sink_12 = (1 + sign(real(X_receive_12)))/2;      % BPSK demodulation
X_sink_13 = (1 + sign(real(X_receive_13)))/2;
X_sink_22 = (1 + sign(real(X_receive_22)))/2;
X_sink_23 = (1 + sign(real(X_receive_23)))/2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%     Error counts      %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% bit error ratio
[bit,ratio] = biterr(X_info,X_receive_11);         % Coventional DF Detection
BER_DF = ratio;
[bit,ratio] = biterr(X_info,X_receive_21);         % New DF Detection
BER_DF_new = ratio;
[bit,ratio] = biterr(X_info,X_sink_12);            % Coventional ZF Detection
BER_ZF = ratio;
[bit,ratio] = biterr(X_info,X_sink_22);            % New ZF Detection
BER_ZF_new = ratio;
[bit,ratio] = biterr(X_info,X_sink_13);            % Coventional SML Detection
BER_SML = ratio;
[bit,ratio] = biterr(X_info,X_sink_23);            % New SML Detection
BER_SML_new = ratio;
[bit,ratio] = biterr(X_info,X_receive_14);         % Coventional JML Detection
BER_JML = ratio;
[bit,ratio] = biterr(X_info,X_receive_24);         % New JML Detection
BER_JML_new = ratio;