%    本程序中得到的收端OFDM信號的頻譜波形，是與其發端信號的排步有關的。在發端的
% 載波安排上，128個載波有前后各32個載波是null載波(如果這前后各32各載波是帶外頻段，
% 那么理論上它們都應該是零！)，中間的64個載波是數據載波。這樣的排步很明顯就是一個
% 兩邊低，中間高的頻譜形式。所以，收端也應該是這個輪廓。
clc;clear all;close all;echo off;tic;
% -------------------------------------------------------------------
%                      Parameter Definition
% --------------------------------------------------------------
Fd    = 1;             % symbol rate (1Hz)
Fs    = 1*Fd;          % number of sample per symbol
M     = 4;             % kind(range) of symbol (0,1,2,3)
Ndata = 1024;          % all transmitted data symbol 
Sdata = 64;            % 64 data symbol per frame to ifft
Slen  = 128;           % 128 length symbol for IFFT 
Nsym  = Ndata/Sdata;   % number of frames -> Nsym frame
GIlen = 144;           % symbol with GI insertion  GIlen = Slen + GI
GI    = 16;            % guard interval length
% ----------------------------------------------------------------
%                        Vector Initialization
% ----------------------------------------------------------------
X  = zeros(Ndata,1);
Y1 = zeros(Ndata,1);
Y2 = zeros(Ndata,1);
Y3 = zeros(Slen,1);
z0 = zeros(Slen,1);
z1 = zeros(Ndata/Sdata*Slen,1);
g  = zeros(GIlen,1);
z2 = zeros(GIlen*Nsym,1);
z3 = zeros(GIlen*Nsym,1);
% random integer generation by M kinds 
X = randint(Ndata, 1, M);
% digital symbol mapped as analog symbol
% Y1 is a Ndata-by-2 matrix, is changed into Y2 by "amodce"
Y1 = modmap(X, Fd, Fs, 'qask', M);
% covert to complex number
Y2 = amodce(Y1,1,'qam');
% figure(1);
% scatterplot(Y2,length(Y2),0,'bo');grid on;
scatterplot(Y2,Fd,0,'bo');grid on;
title('4-QAM Constellation');
Tx_spectrum = zeros(size(Y3));
for j=1:Nsym;
    for i=1:Sdata;
        Y3(i+Slen/2-Sdata/2,1)=Y2(i+(j-1)*Sdata,1);
        Tx_spectrum = Tx_spectrum + abs(Y3);
    end
    z0=ifft(Y3);
    
    for i=1:Slen;    % generate time-domain vector, z1, without GI
        z1(((j-1)*Slen)+i)=z0(i,1);
    end
    %
    for i=1:Slen;
        g(i+16)=z0(i,1);
    end
    
    for i=1:GI;
        g(i)=z0(i+Slen-GI,1);    
    end
    for i=1:GIlen;  % generate time-domain vector, z2, with GI
        z2(((j-1)*GIlen)+i)=g(i,1);
    end
    
end
% graph on time domain
figure(2);
f = linspace(-Sdata,Sdata,length(z1));
plot(f,abs(z1));
   
Y4 = fft(z1);  % FFT operation, at receiver
% if Y4 is under 0.01 Y4=0.001
for j=1:Ndata/Sdata*Slen;
if abs(Y4(j)) < 0.01
    Y4(j)=0.01;
end
end
Y4 = 10*log10(abs(Y4));  %  Y4 is used for spectrum display.
% graph on frequency domain
figure(3);
f = linspace(-Sdata,Sdata,length(Y4));
plot(f,Y4);grid on;
axis([-Slen/2 Slen/2 -20 20]);
title('Received OFDM signal spectrum');
figure(4);
f = linspace(-Sdata,Sdata,length(Y3));
plot(f,abs(Y3),'b-','LineWidth',6);grid on;
axis([-Slen/2 Slen/2 -2 2]);
title('Transmitted OFDM signal spectrum');
simulation_time = toc