%本程序采用迭代的聯合信道估計與符號檢測算法，仿真了不同信噪比條件下，%
%信道估計的性能，包括有迭代和無迭代條件下的誤碼率和均方誤差%
%------------------------------------------------------------------------%
echo off;clear all;
close all;
clc;
fprintf( 'OFDM仿真\n') ;
tic
% ---------------------------------------------%
%                   參數定義                    %
% --------------------------------------------- %
% Initialize the parameters
Num_Symbol=1000;  %仿真OFDM符號數
Bits_per_Symbol=60;  %每符號比特數
NumSubc=128;   %   載波數
Numcp=NumSubc/4;  %cp數約為載波數的1/4
mentor_times=20;  %蒙氏仿真  最大次數
% h_time=create_channel(6);%產生一個6徑信道
%  figure;
%  stem(abs(h_time));
%  CL=length(h_time);  %信道長度
pilots=1:NumSubc/16:NumSubc;  %導頻插入位置 共16個導頻  等間隔插入
PL=length(pilots);   %PL=16； 共16個導頻
%------------------------------------------------%
Bits_Tx = floor(rand(Num_Symbol,Bits_per_Symbol)*2); 
% Generate the random binary stream for transmit test
%-------------------------------------------------%
% 卷積編碼
for i=1:Num_Symbol
    a=Bits_Tx(i,:);
    Bits_convenc(i,:)=cnv(a);   
end
%------------------------------------------------%
%  隨機交織   用系統函數生成
state=32;  %定義交織器初始參數
for i=1:Num_Symbol
    b=Bits_convenc(i,:);
Bits_Interleaved(i,:)=randintrlv(b,state); %隨機交織
end
%---------------------------------------------------%
%  4QAM  星座圖映射  
  for i=1:Num_Symbol 
      c=Bits_Interleaved(i,:);
      QAM_modulated(i,:)=QAM_modu(c); %調用映射函數 
  end
  %---------------------------------------------------%
  %  在第一行插入導頻  作信道估計用    
   QAM_modulated(1,pilots)=1;%  導頻位置處插入導頻1+j
   QAM_modulated(1,pilots+1)=0;%   與導頻相鄰位置處插0
   %----------------------------------------------------%
%  串/并轉換
Symbol_paralleled=QAM_modulated.';
%--------------------------------------------------------%
 %  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 
 
% %-------------------------------------------------------%
%過信道
after_channel=zeros(m,n+Numcp);
 snrtable=zeros(11,3);  %存儲信噪比和誤碼率
 MSE_table=zeros(11,2); %存儲估計的均方誤差
 for snr=0:3:30
     ber_ratio1=0;  %初始化誤碼率
     ber_ratio2=0;  %初始化誤碼率  
     estimation_error=zeros(2,NumSubc); %用以存儲中間過程的估計誤差 
   snrtable(snr/3+1,1)=snr;%第一列存儲信噪比
   for mentor=1:mentor_times
       h_time=create_channel(6);%產生一個6徑信道
       CL=length(h_time);  %信道長度
   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_temp=after_channel.';  %轉置  出來是并聯的
 Symbol_de_cp=zeros(n,m);
 Symbol_de_cp(1:n,:)=Symbol_de_cp_temp((Numcp+1):(n+Numcp),:);  %去掉最上面的Numcp行
 %-----------------------------------------------------%
 %  送入fft解調器
 Symbol_fft=fft(Symbol_de_cp,NumSubc);
  Symbol_fft_temp=Symbol_fft.'; %轉置  變成串聯
 %-----------------------------------------------------%
 %------------------------------------------------------%
 % 信道估計部分
 H_act=fft([h_time zeros(1,NumSubc-CL)],NumSubc);
  H_est_temp=Symbol_fft_temp(1,pilots)./QAM_modulated(1,pilots);  %導頻位置處的信道初始估計
   h1=ifft(H_est_temp,PL);  %導頻點數的ifft
 h2=h1(1:CL);  %取前CL點
 h3=[h2 zeros(1,NumSubc-CL)];  %補零  達到NumSubc位
 H_est1=fft(h3,NumSubc);  %NumSubc點fft的變換  到頻域
 
