%OFDM Channel Estimation Based on Comb Pilot%IFFT_bin_length: IFFT和FFT的點數%carrier_count: 子載波個數%bits_per_symbol: 每符號上的比特數%symbols_per_carrier: 每楨的OFDM符號數%Ｘ：欲發送的二進制比特流clear all;clc;IFFT_bin_length=600;%ifft長度carrier_count=512;%子載波數bits_per_symbol=2;%每符號比特數symbols_per_carrier=6;%一楨符號數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=128%8;             % guard interval lengthN_snr=40;          % 每比特信噪比snr=8;             %信噪比間隔%------------------------------------------------------------% vector initializationX=zeros(1,N_number);%2400個bitX1=[];X2=[];X3=[];X4=[];X5=[];X6=[];X7=[];Y1=[];Y2=[];Y3=[];Y4=[];Y5=[];Y6=[];Y7=[];XX=zeros(1,N_number);%2400dif_bit=zeros(1,N_number);%2400dif_bit1=zeros(1,N_number);%2400dif_bit2=zeros(1,N_number);%2400dif_bit3=zeros(1,N_number);%2400X=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*24t=(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*12training_symbols2=training_symbols1.';%12*1training_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  113if length(pilot)~=Np    pilot=[pilot,carrier_count+Np];%最后一列變成導頻1     9    17    25    33    41    49    57    65    73    81    89    97   105  113   116end%--------------------------------------------------------%串并轉換X2=reshape(X1,carrier_count,symbols_per_carrier).';%12*100，12個復信號符號，100列載波%---------------------------------------------------------%插入導頻signal=1:carrier_count+Np;%1*116signal(pilot)=[];%1*100X3(:,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*128IFFT_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);       endend%---------------------------------------------------------%并串轉換X7=reshape(X6.',1,symbols_per_carrier*(IFFT_bin_length+GI));%12*136先轉置，再變成1*1632復信號流%---------------------------------------------------------%信道模型:帶多普勒頻移的瑞利衰落信道fd=0; %多普勒頻移r=6;   %多徑數a=[0.123 0.3 0.4 0.5 0.7 0.8]; %多徑的幅度d=[2*16 3*16 4*16 5*16 9*16 13*16]; %各徑的延遲% a=[0 0 0 0 0 0 ]; %多徑的幅度% d=[0 0 0 0 0 0]; %各徑的延遲T=1;  %系統采樣周期%th=[90 0 72 144 216 288]*pi./180;%相移h=zeros(1,carrier_count);%1*100hh=[];    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*1632channel1(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+channel5+channel6;%4徑干擾后的數據流1*1632%---------------------------------------------------------------   tyb=awgn(Tx_data,20,'measured');%%%%%%%555555%corrcoef(x,y)for i=1:para*11    c=0;    for k=1:4    for j=1:gl       a(k)=ReData(i+j-1+(k-1)*(para+gl)).*conj(ReData(i+j+para-1+(k-1)*(para+gl)));      r(i)=r(i)+a(k)/2;       c=r(i);    end      c=r(i+k-1)+c;      end    r(i)=c/2;    b(i)=abs(r(i));endr(1:IFFT_bin_length*5)=0;b(1:IFFT_bin_length*5)=0;a(1:4)=0;%corrcoef(x,y)for i=1:IFFT_bin_length*5    for j=1:GI        a=tyb(i+j-1).*conj(tyb(i+j+IFFT_bin_length-1));        r(i)=r(i)+a;    end    b(i)=abs(r(i));end  x=1:IFFT_bin_length*5;figureplot(x,b)xlabel('采樣點(k)');ylabel('循環前綴的相關值Rxy(k)')%---------------------------------------------------------------%----------------------------------------------------------------%加高斯白噪聲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;     Y7=Tx_data;     %-------------------------------------------------------  %串并變換   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];    endelse    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];endBER=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