%%%%%%%%%%%  Developed for 16 QAM %%%%%%%%%%%%%%
clear all
close all
t0=clock;
%%%%%%%  INPUT PARAMETERS %%%%%%%%%%%
%%%%%  N=, Number of FFT Points
N=1024;
%%%%% M=, Number of OFDM Sub Carriers
M=512;
%%%%% L=, Number of Symbols
L=10;
MMMa=16;%%%% Modulation Method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% QPSK 16QAM
MMM1a=log2(MMMa);
%%%%% Tg=, Gurd Interval Ratio %%%%%%%%
Tg=0.1;
W=10;%%% Bandwidth(MHz)
DFF=W/M;
%%%%%% ts=, Effective Sumbol Duration (Micro Sec) %%%%%%%
ts=M/W;
%%%%%% D=, Gurd Interval Duration (Micro Sec) %%%%%%%
D=ts*Tg;
%%%%%% Ts=, Symbol Duration (Micro Sec) %%%%%%%%%%
Ts=ts+D;
st1=ts/N;
%%%%%%% Ng=, Number of Samples for Gurd Interval %%%%%%%
Ng=round(D/(ts/N));
%%%%%%% Nt=, Number of Samples for a Symbol %%%%%%%
Nt=(N+Ng);
%%%%%%% NNt=, Number of Samples for a Flame %%%%%%
NNt=Nt*L;
%%%%%%%%% Sampling frequency (usec) %%%%%%%%%%
tbw=1/st1;
%%%%%%% tsw=, Total Allocated Signal Bandwidth (MHz) %%%%%%%
tsw=(1/ts*M);%%%% W bandwidth 
ddf=tbw/NNt; %%%% Frequency resolution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dt=0.1;%%%%%%  delay spread
ds=st1;%%%%%%  delay spacing
NP=10;  %%%%%% number of delay path 

