% 包含有未經信道估計(figure2),經過信道估計(導頻)(figure3),經過信道估計(滿幀)(figure3)

clear all;
clc;

% ****************************initialize the parametres*********************************

Bit_Num=4;  %每個符號比特數
IFFT_Num=256; %采樣點數
Carrier_Num=150; %子載波數
OFDM_per_Symbol=30; %每幀OFDM信號含有的OFDM信號的個數
SNR_N=20; %信噪比間隔
SNR_Bit=50; %每比特信噪比
GI=20; %保護間隔長度
LI=5 ; %導頻間隔
Np=ceil(Carrier_Num/LI)+1; %導頻數    加1保證最后一行子載波中也有導頻
Carriers=1:Carrier_Num+Np;
N_Num=Bit_Num*Carrier_Num*OFDM_per_Symbol;

% ***********************Generate the random binary stream for transmit test******
 
BitTx=randint(1,N_Num);

% ********************Modulate (generates 16QAM symbols)*****************

SymQAM=tx_16QAM_mod(BitTx);

figure(1)
plot(SymQAM,'b*');
axis([-5 5 -5 5]);

%*******************Generate the pilot signal************
  
train_sym1=randint(1,2*OFDM_per_Symbol);
train_sym2=randint(1,2*OFDM_per_Symbol);
t=2*(train_sym1.*2-1)+(train_sym2.*2-1);
treal=t(1:2:2*OFDM_per_Symbol);
timage=t(2:2:2*OFDM_per_Symbol);
training_symbols1=treal+i.*timage;
training_symbols2=training_symbols1.';
training_symbols=repmat(training_symbols2,1,Np);

pilot=1:LI+1:Carrier_Num+Np;
if length(pilot)~=Np
    pilot=[pilot,Carrier_Num+Np];
end

% **********************P/S transform*********************

SymQAM1=reshape(SymQAM,Carrier_Num,OFDM_per_Symbol).';

% *******************Insert the pilot symbols****************

signal=1:Carrier_Num+Np;
signal(pilot)=[];
SymQAM2(:,pilot)=training_symbols;
SymQAM2(:,signal)=SymQAM1;
SymQAM3=SymQAM2.';

% **********************IFFT transform******************************

% Input:  SymQAM3(Carrier_Num+Np,OFDM_per_Symbol)    Output:  SymTx(1,(IFFT_Num+GI)*OFDM_per_Symbol)
SymIFFT=ifft(SymQAM3,IFFT_Num,1);
SymIFFT_test=ifft(SymQAM1.',IFFT_Num,1);

%加保護間隔(循環前綴)
SymTxtmp=zeros(IFFT_Num+GI,OFDM_per_Symbol);
SymTxtmp(1:GI,:)=SymIFFT(IFFT_Num-GI+1:IFFT_Num,:);
SymTxtmp(GI+1:GI+IFFT_Num,:)=SymIFFT;
SymTx=reshape(SymTxtmp,1,(IFFT_Num+GI)*OFDM_per_Symbol);

% *************************Channel*********************************


fd=100; %多普勒頻移
r=6;   %多徑數
d=[2 3 4 5 9 13]; % 信道延遲
a=[0.123 0.3 0.4 0.5 0.7 0.8]; %多徑的幅度
T=1;  %系統采樣周期
th=[90 0 72 144 216 288]*pi./180;%相移
h=zeros(1,Carrier_Num);
hh=[];
for k=1:r
    h1=a(k)*exp(j*((2*pi*T*fd*d(k)/Carrier_Num)));
    hh=[hh,h1];
end
h(d+1)=hh;



channel1=zeros(size(SymTx));
channel1(1+d(1):length(SymTx))=hh(1)*SymTx(1:length(SymTx)-d(1));
channel2=zeros(size(SymTx));
channel2(1+d(2):length(SymTx))=hh(2)*SymTx(1:length(SymTx)-d(2));
channel3=zeros(size(SymTx));
channel3(1+d(3):length(SymTx))=hh(3)*SymTx(1:length(SymTx)-d(3));
channel4=zeros(size(SymTx));
channel4(1+d(4):length(SymTx))=hh(4)*SymTx(1:length(SymTx)-d(4));
channel5=zeros(size(SymTx));
channel5(1+d(5):length(SymTx))=hh(5)*SymTx(1:length(SymTx)-d(5));
channel6=zeros(size(SymTx));
channel6(1+d(6):length(SymTx))=hh(6)*SymTx(1:length(SymTx)-d(6));

BER_plot=[];
%for snr_db=SNR_Bit:-SNR_N:-SNR_Bit
%for snr_db=SNR_Bit:(-SNR_N):0 % ???????????????????????????????????????????倒序才可產生對應圖  
snr_db=20;
code_power=(norm(SymTx))^2/(length(SymTx));  %符號平均功率
bit_power=code_power/Bit_Num; %比特功率
noise_power=10*log10(bit_power/(10^(snr_db/10))); %噪聲功率
noise=wgn(1,length(SymTx),noise_power,'complex'); %產生高斯白噪聲信號
SymTx=SymTx+noise;

SymTxChannel=SymTx+channel1+channel2+channel3+channel4+channel5+channel6;

% ************************FFT  transform*********************************

%去保護間隔
SymRxtmp=reshape(SymTxChannel,IFFT_Num+GI,OFDM_per_Symbol);
SymRx=SymRxtmp(GI+1:GI+IFFT_Num,:);

SymFFT=fft(SymRx,IFFT_Num,1);
SymFFT=SymFFT(1:Carrier_Num+Np,:);

figure(2)
plot(SymFFT,'x');
axis([-5 5 -5 5]);

%************************Channel Estimation*****************************

SymFFT1=SymFFT.';
SymFFT2=SymFFT1(:,signal);

HLs=ChannelEstimate(SymFFT,OFDM_per_Symbol,Np,training_symbols,Carrier_Num,LI,pilot);
SymFFT3=SymFFT2./HLs;
SymFFT=SymFFT3.';

HLs_test2=ChannelEstimate_test2(SymQAM1,SymFFT2,OFDM_per_Symbol);
SymFFT3_test2=SymFFT2./HLs_test2;
SymFFT_test2=SymFFT3_test2.';

figure(3)
plot(SymFFT,'x');
axis([-5 5 -5 5]);

figure(4)
plot(SymFFT_test2,'x');
axis([-5 5 -5 5]);


% ***********************Demodulate ***************************************

SymDeQAM=rx_16QAM_mod(SymFFT,OFDM_per_Symbol,Carrier_Num,Bit_Num);
SymDeQAM_test2=rx_16QAM_mod(SymFFT_test2,OFDM_per_Symbol,Carrier_Num,Bit_Num);

%***********************Error detection************************************

Error_Num=length(find(BitTx~=SymDeQAM));
BER=Error_Num./(4*N_Num);
fprintf('\n SNR is %1.2f dB ,BER is %e\n',snr_db,BER);

Error_Num_test2=length(find(BitTx~=SymDeQAM_test2));
BER_test2=Error_Num_test2./(4*N_Num);
fprintf('\n SNR is %1.2f dB ,BER_test2 is %e\n',snr_db,BER_test2);

BER_plot=[BER_plot,BER];
%end