%ofdm.m
%QPSK仿真實現OFDM傳輸系統
%---------定義參數--------
para=128;     %并行傳輸的信道數(點數)
fftlen=128;    %FFT長度
noc=128;      %子載波數
nd=6;         %每次循環含有的OFDM 符號數
ml=2;         %調制電平:QPSK(二電平)
sr=4800;      %符號速率
br=sr.*ml;    %每載波的比特率
gilen=32;     %保護間隔(點數)
%------------主循環部分---------
nloop2=100;   %仿真循環次數
noe=0;        %誤碼率
nod=0;        %傳輸的數據數量
eop=0;        %誤組數
nop=0;        %傳輸的分組數
snr=0;        %信噪比

nloop1=input('nloop1='); 

for jjj=1:nloop1
    
    snr=snr+1;
    
  for iii=1:nloop2
    %--------------發送部分--------
    %--------------產生數據--------
    seldata=rand(1,para*ml)>0.5;%(均勻分布的隨機矩陣函數；para*ml每次產生的數據個數)
   %seldata得到的實際上是右邊那個邏輯表達式的結果,rand產生的值如果大于0.5,則seldata為1,如果小于0.5,則seldata為0
   %---------------------串并變換-------------
   paradata=reshape(seldata,para,ml);
   %變換方式:前128bit變為第一列,隨后的128bit為第二列,依此類推
   %----------------------QPSK調制-------------
   x=qpskmod(paradata,para,ml);  % 實現QPSK調制，并規一化  
   
   %-------------IFFT--------------
   y=ifft(x,fftlen) ;
   
   %-----------并串轉換－－－－－－－－
   paradata1=reshape(y,1,para);
   ich2=real(paradata1);
   qch2=imag(paradata1);
  
   %------------插入保護間隔----------
   [ich3,qch3]=giins(ich2,qch2,para,gilen);
   para1=para+gilen;
   paradata2=ich3+qch3.*i;
  
   %--------------加入加性高斯白噪聲--------------
   
   paradata3=awgn(paradata2,snr,'measured');
   %--------------接收部分-------------
   %-----------------移除保護間隔---------------
   ich4=real(paradata3);
   qch4=imag(paradata3);
   [ich5temp,qch5temp]=girem(ich4,qch4,para1,gilen);
  
   
   %-----------------串并轉換－－－－－－－－
   ich5=reshape(ich5temp,para,1);
   qch5=reshape(qch5temp,para,1);
   %----------------FFT----------------
   rx=ich5+qch5.*i;
   ry=fft(rx,fftlen);
   ich6=real(ry);
   qch6=imag(ry);
  
   %---------------QPSK解調-----------
   [demodata]=qpskdemod(ich6,qch6,para);
   %------------------并串轉換-------------
   demodata1=reshape(demodata,1,para*ml);
   %-------------------誤比特率(BER)-------------
   %即使的誤碼和數據
   noe2=sum(abs(demodata1-seldata));
   %將接受解調后的信號與原始信號相比較,累計不一樣的
   nod2=length(seldata);   %發送信號總長度
   %累計誤碼數和總的數據
   noe=noe+noe2;
   %因為是進行了nloop次的循環,所以把每次循環的數據和誤碼數積累起來
   nod=nod+nod2;
   %計算誤組率(PER)(每次循環作為一個分組)
   end
 ber=noe/nod %總的誤碼率
 noe=0;
 nod=0;
