% 本程序用以仿真塊狀導頻時不同信噪比條件下的誤碼率和均方誤差

echo off;clear all;
close all;
clc;
fprintf( 'OFDM仿真\n') ;
tic
% ---------------------------------------------%
%                   參數定義                    %
% --------------------------------------------- %
% Initialize the parameters
Num_frame=1000;  %仿真OFDM符號數
Bits_per_frame=60;  %每符號比特數
NumSubc=128;   %   載波數
Numcp=NumSubc/4;  %cp數約為載波數的1/4
mentor_times=20;  %蒙氏仿真  最大次數
%------------------------------------------------%
Bits_Tx = floor(rand(Num_frame,Bits_per_frame)*2); 
% Generate the random binary stream for transmit test
%------------------------------------------------%
% 信道編碼   卷積碼  對符號(即每行)進行編碼
for i=1:Num_frame
    a=Bits_Tx(i,:);
    Bits_convenc(i,:)=cnv(a);   
end
%------------------------------------------------%
%  隨機交織  對每符號進行交織   用系統函數生成
state=32;  %定義交織器初始參數
for i=1:Num_frame
    b=Bits_convenc(i,:);
Bits_Interleaved(i,:)=randintrlv(b,state); %隨機交織
end
%-------------------------------------------------%
%  4QAM  星座圖映射  對每幀進行映射
  for i=1:Num_frame 
      c=Bits_Interleaved(i,:);
      QAM_modulated(i,:)=QAM_modu(c); %調用映射函數 
  end
%----------------------------------------------------%
%  插入導頻序列
[m,n]=size(QAM_modulated);
Pilot_seq=ones(1,n);  %導頻序列 
Symbol_tx=[Pilot_seq;QAM_modulated];  %將其插入信息符號的前部， 成為發送符號 
%----------------------------------------------------%
%  串/并轉換
 Symbol_paralleled=Symbol_tx.';
 %---------------------------------------------------%
 %  IFFT變換
 Symbol_ifft_temp=ifft( Symbol_paralleled,NumSubc);  %已變換
 %---------------------------------------------------%
 %  并/串轉換
 Symbol_ifft= Symbol_ifft_temp.';
 %----------------------------------------------------%
 %  加cp前綴 
 [m,n]=size(Symbol_ifft);
 
 Symbol_cp=zeros(m,n+Numcp);
 Symbol_cp(:,(Numcp+1):(n+Numcp))=Symbol_ifft(:,1:n);%先把Symbol_ifft整體復制到Symbol_cp_temp的后半部分去
 Symbol_cp(:,1:Numcp)=Symbol_ifft(:,(n-Numcp+1):n);  %加cp  

 %-----------------------------------------------------%
%------------------------------------------------------%
 %  送入信道 
 h_time=[0.2901+0.6129j 0.4061+0.3881j 0.2103+0.1637j];  %假設的信道沖擊響應
 L=length(h_time);  %信道長度
 H_act=fft([h_time zeros(1,NumSubc-L)],NumSubc);  %信道頻率響應
 
 %-------------------------------------------------------%
 after_channel=zeros(m,n+Numcp);
 %  加高白噪聲
 snrtable=zeros(11,2);
 MSE_table=zeros(11,1);
 for snr=0:2:20
     ber_ratio=0;  %初始化誤碼率
   snrtable(snr/2+1,1)=snr;
   estimation_error=zeros(1,NumSubc); %用以存儲中間過程的估計誤差 
   for mentor=1:mentor_times
 for i=1:m
     d=Symbol_cp(i,:);
 after_channel_temp=filter(h_time,[1],d);  %   過信道  用濾波器實現卷積
 after_channel(i,:)=awgn(after_channel_temp,snr,'measured');  %加噪聲
 end
 %------------------------------------------------------% 
 %-----------------------------------------------------%
 % 去cp前綴
 Symbol_de_cp_temp1=after_channel.';  %串/并轉換
 Symbol_de_cp=zeros(n,m);
 Symbol_de_cp(1:n,:)=Symbol_de_cp_temp1((Numcp+1):(n+Numcp),:);  %即要把Symbol_de_cp_temp1最上面的NumCP行去掉 變成n*m矩陣 已經作好串/并轉換
 %-----------------------------------------------------%
 %  送入fft解調器
 Symbol_fft=fft(Symbol_de_cp,NumSubc);  %出來是并聯的
 %-----------------------------------------------------%
 % 信道估計部分 
 Y1_Pilot_seq=Symbol_fft(:,1);  %接收到的導頻序列
 X1_Pilot_seq=Pilot_seq(1,:).'; %這是發送端的導頻序列
 H_p=zeros(NumSubc,1); 
 H_p(:)=Y1_Pilot_seq./X1_Pilot_seq; %由導頻得到的初始信道估計
  H_temp1=ifft(H_p,NumSubc);  %反變換到時域
  H_temp2=H_temp1.';
  H_temp3=[H_temp2(:,1:L) zeros(1,NumSubc-L)];  %補零  
  H_est=fft(H_temp3,NumSubc).';  %變換到頻域
  error=(abs(H_act-H_est.').^2);  %估計誤差
  estimation_error(1,:)=estimation_error(1,:)+error;  %累計誤差
%--------------------------------------------------------%
 for i=1:Num_frame     
     X(:,i)= Symbol_fft(:,i+1)./H_est;  %得到估計值
 end
 for i=1:Num_frame
     b=X(:,i).';
     transmit_bits_demapped(i,:)=QAM_demodu(b);  %解調
 end
 for i=1:Num_frame
     c=transmit_bits_demapped(i,:);
 transmit_bits_deinterleaved(i,:)=randdeintrlv(c,state);  %解交織
 end
  for i=1:Num_frame
      e=transmit_bits_deinterleaved(i,:);
      [transmit_bits_decoded,survivor_state,cumulated_metric]=viterbi(e);  %譯碼
      decode(i,:)=transmit_bits_decoded;
  end
  %-----------------------------------------------%
  %  誤碼率計算
  [number,ratio] = biterr(Bits_Tx,decode);
ber_ratio=ber_ratio+ratio;
   end
 snrtable(snr/2+1,2)=ber_ratio/mentor_times;
 estimation_error=estimation_error./mentor_times;
 MSE_table(snr/2+1,1)=sum(estimation_error(1,:))/NumSubc; %賦值
 
%---------------------------------
%各信噪比條件下估計值與真實值的比較
  figure(snr/2+1);
 plot(abs(H_act),'k-');
 grid on;
hold on;
 plot(abs(H_est),'k-.');
 legend('真實響應','估計的響應');
 title('擬合效果');
 hold off;
%-----------------------------------
 end

snrtable
 MSE_table
%%-----------------------------------------------%
figure;
semilogy(snrtable(:,1),snrtable(:,2),'k*-');
grid on;
xlabel('SNR/dB');
ylabel('BER');
title('信噪比/誤碼率曲線')

figure;
semilogy(snrtable(:,1),MSE_table(:,1),'k*-');
xlabel('SNR/dB');
ylabel('MSE');
title('信噪比/均方誤差曲線');
grid on;
time_of_sim = toc
echo on;