%------------------------------------------------------%
%---------------------------------------------------------%
 for i=1:Num_Symbol-1  
     X(i,:)= Symbol_fft_temp(i+1,:)./H_est1;  %接收信號除以估計的信道頻率響應
 end
 for i=1:Num_Symbol-1
     b=X(i,:);
     transmit_bits_demapped(i,:)=QAM_demodu(b);  %解映射
 end
  for i=1:Num_Symbol-1
     c=transmit_bits_demapped(i,:);
 transmit_bits_deinterleaved(i,:)=randdeintrlv(c,state);  %解交織
  end
  for i=1:Num_Symbol-1
      e=transmit_bits_deinterleaved(i,:);
      [transmit_bits_decoded,survivor_state,cumulated_metric]=viterbi(e);  %譯碼
      decode(i,:)=transmit_bits_decoded;
  end
  %--------------------------------------------
  Bits_Tx_temp=Bits_Tx(2:Num_Symbol,:);
  [number,ratio] = biterr(Bits_Tx_temp,decode);   %計算誤碼率
  ber_ratio1=ber_ratio1+ratio;  %累計誤碼率
  %--------------------------------------------------------%
 %---------------------------------------------------------%
 
 
 
 
   error1=(abs(H_act(1,:)-H_est1(1,:)).^2);  %初始估計誤差
  estimation_error(1,:)=estimation_error(1,:)+error1;  %累計誤差
 
 RX=Symbol_fft_temp(2,:); %用接收的第二行數據做測試數據用
 X_est1=RX./H_est1;  %接收到的數據信號除以信道頻率響應為發送信號的估計值
 x_demap1=QAM_demodu(X_est1);  %解調
 
 x_deinterleaved=randdeintrlv(x_demap1,state);  %解交織
 [x_decoded,survivor_state,cumulated_metric]=viterbi(x_deinterleaved);  %譯碼
 
 x_coded=cnv(x_decoded);  %編碼
 x_interleaved=randintrlv(x_coded,state);%交織
 x_mapped=QAM_modu(x_interleaved);  %映射
 H_est2_temp1=RX./x_mapped;  
 H_est2_temp2=ifft(H_est2_temp1,NumSubc);
 H_est2=fft([H_est2_temp2(1,1:CL) zeros(1,NumSubc-CL)],NumSubc);
 X_est2=RX./H_est2;
 x_demap2=QAM_demodu(X_est2);
 %-----------------------------------------
 %迭代信道估計部分
 while(x_demap2~=x_demap1)  %判斷兩次估計值的是否相同，相同則表示已收斂，無需繼續迭代，否則繼續迭代
     x_deinterleaved=randdeintrlv(x_demap2,state);  %解交織
     [x_decoded,survivor_state,cumulated_metric]=viterbi(x_deinterleaved);  %譯碼
     x_coded=cnv(x_decoded);  %編碼
     x_interleaved=randintrlv(x_coded,state);%交織
     x_mapped=QAM_modu(x_interleaved);  %映射
     %以下三行是估計
     H_est2_temp1=RX./x_mapped;  
     H_est2_temp2=ifft(H_est2_temp1,NumSubc);  %變換到時域
     H_est2=fft([H_est2_temp2(1,1:CL) zeros(1,NumSubc-CL)],NumSubc);%補零后再變換到頻域
     
     X_est1=RX./H_est2;  %接收信號除以信道響應得到發送信號的估計值
     x_demap1=QAM_demodu(X_est1);  %解映射
     %以下三行在交換，用以滿足繼續迭代的條件
     temp=x_demap1;
     x_demap1=x_demap2;
     x_demap2=temp;
 end
   H_est=H_est2; %得到估計值   
 
%-------------------------------------------------------------%
%---------------------------------------------------------------%
      %迭代后的接收判決部分

 for i=1:Num_Symbol-1  
     X(i,:)= Symbol_fft_temp(i+1,:)./H_est;  %接收信號除以估計的信道頻率響應
 end
 for i=1:Num_Symbol-1
     b=X(i,:);
     transmit_bits_demapped(i,:)=QAM_demodu(b);  %解映射
 end
  for i=1:Num_Symbol-1
     c=transmit_bits_demapped(i,:);
 transmit_bits_deinterleaved(i,:)=randdeintrlv(c,state);  %解交織
  end
  for i=1:Num_Symbol-1
      e=transmit_bits_deinterleaved(i,:);
      [transmit_bits_decoded,survivor_state,cumulated_metric]=viterbi(e);  %譯碼
      decode(i,:)=transmit_bits_decoded;
  end
  %--------------------------------------------
  Bits_Tx_temp=Bits_Tx(2:Num_Symbol,:);
  [number,ratio] = biterr(Bits_Tx_temp,decode);   %計算誤碼率
  ber_ratio2=ber_ratio2+ratio;  %累計誤碼率
  error2=(abs(H_act(1,:)-H_est(1,:)).^2);  %估計誤差
  estimation_error(2,:)=estimation_error(2,:)+error2;  %累計誤差
  %------------------------------ -----------------------
  %------------------------------------------------------%
   end
   snrtable(snr/3+1,2)=ber_ratio1/mentor_times;
   snrtable(snr/3+1,3)=ber_ratio2/mentor_times;
   estimation_error=estimation_error./mentor_times;
    MSE_table(snr/3+1,1)=sum(estimation_error(1,:))/NumSubc; %賦值
    MSE_table(snr/3+1,2)=sum(estimation_error(2,:))/NumSubc; %賦值
    %-------------------------------------------------------------%
%     figure;
%  plot(abs(H_act));
%  hold on;
%  plot(abs(H_est));
%  hold off;
 end
  snrtable
  MSE_table
 figure;
semilogy(snrtable(:,1),snrtable(:,2),'k*-');
hold on;
semilogy(snrtable(:,1),snrtable(:,3),'ko-');
legend('無迭代','有迭代');
grid on;
xlabel('SNR/dB');
ylabel('BER');
title('誤碼率/信噪比曲線');
hold off;

figure;
semilogy(snrtable(:,1),MSE_table(:,1),'k*-');
hold on;
semilogy(snrtable(:,1),MSE_table(:,2),'ko-');
grid on;
legend('無迭代','有迭代');
xlabel('SNR/dB');
ylabel('MSE');
title('均方誤差/信噪比曲線');
hold off;

time_of_sim = toc
echo on;