dvb_t_tongbu_1.m

基于DVB-T标准,COFDM调制系统的利用导频信号进行符号粗同步,仿真条件：68个OFDM符号

%###### 基于DVB-T标准,COFDM调制系统的利用导频信号进行符号粗同步 ######
%######  仿真条件：68个OFDM符号，64QAM调制，莱斯信道  ######
% 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=6;
symbols_per_carrier=74; % 总的OFDM符号数(4=10-3-3)
LI=12;                  % 导频之间的间隔
GI=2048/4;              % 保护间隔长度
N_number=carrier_count*symbols_per_carrier*bits_per_symbol;

%--------------------------------------------------------------------------
%矩阵初始化
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）

%--------------------------------------------------------------------------
%% 64QAM调制 %%
[sreal simage]=QAM_64(X);        %调用64QAM调制函数
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=zeros(symbols_per_carrier,carrier_count);
train_sym=zeros(symbols_per_carrier,carrier_count);
for i=1:l+1
    m=mgen(g,state,N);                      %调用mgen函数，生成二进制随机序列
    train(i,:)=m(1:carrier_count);
end

for i=1:l+1
    train_sym(i,pilot(i,:))=(4/3)*2.*(1/2-train(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

%% 插入连续导频 %%
CP=[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]+1;                    %45个连续导频位置

for i=1:length(CP)
    train_sym(:,CP(i))=(4/3)*2.*(1/2-train(:,CP(i)));%连续导频值
end
for i=1:length(CP)
    X3(:,CP(i))=train_sym(:,CP(i));                  %插入连续导频
end

%% 插入传数参数信令 %%
TPS=[34 50 209 346 413 569 595 688 790 901 1073 1219 1262 1286 1469 1594 1687]+1; %17个传数参数信令位置

for i=1:length(TPS)
    train_sym(:,TPS(i))=randint;
    train_sym(:,TPS(i))=2*(1/2-train_sym(:,TPS(i))); %传数参数信令值 
end
for i=1:length(TPS)
    X3(:,TPS(i))=train_sym(:,TPS(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));

% %--------------------------------------------------------------------------
% %% 加入基于DVB-T的莱斯多径衰落信道(固定接受模式) %%
% [Tx_data h]=dvb_channel(X7);     %调用多径加入函数
Tx_data=X7;

%--------------------------------------------------------------------------
%% 加高斯白噪声 %%  
snr_db=15;                                          %信噪比
X2=X2(4:symbols_per_carrier-3,:);                   %取中间连续四个符号
    
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;


N=IFFT_bin_length;
Ng=GI;
r=Y7(20001:30000);           %任意选取一段数据进行估计
% r=Y7;

%--------------------------------------------------------------------------
% 最大似然法
N=IFFT_bin_length;
G=GI;
T=1:3*N+2*G;
for i=1:length(T)            % FFT窗移动的距离
    data1=r(i:G+i-1);        % 取GI长个数据存data1
    data2=r(N+i:N+G+i-1);    % 取相距IFFT_bin_length的GI长个数据存data2
    s(i)=data1*data2';       % 求互相关值    
    cor(i)=abs(s(i)); 
%     cor(i)=abs(real(data1))*(abs(real(data2)))'+abs(imag(data1))*(abs(imag(data2)))';
end

%--------------------------------------------------------------------------
% 作图
subplot(211);
plot(T,cor);
xlabel('载波数');
ylabel('相关值');
grid on;

[Max Max_i]=max(cor(1:N));                  % Max_i对应一个符号的第一个数据
start_place=Max_i;
for i=1:length(T)
    deltaF(i)=1/(2*pi)*angle(conj(s(i)));
    r1=r(i:G+i-1);
    r2=r(i+N:G+i+N-1);
    e(i)=1/(2*pi)*atan(imag(r1*r2')/(real(r1*r2')));
end
% deltaF=1/(2*pi)*angle(conj(s(Max_i)));      % 分数频偏估计值
% r1=r(Max_i:G+Max_i-1);
% r2=r(Max_i+N:G+Max_i+N-1);
% e=1/(2*pi)*atan(imag(r1*r2')/(real(r1*r2')));
% disp(start_place);disp(e);
subplot(212);  
plot(T,e);
xlabel('载波数');
ylabel('粗频偏估计值');
axis(1,length(T),-0.5,0.5);
grid on;