%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=128;carrier_count=100;bits_per_symbol=2;symbols_per_carrier=12;LI=7 ; %導頻之間的間隔Np=ceil(carrier_count/LI)+1;%導頻數 %加1的原因：使最后一列也是導頻N_number=carrier_count*symbols_per_carrier*bits_per_symbol;%（插入導頻后的）？？/總比特數carriers=1:carrier_count+Np;%？？？GI=8;             % guard interval length，保護間隔長度N_snr=40;          % 每比特信噪比snr=4;             %信噪比間隔？？%------------------------------------------------------------% vector initialization數據初始化X=zeros(1,N_number);X1=[];X2=[];X3=[];X4=[];X5=[];X6=[];X7=[];Y1=[];Y2=[];Y3=[];Y4=[];Y5=[];Y6=[];Y7=[];XX=zeros(1,N_number);dif_bit=zeros(1,N_number);dif_bit1=zeros(1,N_number);dif_bit2=zeros(1,N_number);dif_bit3=zeros(1,N_number);X=randint(1,N_number);%產生二進制隨即序列（非0即1）%--------------------------------------------------------%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;%？？？？產生1200個數據%bpsk調制for i=1:length(X)       if  X(i)==0           mod_out(i)=-1;       else mod_out(i)=1;%完成星座圖的映射 (mapping)       end    end%---------------------------------------------------------%產生隨機導頻信號%--------------------------------------------------------train_sym=randint(1,4*symbols_per_carrier);t=(train_sym.*2-1)/sqrt(2);treal=t(1:2:4*symbols_per_carrier);timage=t(2:2:4*symbols_per_carrier);training_symbols1=treal+j.*timage;training_symbols2=training_symbols1.';training_symbols=repmat(training_symbols2,1,Np);%產生的16個導頻信號%disp(training_symbols)pilot=1:LI+1:carrier_count+Np;if length(pilot)~=Np    pilot=[pilot,carrier_count+Np];%16個導頻在子100個子載波上的插入位置（共116，子載波100外加16個導頻）end%--------------------------------------------------------%串并轉換X2=reshape(mod_out,carrier_count,4*symbols_per_carrier/2).';%---------------------------------------------------------%插入導頻signal=1:carrier_count+Np;signal(pilot)=[];X3(:,pilot)=training_symbols;%把X3的第polit列全部賦值為training_symbolsX3(:,signal)=X2;%X3=cat(1,training_symbols,X2);%IFFT快速傅立葉逆變換IFFT_modulation=zeros(4*symbols_per_carrier/2,IFFT_bin_length);IFFT_modulation(:,carriers)=X3;%IFFT_modulation(:,conjugate_carriers)=conj(X3);X4=ifft(IFFT_modulation,IFFT_bin_length,2);%X5=X4.';%加保護間隔（循環前綴）for k=1:4*symbols_per_carrier/2;   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,4*symbols_per_carrier/2*(IFFT_bin_length+GI));%---------------------------------------------------------%信道模型:帶多普勒頻移的瑞利衰落信道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);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;%noise=randn(1,length(X7))+j.*randn(1,length(X7)); %--------------------------------------------------------channel1=zeros(size(X7));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;%---------------------------------------------------------------%---------------------------------------------------------------%----------------------------------------------------------------%加高斯白噪聲Error_ber=[];%誤比特率Error_ber1=[];Error_ber2=[];%誤比特率Error_ber3=[];%Error_ser=[];%誤符號率BER=zeros(1,length(0:snr:N_snr));BER1=zeros(1,length(0:snr:N_snr));j=1;for snr_db=0:snr:N_snr    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,4*symbols_per_carrier/2).';     %去保護間隔    for k=1:4*symbols_per_carrier/2;       for i=1:IFFT_bin_length;           Y5(k,i)=Y6(k,i+GI);       end    end    %FFT，傅立葉變換     Y4=fft(Y5,IFFT_bin_length,2);     Y3=Y4(:,carriers); %-------------------------------------------------------------    %LS信道估計,%信道估計準則：最小均方誤差（MMSE）最大似然估計（MLE） 最小平方（LS）  H=[];  Y2=Y3(:,signal);  Rx_training_symbols=Y3(:,pilot);  Rx_training_symbols0=reshape(Rx_training_symbols,4*symbols_per_carrier/2*Np,1);    training_symbol0=reshape(training_symbols,1,4*symbols_per_carrier/2*Np);  training_symbol1=diag(training_symbol0);  %disp(training_symbols)  training_symbol2=inv(training_symbol1);  Hls=training_symbol2*Rx_training_symbols0;  Hls1=reshape(Hls,4*symbols_per_carrier/2,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,4*symbols_per_carrier,1);%將導頻擴展成4*symbols_per_carrier*carrier_count矩陣  Y1=Y2./HLs;      %-------------------------------------------------------------  %并串變換  YY=reshape(Y2.',1,N_number);%沒有做信道估計  YY1=reshape(Y1.',1,N_number);%做了信道估計       %BPSK modulation    demod_out1=zeros(1,length(YY1));    demod_out=zeros(1,length(YY));     for i=1:length(YY1)       if abs( YY1(i)-1)-abs(YY1(i)+1)<=0    %demod_in(i)-1:接收點到1的距離;    demod_in(i)+1:接收點到-1的距離          demod_out1(i)=1;       else           demod_out1(i)=0;       end     end         for i=1:length(YY)       if abs( YY(i)-1)-abs(YY(i)+1)<=0    %demod_in(i)-1:接收點到1的距離;    demod_in(i)+1:接收點到-1的距離          demod_out(i)=1;       else           demod_out(i)=0;       end     end    %-------------------------------------------------------------%計算在不同信噪比下的誤比特率并作圖;%沒做信道估計BER(j)=length(find((X(1:length(demod_out)) == demod_out) == 0 ))/N_number;%沒做信道估計BER1(j)=length(find((X(1:length(demod_out1)) == demod_out1) == 0 ))/N_number;%做了信道估計j=j+1;end%------------------------------------------------------------- i=0:snr:N_snr;semilogy(i,BER,'-*r');%沒做信道估計hold on;semilogy(i,BER1,'-*g');%做了信道估計hold on;grid ontitle('信噪比和誤碼率曲線');xlabel('信噪比');ylabel('誤碼率');legend('沒做信道估計','LS信道估計');hold off