%基于梳妝導頻的OFDM信道估計 
clear all; 
clc; 
IFFT_bin_length=128;%ifft長度 
carrier_count=100;%子載波數 
bits_per_symbol=2;%每符號比特數 
symbols_per_carrier=12;%一楨符號數 
LI=7 ; %導頻之間的間隔 
Np=ceil(carrier_count/LI)+1;%16 導頻數加1的原因：使最后一列也是導頻 
N_number=carrier_count*symbols_per_carrier*bits_per_symbol;%2400一禎比特數 
carriers=1:carrier_count+Np;%116 子載波加導頻 
GI=8;             % guard interval length 
N_snr=40;          % 每比特信噪比 
snr=8;             %信噪比間隔 
%------------------------------------------------------------ 
% 初始化 
X=zeros(1,N_number);%2400個bit 
X1=[]; 
X2=[]; 
X3=[]; 
X4=[]; 
X5=[]; 
X6=[]; 
X7=[]; 
Y1=[]; 
Y2=[]; 
Y3=[]; 
Y4=[]; 
Y5=[]; 
Y6=[]; 
Y7=[]; 
XX=zeros(1,N_number);%2400 
dif_bit=zeros(1,N_number);%2400 
dif_bit1=zeros(1,N_number);%2400 
dif_bit2=zeros(1,N_number);%2400 
dif_bit3=zeros(1,N_number);%2400 
X=randint(1,N_number);%產生二進制隨即序列（非0即1）2400 
%-------------------------------------------------------- 
%QPSK調制：(1 1)->pi/4;(0 1)->3*pi/4;(0 0)->-3*pi/4;(1,0)->-pi/4; 
 
s=(X.*2-1)/sqrt(2); 
sreal=s(1:2:N_number); 
simage=s(2:2:N_number); 
X1=sreal+j.*simage;%已調信號bit流0.7071 - 0.7071i   0.7071 - 0.7071i   0.7071 + 0.7071i。。。。。（1*1200） 
 
%--------------------------------------------------------- 
%產生隨機導頻信號 
%-------------------------------------------------------- 
train_sym=randint(1,2*symbols_per_carrier);%1*24 
t=(train_sym.*2-1)/sqrt(2); 
treal=t(1:2:2*symbols_per_carrier); 
timage=t(2:2:2*symbols_per_carrier); 
training_symbols1=treal+j.*timage;% 0.7071 - 0.7071i  -0.7071 - 0.7071i  -0.7071 - 0.7071i   1*12 
training_symbols2=training_symbols1.';%12*1 
training_symbols=repmat(training_symbols2,1,Np);%12*16 復制第一列變成16列 
 
%disp(training_symbols) 
pilot=1:LI+1:carrier_count+Np;%導頻插入位置序號1     9    17    25    33    41    49    57    65    73    81    89    97   105  113 
if length(pilot)~=Np 
    pilot=[pilot,carrier_count+Np];%最后一列變成導頻1     9    17    25    33    41    49    57    65    73    81    89    97   105  113   116 
end 
%-------------------------------------------------------- 
%串并轉換 
X2=reshape(X1,carrier_count,symbols_per_carrier).';%12*100，12個復信號符號，100列載波 
%--------------------------------------------------------- 
%插入導頻 
signal=1:carrier_count+Np;%1*116 
signal(pilot)=[];%1*100 
X3(:,pilot)=training_symbols;%先放入16列導頻 
X3(:,signal)=X2;%再放入12*100，100列子載波，共12*116 
%X3=cat(1,training_symbols,X2); 
IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);%12*128 
IFFT_modulation(:,carriers)=X3;%116列后補12列全0子載波，12*128 
%IFFT_modulation(:,conjugate_carriers)=conj(X3); 
X4=ifft(IFFT_modulation,IFFT_bin_length,2);%每個符號128點ifft 
%X5=X4.'; 
%加保護間隔（循環前綴） 
for k=1:symbols_per_carrier; 
   for i=1:IFFT_bin_length; 
      X6(k,i+GI)=X4(k,i); 
   end 
   for i=1:GI; 
      X6(k,i)=X4(k,i+IFFT_bin_length-GI);     
   end 
