dvb_t_ce.m

基于DVB-T标准,COFDM调制系统的信道估计,仿真条件：4个OFDM符号

%###### 基于DVB-T标准,COFDM调制系统的信道估计的原始模型 ######
%######  仿真条件：4个OFDM符号，QPSK调制，二维内插，信道模型  ######
% IFFT_bin_length: IFFT和FFT的点数
% carrier_count: 子载波个数
% bits_per_symbol: 每个载波含有的比特数
% symbols_per_carrier: 不重复的OFDM符号数
% Ｘ：欲发送的二进制比特流

clear all;
close all;

IFFT_bin_length=2048;
carrier_count=1705;
bits_per_symbol=2;
symbols_per_carrier=4;
LI=12;                  % 导频之间的间隔
GI=2048/4;              % 保护间隔长度
N_number=carrier_count*symbols_per_carrier*bits_per_symbol;

%--------------------------------------------------------------------------
%矩阵初始化
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=1-X.*2;
sreal=s(1:2:N_number);
simage=s(2:2:N_number);
X1=sreal+j.*simage;

%--------------------------------------------------------------------------
%% 串并转换 %%
X2=reshape(X1,carrier_count,symbols_per_carrier).';

%--------------------------------------------------------------------------
%产生随机导频信号
Kmin=0;
Kmax=1704;
pilot(:,1)=Kmin;               %使第一列是导频
for l=0:symbols_per_carrier-1
    p=0;
    k=0;
    while k<=Kmax-12
       k=Kmin+3*mod(l,4)+12*p;
       pilot(l+1,p+2)=k;       %分散导频位置
       p=p+1;
    end
end
pilot(:,p+2)=Kmax;             %使最后一列也是导频
pilot=pilot+1;

g=bin2dec('100000000101');     %伪随机二进制序列生成多项式:x11+x2+1
state=bin2dec('11111111111');  %伪随机二进制序列生成寄存器初始状态
N=2^11-1;                      %生成二进制序列长度
train_sym=zeros(symbols_per_carrier,carrier_count);
m=mgen(g,state,N);             %调用mgen函数，生成二进制随机序列
for i=1:l+1
    train_sym(i,:)=m(1:carrier_count);
end
for i=1:l+1
    train_sym(i,pilot(i,:))=(4/3)*2.*(1/2-train_sym(i,pilot(i,:))); %导频值
end

%--------------------------------------------------------------------------
%% 插入导频 %%
signal=1:carrier_count;
X3(:,signal)=X2;
for i=1:l+1
    X3(i,pilot(i,:))=train_sym(i,pilot(i,:));
end

%--------------------------------------------------------------------------
%% IFFT变换 %%
IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);
IFFT_modulation(:,signal)=X3;
X4=ifft(IFFT_modulation,IFFT_bin_length,2);

%--------------------------------------------------------------------------
%% 加循环前缀保护间隔 %%
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));

%--------------------------------------------------------------------------
%信道模型:带多普勒频移的瑞利衰落信道
fd=200;                           %多普勒频移
r=6;                              %多径数
a=[0.123 0.3 0.4 0.5 0.7 0.8];    %多径的幅度
d=[20 30 40 50 90 130];           %各径的延迟
T=1;                              %系统采样周期
th=[90 0 72 144 216 288]*pi./180; %相移
h=zeros(1,carrier_count);
hh=[];
for k=1:r
    h1=a(k)*exp(j*((2*pi*T*fd*d(k)/carrier_count)));
    hh=[hh,h1];
end
h(d+1)=hh;                        %信道冲激响应
H0=fft(h);

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;