xsnr(jjj)=snr;
yber(jjj)=ber;
end 
%-----------調制后信號的實部，虛部圖-----------
   figure(1); 
   xr=real(x);
   xi=imag(x);
   subplot(2,1,1);stem([1:length(xr)],xr.');
   axis([0,30,-1,1]);
   title('QPSK實部圖形');
   subplot(2,1,2);stem([1:length(xi)],xi.');
   axis([0,30,-1,1]);
   title('QPSK虛部圖形');
     %------調制后的信號頻譜圖——————
   figure(2);
   plot([1:256],abs(fft(x,256)));
   axis([0 256 0 30]);
   title('QPSK信號頻譜圖');  
    %------FFT后的星座圖————————
   figure(3);
   ich6temp=ich6*sqrt(2);
   qch6temp=qch6*sqrt(2);
   plot(ry,'o');
   axis([-2,2,-2,2]);
   title('QPSK星座圖');
%------------------------結束---------------
   
   
  %--------------BPSK仿真實現OFDM傳輸系統-------------
                        
snr=0;

for jjj=1:nloop1
    
    snr=snr+1;
    
  for iii=1:nloop2
    %--------------發送部分--------
    %--------------產生數據--------
    seldata=rand(1,para)>0.5;%(均勻分布的隨機矩陣函數；para*ml每次產生的數據個數)
   %seldata得到的實際上是右邊那個邏輯表達式的結果,rand產生的值如果大于0.5,則seldata為1,如果小于0.5,則seldata為0
   %---------------------串并變換-------------
   paradata=reshape(seldata,para,1);
   %變換方式:前128bit變為第一列,隨后的128bit為第二列,依此類推
   %----------------------BPSK調制-------------
   x=bpskmod(paradata,para);  % 實現BPSK調制，并規一化  
   %-------------IFFT--------------
   y=ifft(x,fftlen) ;
   %-----------并串轉換－－－－－－－－
   paradata1=reshape(y,1,para);
   ich2=real(paradata1);
   qch2=imag(paradata1);
  
   %------------插入保護間隔----------
   [ich3,qch3]=giins(ich2,qch2,para,gilen);
   para1=para+gilen;
   paradata2=ich3+qch3.*i;
   
  
   %--------------加入加性高斯白噪聲--------------
   
   paradata3=awgn(paradata2,snr,'measured');
   %--------------接收部分-------------
   %-----------------移除保護間隔---------------
   ich4=real(paradata3);
   qch4=imag(paradata3);
   [ich5temp,qch5temp]=girem(ich4,qch4,para1,gilen);
  
%-----------------串并轉換－－－－－－－－
   ich5=reshape(ich5temp,para,1);
   qch5=reshape(qch5temp,para,1);
   %----------------FFT----------------
   rx=ich5+qch5.*i;
   ry=fft(rx,fftlen);
   ich6=real(ry);
   qch6=imag(ry);
  
   %---------------BPSK解調-----------
   [demodata]=bpskdemod(ich6,para);
   %------------------并串轉換-------------
   demodata1=reshape(demodata,1,para);
   %-------------------誤比特率(BER)-------------
   %即使的誤碼和數據
   noe2=sum(abs(demodata1-seldata));
   %將接受解調后的信號與原始信號相比較,累計不一樣的
   nod2=length(seldata);   %發送信號總長度
   %累計誤碼數和總的數據
   noe=noe+noe2;
   %因為是進行了nloop次的循環,所以把每次循環的數據和誤碼數積累起來
   nod=nod+nod2;
   %計算誤組率(PER)(每次循環作為一個分組)
   end
 ber=noe/nod %總的誤碼率
 noe=0;
 nod=0;
xsnr(jjj)=snr;
yber(jjj)=ber;
end 

%------星座圖————————
   figure(4);
   plot(ich6,qch6,'o');
   axis([-2,2,-2,2]);
   title('BPSK星座圖'); 
    %------BPSK信號頻譜圖——————
   figure(5);
   plot([1:256],abs(fft(x,256)));
   axis([0 256 0 30]);
   title('BPSK信號頻譜圖');  
%------------------------結束--------------- 
%仿真實現OFDM傳輸系統



%------------使用16QAM調制
snr=0;      %信噪比

for jjj=1:nloop1
    snr=snr+1; 
for iii=1:nloop2
%--------------發送部分--------
   %--------------產生數據--------
    BitsTx = floor(rand(1,para*4)*2);
    %---------------------串并變換-------------
    SymQAMtmp=zeros(para,4);
    aa=1;
    for ii=1:para
        SymQAMtmp(aa,1)=BitsTx(4*aa-3);
        SymQAMtmp(aa,2)=BitsTx(4*aa-2);
        SymQAMtmp(aa,3)=BitsTx(4*aa-1);
        SymQAMtmp(aa,4)=BitsTx(4*aa);
        aa=aa+1;
    end
    
    %----------------------16QAM調制-------------
    SymQAMtmptmp = bi2de(SymQAMtmp,2,'right-msb');
    x=qammod(SymQAMtmptmp,16);
    
   %-------------IFFT--------------
   y=ifft(x,fftlen) ;
   %-----------并串轉換－－－－－－－－
   paradata1=reshape(y,1,para);
   ich2=real(paradata1);
   qch2=imag(paradata1);
  
   %------------插入保護間隔----------
   [ich3,qch3]=giins(ich2,qch2,para,gilen);
   para1=para+gilen;
   paradata2=ich3+qch3.*i;
   
   %--------------加入加性高斯白噪聲--------------
   
   paradata3=awgn(paradata2,snr,'measured');

   %--------------接收部分-------------
   %-----------------移除保護間隔---------------
    ich4=real(paradata3);
   qch4=imag(paradata3);
   [ich5temp,qch5temp]=girem(ich4,qch4,para1,gilen);
   %------------串并轉換—————————
   ich5=reshape(ich5temp,para,1);
   qch5=reshape(qch5temp,para,1);
   %----------------FFT----------------
   rx=ich5+qch5.*i;
   ry=fft(rx,para);
   ich6=real(ry);
   qch6=real(ry);
   %---------------16QAM解調----------
   [demodata]=qamdemod(ry,16);
   %------------------并串轉換------------
   
   demodata1=de2bi(demodata).';
 demodata2=reshape(demodata1,1,para*4);
           
   
   %-------------------誤比特率(BER)-------------
   %即使的誤碼和數據
   noe2=sum(abs(demodata2-BitsTx));
   %將接受解調后的信號與原始信號相比較,累計不一樣的
   nod2=length(BitsTx);   %發送信號總長度
   %累計誤碼數和總的數據
   noe=noe+noe2;
   %因為是進行了nloop次的循環,所以把每次循環的數據和誤碼數積累起來
   nod=nod+nod2;
end
ber=noe/nod    %總的誤碼率
noe=0;
nod=0;
xsnr(jjj)=snr;
yber(jjj)=ber;
end
 %------IFFT后的信號頻譜圖——————
   figure(6);
   plot([1:256],abs(fft(x,256)));
   axis([0 256 0 100]);
   title('16-QAM信號頻譜圖');  
    %------FFT后的星座圖————————
   figure(7);

   plot(ry,'.');
   axis([-8,8,-8,8]);
   title('16-QAM星座圖');

%grid;