%%% Power Delay profile with  Exponential decay constant 
amp(1)=0.7854;
amp(2)=0.4861;
amp(3)=0.3009;
amp(4)=0.1863;
amp(5)=0.1153;
amp(6)=0.0714;
amp(7)=0.0442;
amp(8)=0.0273;
amp(9)=0.0169;
amp(10)=0.0105;
dts1(1:NP)=0:9;
%[amp,dst,dts1]=f_MUL(NP,dt,ds,st1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cn=30;%%% dB
cn1=10^(-cn/10);
 
LLL=1;%%% Number of trial
for KLL=1:LLL;
 if mod(KLL,1)==0,KLL
  end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  zxn1=zeros(1,NP); 
  zyn1=zeros(1,NP);
  zn1=zeros(1,NP);
  zn2=zeros(1,NP);
  phn=zeros(1,NP);
  amp1=zeros(1,NP);
  qd3=zeros(1,NNt);
  %%% Generation of Rayleigh fading
  %%% Amplitude
  zxn1(1:NP)=randn(1,NP);
  zyn1(1:NP)=randn(1,NP);
  zn1(1:NP)=zxn1(1:NP)+j*zyn1(1:NP);
  zn2(1:NP)=abs(zn1(1:NP))/sqrt(2);
  %%% Phase
  phn(1:NP)=atan2(zyn1(1:NP),zxn1(1:NP));
  amp1(1:NP)=zn2(1:NP).*amp(1:NP).*exp(j*phn(1:NP)); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ya=zeros(L,M);
x3=zeros(L,N);
L1=N/2+1-M/2;
L2=N/2+M/2;
%%%%%%%%%%%%%Data generation %%%%%%%%%%%%%%
ya(1:L,1:M)=randint(L,M,MMMa);
%%%%%%%%%%%%% Modulation %%%%%%%%%
x3(1:L,L1:L2)=f_EnMQAM(ya(1:L,1:M),MMMa);
%%%%%%%%%%%Converting from Frequency to Time domain signal %%%%%%%%%
 for k=1:L; 
 zz1(k,1:N)=ifft(fftshift(x3(k,1:N)))*sqrt(N);
 %%%% Adding Guard Interval 
 z3(k,1:Nt)=[zz1(k,N-Ng+1:N) zz1(k,1:N)];
end
  q13=z3';
  q73=q13(:)';
  %%%%%%%%%% Evaluation of Power Spectrum under AWGN condition 
  q73a=fftshift(fft(q73));
  pow1=q73a.*conj(q73a);
  pow2=10*log10(pow1/max(pow1));
  %%%%%%%%%%%
  po3=mean(q73.*conj(q73));
  c3=po3/2*cn1;
  %%%%%% Noise Generation %%%%%
     zxc(1:L,1:Nt)=sqrt(c3).*randn(L,Nt);
     zyc(1:L,1:Nt)=sqrt(c3).*randn(L,Nt);
     znc(1:L,1:Nt)=zxc(1:L,1:Nt)+j*zyc(1:L,1:Nt);
     znn1c=znc';
     znnc=znn1c(:)';
     
     qd3(1:NNt)=amp1(1).*q73(1:NNt);
     %%%% Generation of Multi-path fading signal
  for k=2:NP;
     qu23(k,1:NNt)=[q73(NNt+1-dts1(k):NNt) q73(1:NNt-dts1(k))];
     qu23(k,1:NNt)=amp1(k).*qu23(k,1:NNt);
     qd3(1:NNt)=qd3(1:NNt)+qu23(k,1:NNt);
  end
  %%%%%%%%%% Evaluation of Power Spectrum under multi-path fading condition 
  qd3a=fftshift(fft(qd3));
  pow3=qd3a.*conj(qd3a);
  pow4=10*log10(pow3/max(pow3));
  %%%% Adding AWGN
  q73x=qd3+znnc;
  q73y=q73+znnc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
imax=1;
%%% Removing Guard interval
for k=1:L;
     p93(k,1:Nt)=q73x(imax-1+(k-1)*Nt+(1:Nt));
     p993(k,1:N)=[p93(k,Ng+1:N+Ng)];
     %%%%%%%%%%  FFT %%%%%
     p813(k,1:N)=fftshift(fft(p993(k,1:N))/sqrt(N));
     HH3(1:M)=p813(1,L1:L2)./x3(1,L1:L2);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     p8x3(k,1:M)=p813(k,L1:L2)./HH3(1:M);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     p93y(k,1:Nt)=q73y(imax-1+(k-1)*Nt+(1:Nt));
     p993y(k,1:N)=[p93y(k,Ng+1:N+Ng)];
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%
     p813y(k,1:N)=fftshift(fft(p993y(k,1:N))/sqrt(N));
     p8y3(k,1:M)=p813y(k,L1:L2);
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
    
    %%% Demodulation of signal under AWGN
  pb1(1:L,1:M)=f_DeMQAM(p8x3(1:L,1:M),MMMa);%%%
  %%% evaluation of bit error rate performance
  per1(KLL)=biterr(pb1(2:L,1:M),ya(2:L,1:M))/(M*(L-1)*MMM1a);
  %%% Demodulation of signal under multi-path fading 
  pb2(1:L,1:M)=f_DeMQAM(p8y3(1:L,1:M),MMMa);%%%
  per2(KLL)=biterr(pb2(2:L,1:M),ya(2:L,1:M))/(M*(L-1)*MMM1a);
end
%%%% Final BER
 ERR1=mean(per1(1:LLL));
 ERR2=mean(per2(1:LLL));
 
 %%%%%%%%%%%%%%%%%%%%%%%%
 figure(1)%%%%%% Information Real Data in Freq
 stem(real(x3(2,1:N)))
 grid on
 figure(2)%%%%%% Information Imag Data in Freq
 stem(imag(x3(2,1:N)))
 grid on
 figure(3)%%%%%% Time domain signal
 plot(abs(z3(1,1:Nt)))
 grid on
 figure(4)
 plot(abs(z3(1,1:Ng)))
 grid on
 figure(5)
 plot(abs(z3(1,Nt-Ng:Nt)))
 grid on
 figure(6)
 plot(pow2)
 grid on
 figure(7)
 plot(pow4)
 grid on
 figure(8)
 plot(abs(HH3))
 grid on
 figure(9)
 stem(real(p8x3(2,1:M)))
 grid on
 figure(10)
 stem(imag(p8x3(2,1:M)))
 grid on
 figure(11)
 plot(p8x3(2:L,1:M),'x')
 grid on
 figure(12)
 plot(p8y3(2:L,1:M),'x')
 grid on
 
 

etime(clock,t0)/60