end 
%--------------------------------------------------------- 
%并串轉換 
X7=reshape(X6.',1,symbols_per_carrier*(IFFT_bin_length+GI));%12*136先轉置，再變成1*1632復信號流 
 
%--------------------------------------------------------- 
%信道模型:帶多普勒頻移的瑞利衰落信道 
fd=100; %多普勒頻移 
r=6;   %多徑數 
a=[0.123 0.3 0.4 0.5 0.7 0.8]; %多徑的幅度 
d=[2 3 4 5 9 13]; %各徑的延遲 
T=1;  %系統采樣周期 
th=[90 0 72 144 216 288]*pi./180;%相移 
h=zeros(1,carrier_count);%1*100 
hh=[]; 
    for k=1:r 
        %deta=[zeros(1,d(k)-1),1,zeros(1,carrier_count-d(k))]; 
        h1=a(k)*exp(j*((2*pi*T*fd*d(k)/carrier_count))); 
        %h1=a(k)*exp(j*((2*pi*T*fd*d(k)/carrier_count))); 
        hh=[hh,h1]; 
    end 
    h(d+1)=hh;%3 4 5 6 10 14處有多徑效應 
%noise=randn(1,length(X7))+j.*randn(1,length(X7));  
%-------------------------------------------------------- 
 
channel1=zeros(size(X7));%1*1632 
channel1(1+d(1):length(X7))=hh(1)*X7(1:length(X7)-d(1)); 
channel2=zeros(size(X7)); 
channel2(1+d(2):length(X7))=hh(2)*X7(1:length(X7)-d(2)); 
channel3=zeros(size(X7)); 
channel3(1+d(3):length(X7))=hh(3)*X7(1:length(X7)-d(3)); 
channel4=zeros(size(X7)); 
channel4(1+d(4):length(X7))=hh(4)*X7(1:length(X7)-d(4)); 
channel5=zeros(size(X7)); 
channel5(1+d(5):length(X7))=hh(5)*X7(1:length(X7)-d(5)); 
channel6=zeros(size(X7)); 
channel6(1+d(6):length(X7))=hh(6)*X7(1:length(X7)-d(6)); 
%--------------------------------------------------------------- 
Tx_data=X7+channel1+channel2+channel3+channel4;%4徑干擾后的數據流1*1632 
%---------------------------------------------------------------- 
%加高斯白噪聲 
Error_ber=[];%誤比特率 
Error_ber1=[]; 
Error_ber2=[];%誤比特率 
Error_ber3=[]; 
%Error_ser=[];%誤符號率 
for snr_db=0:snr:N_snr %0:8:40 
 
    code_power=0; 
    code_power=[norm(Tx_data)]^2/(length(Tx_data));%信號的符號功率 
    %bit_power=var(Tx_data); 
    bit_power=code_power/bits_per_symbol;%比特功率  
    noise_power=10*log10((bit_power/(10^(snr_db/10))));%噪聲功率 
    noise=wgn(1,length(Tx_data),noise_power,'complex');%產生GAUSS白噪聲信號 
    Y7=Tx_data+noise; 
%------------------------------------------------------- 
  %串并變換 
   Y6=reshape(Y7,IFFT_bin_length+GI,symbols_per_carrier).';%先變成136*12，再轉置成12*136，恢復成行為符號，列為載波   
  %去保護間隔 
    for k=1:symbols_per_carrier; 
       for i=1:IFFT_bin_length; 
           Y5(k,i)=Y6(k,i+GI); 
       end 
    end 
     Y4=fft(Y5,IFFT_bin_length,2);%每行的符號進行128點fft  12*128 
     Y3=Y4(:,carriers);%去掉尾部12列原補零點， 12*116 
 %-------------------------------------------------------------    
 %LS信道估計 
  H=[]; 
  Y2=Y3(:,signal); 
  Rx_training_symbols=Y3(:,pilot); 
  Rx_training_symbols0=reshape(Rx_training_symbols,symbols_per_carrier*Np,1); 
   
  training_symbol0=reshape(training_symbols,1,symbols_per_carrier*Np); 
  training_symbol1=diag(training_symbol0); 
  %disp(training_symbols) 
  training_symbol2=inv(training_symbol1); 
  Hls=training_symbol2*Rx_training_symbols0; 
  Hls1=reshape(Hls,symbols_per_carrier,Np); 
  HLs=[]; 
  HLs1=[]; 
 if ceil(carrier_count/LI)==carrier_count/LI 
     for k=1:Np-1 
        HLs2=[]; 
           for t=1:LI 
           HLs1(:,1)=(Hls1(:,k+1)-Hls1(:,k))*(t-1)./LI+Hls1(:,k); 
           HLs2=[HLs2,HLs1]; 
        end 
       HLs=[HLs,HLs2]; 
    end 
else 
    for k=1:Np-2 
        HLs2=[]; 
        for t=1:LI 
           HLs1(:,1)=(Hls1(:,k+1)-Hls1(:,k))*(t-1)./LI+Hls1(:,k); 
           HLs2=[HLs2,HLs1]; 
        end 
       HLs=[HLs,HLs2]; 
    end 
    HLs3=[]; 
    for t=1:mod(carrier_count,LI) 
        HLs1(:,1)=(Hls1(:,Np)-Hls1(:,Np-1))*(t-1)./LI+Hls1(:,Np-1); 
        HLs3=[HLs3,HLs1]; 
    end; 
    HLs=[HLs,HLs3]; 
end 
  %Hls1=Hls.'; 
  %H=repmat(Hls1,symbols_per_carrier,1);%將導頻擴展成symbols_per_carrier*carrier_count矩陣 
  Y1=Y2./HLs;  
%------------------------------------------------------------------- 
  %并串變換 
  YY=reshape(Y2.',1,N_number/bits_per_symbol); 
  YY1=reshape(Y1.',1,N_number/bits_per_symbol); 
%------------------------------------------------------------ 
%QPSK解調 
   y_real=sign(real(YY)); 
   y_image=sign(imag(YY)); 
   y_re=y_real./sqrt(2); 
   y_im=y_image./sqrt(2);  
   y_real1=sign(real(YY1)); 
   y_image1=sign(imag(YY1)); 
   y_re1=y_real1./sqrt(2); 
   y_im1=y_image1./sqrt(2);  
  
   r00=[]; 
   r01=[]; 
   r10=[]; 
   r11=[]; 
  
  for k=1:length(y_real); 
     r00=[r00,[y_real(k),y_image(k)]]; 
  end; 
  for k=1:length(y_real1); 
     r10=[r10,[y_real1(k),y_image1(k)]]; 
  end; 
  
 for k=1:length(y_re); 
     r01=[r01,[y_re(k),y_im(k)]]; 
 end; 
 for k=1:length(y_re1); 
     r11=[r11,[y_re1(k),y_im1(k)]]; 
 end; 
 
    XX(find(r01>0))=1; 
%------------------------------------------------------------- 
%計算在不同信噪比下的誤比特率并作圖 
 
 dif_bit=s-r01;  
 dif_bit1=s-r11;  
 
 ber_snr=0;    %紀錄誤比特數 
    for k=1:N_number; 
       if dif_bit(k)~=0; 
         ber_snr=ber_snr+1; 
       end 
   end; 
 ber_snr1=0;    %紀錄誤比特數 
    for k=1:N_number; 
       if dif_bit1(k)~=0; 
          ber_snr1=ber_snr1+1; 
      end 
    end 
  
 Error_ber=[Error_ber,ber_snr]; 
 Error_ber1=[Error_ber1,ber_snr1]; 
end 
 
BER=zeros(1,length(0:snr:N_snr)); 
BER1=zeros(1,length(0:snr:N_snr)); 
 
BER=Error_ber./N_number; 
BER1=Error_ber1./N_number; 
%------------------------------------------------------------- 
%------------------------------------------------------------- 
 i=0:snr:N_snr; 
semilogy(i,BER,'-*r'); 
 
hold on; 
semilogy(i,BER1,'-og'); 
hold on; 
 
grid on; 
 
legend('No Channel Estimation','LS Channel Estimation'); 
hold off