% Lab#6: OFDM Principle
% by D.C.Chang NCU/CE  2021/4/6
%
clear all

channel_model=0;  % 0 for ideal channel
snr=30;           % dB
cp=16;

symbol_num=2000;
ofdm_length=64+cp;
tx=zeros(1,ofdm_length*symbol_num);
qam=zeros(1,symbol_num*52);
qamd=zeros(1,symbol_num*52);
upper_fft_id=[1:26]+1;
lower_fft_id=[38:63]+1;  

%% 64-QAM Generation
      qam_f=7;
      QAM64set=[-qam_f:2:qam_f];
      for i=1:symbol_num*52
	        qam_i=QAM64set(round(qam_f*rand)+1);
	        qam_q=QAM64set(round(qam_f*rand)+1);
            qam(i)=qam_i+j*qam_q;
      end
%%

%% IFFT Modulation    
    for i=1:symbol_num
       ifft_in=zeros(1,64);      
       ifft_in(upper_fft_id)=qam([52*(i-1)+1:52*(i-1)+26]);
       ifft_in(lower_fft_id)=qam([52*(i-1)+27:52*i]);
       tx_raw=ifft(ifft_in);
       tx_cp=[tx_raw(64-cp+1:64) tx_raw];
       tx(ofdm_length*(i-1)+1:ofdm_length*i)=tx_cp;
    end    
%%
    
    %% Baseband Channel
    h0=1;                                                   % ideal channel
    %--------------channel model--------------------    
    h1=                        % channel length=15
    h2=                        % channel length=19
    %------------------------------------------------
    
    if channel_model==0
       ch=h0;
    elseif channel_model==1 
       ch=h1;  
    else     
       ch=h2;
    end   
    rxt=conv(tx,ch);
    Heq=fft(ch,64); 
   %%
   
   %% Adding Noise
    Ps=tx*tx'/length(tx);
    tx=tx/sqrt(Ps);      % tx signal of normalized power 
    Pn=1/10^(snr/10);  % noise power    
    noise=sqrt(Pn)*crandn(1,length(rxt))/sqrt(2);
    rx=rxt+noise;
  %%   
  
 %% Receiver FFT and FEQ
    for i=1:symbol_num
        %--------------DFT output: FFT +FEQ
      
        %-----------------------------------
        qamd(52*(i-1)+1:52*i)=[fft_out(upper_fft_id) fft_out(lower_fft_id)];
    end        
    
 %%   
  
    figure(1)
    n=2048;
    [X,F]=pwelch(tx,[],[],n,20*10^6);
    plot((F-10*10^6)/10^6,10*log10([X(n/2+1:n)' X(1:n/2)']),'b') 
    xlabel('frequency (MHz)')
    ylabel('Power Spectrum Density (dB)')
   
    figure(2)
    plot(qamd,'.r')
    xlabel('I-channel')
    ylabel('Q-channel')
    axis([-10 10 -10 10])