% By Alger 2007
% Basic BPSK+OFDM system 1Tx+1Rx in rayleigh fading channel using Jakes' model

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%      Initialing       %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function BER = OFDMs(SNRindB, fm, d, M)
%clear;SNRindB = 10;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(Ncarr,1);                        % binary bits
X_send = 2*X_info - 1;                            % BPSK modulation

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

% ifft with length Ncarr
z_ifft = (sqrt(Ncarr)*ifft(X_send,Ncarr)).';      % sqrt(Ncarr) should be mutiplied to maintain the Power

% adding cyclic prefix
Z_ofdm = [z_ifft(Ncarr-CPlen+1:Ncarr) z_ifft];    % final OFDM sampled signals with CP

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

% generate random channel matrix "h" with rayleigh distribution and Doppler shift by Jakes' model
h = 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)) );
for i = 1:M
    h(i,:) = TapGain * jakes(Ncarr,fm,fs);
    r = r + h(i,:) .* Z_ofdm( (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

% Calculate the frequecy domain response
h_est = sum(h,2)/Ncarr;
H = fft([h_est;  zeros(Ncarr-M,1)]);

% equlization in frequence domain
X_receive = z_fft.'./H;                             % the received symbols

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

% demodulate the receive symbols
X_sink = (1 + sign(real(X_receive)))/2;

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

% bit error ratio
[bit,ratio] = biterr(X_info,X_sink);
BER = ratio;