%ofdm1.m
%QPSK調制、AWGN信道下仿真程序

%**********************初始參數***************************
para=52;%并行信道數
fftlen=64;%FFT長度
noc=52;%載波數
nd=6;%每循環中OFDM符號數
ml=2;%調制水平:QPSK
gilen=16;%保護間隔長度(points)
ebn0=2;%信噪比
sr=250000;%OFDM symbol rate(250 ksyombol/s)
br=sr.*ml;%Bit rate per carrier

%**************************主循環部分**************************
nloop=2;%仿真的循環數
noe1=0;%信道解碼前錯誤數據數
nod1=0;%信道解碼前傳輸數據數
noe2=0;%信道解碼后錯誤數據數
nod2=0;%信道解碼后傳輸數據數

%**************************發射機*********************************
for iii=1:nloop
    
%**************************數據產生****************************
%信源編碼
t=[0:pi/25:2*pi];
xx=sin(t);
init=[-1:.1:1];
partition=[-1:.1:.9];
predictor=[0 1];
encode=dpcmenco(xx,init,partition,predictor);
encode2=reshape(encode,51,1);
bin=deci22bin(encode2,6);
recode=reshape(bin,1,306);

%信道編碼
k0=1;
G=[1 0 1 1 0 1 1;1 1 1 1 0 0 1];
channelencode=cnv_encd(G,k0,recode);

%******************串并轉換***********************
paradata=reshape(channelencode,para,nd*ml);%reshape:內建功能

%**************************QPSK調制*****************************
[ich,qch]=qpskmod1(paradata,para,nd,ml);
kmod=1/sqrt(2);%sqrt:內建功能
ich=ich.*kmod;
qch=qch.*kmod;

%數據映射
[ich1,qch1]=crmapping(ich,qch,fftlen,nd);  %pading zeros

%*******************IFFT************************
x=ich1+qch1.*i;
y=ifft(x);%ifft:內建功能
ich2=real(y);%real:內建功能
qch2=imag(y);%imag:內建功能

%*********插入保護間隔**********
[ich3,qch3]=giins1(ich2,qch2,fftlen,gilen,nd);

fftlen2=fftlen+gilen;

%----------Attenuation Calculation----------
spow=sum(ich3.^2+qch3.^2)/nd./para;%
attn=0.5*spow*sr/br*10.^(-ebn0/10);
attn=sqrt(attn);

%---------------AWGN addition---------------
[ich4,qch4]=comb(ich3,qch3,attn);

%******************去除保護間隔*********
[ich5,qch5]=girem1(ich4,qch4,fftlen2,gilen,nd);

%******************FFT******************
rx=ich5+qch5.*i;
ry=fft(rx);
ich6=real(ry);
qch6=imag(ry);

%載波逆映射
[ich7,qch7]=crdemapping(ich6,qch6,fftlen,nd);

%*****************解調*******************
ich7=ich7./kmod;
qch7=qch7./kmod;
[demodata]=qpskdemod1(ich7,qch7,para,nd,ml);

%**************并串轉換*****************
demodata1=reshape(demodata,1,para*nd*ml);

%信道解碼
[channeldecode]=viterbi(G,k0,demodata1);

%信源解碼
reshapechanneldecode=reshape(channeldecode,51,6);
deci=bin22deci(reshapechanneldecode);
redeci=reshape(deci,1,51);
codebook=[-1:.1:5.3];
decode=dpcmdeco(redeci,codebook,predictor);

%**************************Bit Error Rate(BER)*************************
noe10=sum(abs(demodata1-channelencode));
nod10=length(channelencode);
noe20=sum(abs(channeldecode-recode));
nod20=length(recode);
noe1=noe10+noe1;
nod1=nod10+nod1;
noe2=noe20+noe2;
nod2=nod20+nod2;
fprintf('%d\t%e\t%e\n',iii,noe10/nod10,noe20/nod20);

end
ber1=noe1/nod1;
ber2=noe2/nod2;
%**********************輸出結果***************************
fprintf('%f\t%e\t%e\t%d\t\n',ebn0,ber1,ber2,nloop);
%************************end of file***********************************