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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%      Initialing       %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function BER = OFDM(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 = ones(M , GIlen+M);                            % channel matrix
r = zeros(1 , GIlen+M);                           % 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 , i:(i-1+GIlen) ) = jakes(GIlen,fm,fs);
    h(i,:) = TapGain * h(i,:);
    r = r + h(i,:).*[ones(1,i-1) Z_ofdm ones(1,M-i+1)];
    TapGain = TapGain * exp(-1/d/2);
end

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

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

% reduce the cyclic prefix
y_fft = R(1+CPlen:GIlen).';

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

% Calculate the frequecy domain response
h_est = sum(h(:,CPlen+1:CPlen+Ncarr),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;