%--------------------------------------------------------------------------
%% 加高斯白噪声 %%  
    snr_db=15;                                          %信噪比
    code_power=[norm(Tx_data)]^2/(length(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).';
    
%--------------------------------------------------------------------------
%% 去保护间隔 %%
    for k=1:symbols_per_carrier
       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(:,signal);
  
%--------------------------------------------------------------------------
%% 信道估计 %%
    H=zeros(size(Y3));
    for i=1:l+1 
        H(i,pilot(i,:))=Y3(i,pilot(i,:))./X3(i,pilot(i,:));  %导频处信道响应
    end
    
%--------------------------------------------------------------------------
%% 内插 %%
    %导频间隔为12的位置处的内插
    for i=1:l+1     
        for j=pilot(i,2):12:(pilot(i,p+2)-12)
            for n=1:11
                H(i,j+n)=(1-n/12)*H(i,j)+n/12*H(i,j+12);
            end
        end
    end
    %导频间隔不为12的位置的内插
    i=2;
    for n=1:2
        H(i,pilot(i,1)+n)=(1-n/3)*H(i,j)+n/3*H(i,pilot(i,2));
    end
    for n=1:8
        H(i,pilot(i,p+1)+n)=(1-n/9)*H(i,pilot(i,p+1))+n/9*H(i,pilot(i,p+2));
    end
    
    i=3;
    for n=1:5
        H(i,pilot(i,1)+n)=(1-n/6)*H(i,j)+n/6*H(i,pilot(i,2));
    end
    for n=1:5
        H(i,pilot(i,p+1)+n)=(1-n/6)*H(i,pilot(i,p+1))+n/6*H(i,pilot(i,p+2));
    end
    
    i=4;
    for n=1:8
        H(i,pilot(i,1)+n)=(1-n/9)*H(i,j)+n/9*H(i,pilot(i,2));
    end
    for n=1:2
        H(i,pilot(i,p+1)+n)=(1-n/3)*H(i,pilot(i,p+1))+n/3*H(i,pilot(i,p+2));
    end
    
%--------------------------------------------------------------------------
%% 滤波 %%
    IFFT_H=ifft(H,carrier_count,2);
    h_filter=[IFFT_H(:,1:GI),zeros(symbols_per_carrier,carrier_count-GI)];
    H_ce=fft(h_filter,carrier_count,2);
    
%--------------------------------------------------------------------------
%% 均衡 %%
    Y2=Y3./H_ce;
             
%--------------------------------------------------------------------------
%% 并串变换 %%
    Y1=reshape(Y2.',1,carrier_count*symbols_per_carrier);
    Y0=reshape(Y3.',1,carrier_count*symbols_per_carrier);
    H1=reshape(H_ce.',1,carrier_count*symbols_per_carrier);
    H1=H1(1:carrier_count);%信道响应
    
%--------------------------------------------------------------------------
%% 作图 %%
    figure;
    plot(abs(H0),'r:');
    hold on;
    plot(abs(H1));
    hold on;
    legend('No Channel Estimation','Channel Estimation');
    title('信道估计前后传输函数比较');
    hold off;	
   
%--------------------------------------------------------------------------
%% QPSK解调 %%
    y_re1=sign(real(Y1));
    y_im1=sign(imag(Y1));
    
    y_re0=sign(real(Y0));
    y_im0=sign(imag(Y0));
 
    r1=zeros(1,N_number);
    r1(1:2:N_number)=y_re1;
    r1(2:2:N_number)=y_im1;
    
    r0=zeros(1,N_number);
    r0(1:2:N_number)=y_re0;
    r0(2:2:N_number)=y_im0;
    
    Y=(1-r01.*sqrt(2))./2;            %输出二进制序列

% %%%%% 导频位置及导频个数统计程序段 %%%%%
% % 得到每个OFDM符号的分散导频位置
% Kmin=0;
% Kmax=1704;
% for l=0:63
%     p=0;
%     k=0;
%     while k<=Kmax-12
%        k=Kmin+3*mod(l,4)+12*p;
%        m(l+1,p+1)=k;
%        p=p+1;
%     end
% end
% 
% % 比较分散导频位置和连续导频位置，计算重复位置个数
% n=0;
% pilot=[0 48 54 87 141 156 192 201 255 279 282 333 432 450 483 525 531 618 636 714 759 765 780 804 873 888 918 939 942 969 984 1050 1101 1107 1110 1137 1140 1146 1206 1269 1323 1377 1491 1683 1704];
% for j=1:45
%     pil=pilot(j);
%     for i=1:p
%         if m(2,i)==pil;
%             n=n+1;
%         end
%     end
% end