bpsk.m

基于梳状导频的ofdm系统LS信道分析

%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_symbols
X3(:,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);    
   end
end
%---------------------------------------------------------
%并串转换
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];
    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,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 on

title('信噪比和误码率曲线');
xlabel('信噪比');
ylabel('误码率');
legend('没做信道估计','LS信道估计');
hold off