基于DVB-T标准,COFDM调制系统的利用导频信号进行符号粗同步

原创

fpga和matlab 2022-10-10 15:58:07 博主文章分类：MATLAB ©著作权

文章标签 DVB-T COFDM调制系统数据信令插值 文章分类 虚拟化云计算

©著作权归作者所有：来自51CTO博客作者fpga和matlab的原创作品，请联系作者获取转载授权，否则将追究法律责任

%###### 基于DVB-T标准,COFDM调制系统的利用导频信号进行符号粗同步 ######
%###### 仿真条件：68个OFDM符号，64QAM调制，莱斯信道 ######
% IFFT_bin_length: IFFT和FFT的点数
% carrier_count: 子载波个数
% bits_per_symbol: 每个载波含有的比特数
% symbols_per_carrier: 参与估计的总的OFDM符号数
% Ｘ：欲发送的二进制比特流

clear all;
close all;
rng(2);
IFFT_bin_length=2048;
carrier_count=1705;
bits_per_symbol=4;
symbols_per_carrier=68; % 总的OFDM符号数
LI=12; % 导频之间的间隔
GI=2048/4; % 保护间隔长度
N_number=carrier_count*symbols_per_carrier*bits_per_symbol;
% 格雷码特性设定
EbNo = 4.5:.5:7; linEbNo = 10.^(EbNo(:).*0.1);
ModulateIndex = 16; codeRate = 1/2; constlen = 7; k = log2(ModulateIndex); codegen = [171 133];
tblen = 32; % traceback length
trellis = poly2trellis(constlen, codegen);
dspec = distspec(trellis, 7);
expVitBER = bercoding(EbNo, 'conv', 'hard', codeRate, dspec);
%--------------------------------------------------------------------------
%矩阵初始化
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=rand(1,N_number); %产生二进制随机序列（非0即1）
X=randint(carrier_count*symbols_per_carrier,1,16);
%--------------------------------------------------------------------------
%% 16QAM调制 %%
numSymb = 1512*4*2;
hStr = RandStream('mt19937ar', 'Seed', 654321);
msg_orig = randi(hStr, [0 1], numSymb, 1);
msg_enc = convenc(msg_orig, trellis);
hMod = modem.qammod('M', ModulateIndex, 'PhaseOffset', 0, ...
'SymbolOrder', 'Gray', 'InputType', 'Bit');
msg_tx = modulate(hMod, msg_enc);

X1=msg_tx;

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

%--------------------------------------------------------------------------
%产生随机导频信号
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;

% 16QAM的平均幅度
A_avg=(3*sqrt(2)*4+sqrt(2)*4+sqrt(10)*8)/16;
Burst=1*A_avg*4/3;
%--------------------------------------------------------------------------
%% 插入分散导频 %%
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=round(rand(1,carrier_count)); %改成随机数，不用上列PRBS生成法。简化
train(i,:)=m(1:carrier_count);
end

for i=1:l+1
train_sym(i,pilot(i,:))=Burst*2.*(1/2-train(i,pilot(i,:))); %分散导频值
end
signal=1:carrier_count;
X3(:,signal)=0;
for i=1:l+1
X3(i,pilot(i,:))=train_sym(i,pilot(i,:)); %插入分散导频
end
X3_SPCP=X3(1:4,:); %保留原始插入分散导频
X3_SPCP(2:4,1)=0;X3_SPCP(2:4,1705)=0;
ScPilotX=X3(1:4,:);
%% 插入连续导频 %%
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))=Burst*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))=A_avg*2*(1/2-train_sym(:,TPS(i))); %传数参数信令值
end
for i=1:length(TPS)
X3(:,TPS(i))=train_sym(:,TPS(i)); %插入传数参数信令
end
%--------------------------------------------------------------------------
%% 插入数据
Data_index=zeros(4,1705);
X_data=X2;
X_out=X3(1:4,:);
for i=1:4
m=1;
for j1=1:1705
if abs(X_out(i,j1))<0.1
X_out(i,j1)=X_data(i,m);
Data_index(i,j1)=3; % 记忆有效数据点处为3，导频和TPS为0，为接收端提取数据用。
m=m+1; %只有1512个数据，最后一个m多加了1.
end
end
M(i)=m;
end
Data_index=[Data_index;Data_index;Data_index;Data_index];
%组成68个OFDM符号，4个FODM符号为一循环
X41=[X_out;X_out;X_out;X_out;X_out;X_out;X_out;X_out];
X42=[X41;X41;X_out];
CP_pilot=X42(1:6,:);
% ------------------------------------------------------------------------
%% IFFT变换 %%
IFFT_modulation=zeros(symbols_per_carrier,IFFT_bin_length);
IFFT_modulation(:,signal)=X42;
X4=ifft(IFFT_modulation,IFFT_bin_length,2); %2为行运算，1为列运算。
% 下面的方法与前一行运算等同。复数矩阵行列调换按如下方式进行，不能用reshape
%ifft_x=(IFFT_modulation).';
%X41=ifft(ifft_x,2048);
%X4=(X41).';

%--------------------------------------------------------------------------
%% 加循环前缀保护间隔 %%
X10=zeros(68,2560);
for j1=1:68
X10(j1,1:GI)=X4(j1,2048-512+1:end);
X10(j1,GI+1:end)=X4(j1,1:end);
end

%--------------------------------------------------------------------------
%% 并串转换 %% reshape按列重排，X6是68x2560，必须先转成为2560x68，再重排
X7=reshape(X10.',1,symbols_per_carrier*(IFFT_bin_length+GI));

% %--------------------------------------------------------------------------
% %% 加入基于DVB-T的瑞利多径衰落信道(固定接受模式) %%
% 加衰减信道

chan = rayleighchan; % Static Rician channel
chan.PathDelays=[0 0.05 0.4 1.45 2.3 2.8]*1e-6; % 长时延
chan.AvgPathGaindB=[-2.8 0 -3.8 -0.1 -2.6 -1.3];
%chan.PathDelays=[0 5 14 35 54 75]*1e-6; % 长时延
%chan.AvgPathGaindB=[0 -9 -22 -25 -27 -28];
%chan.PathDelays=[0 0.05 0.4 1.45 2.3 2.8]*1e-6; %短时延
%chan.AvgPathGaindB=[-2.8 0 -3.8 -0.1 -2.6 -1.3];
chan.InputSamplePeriod=1e-7;
fadedSig = filter(chan, X7); % Apply channel effects
Tx_data = awgn(fadedSig,16,'measured'); % Add Gaussian noise.

Y7=[zeros(1,400) Tx_data(2560*3+1:end)]; % Tx_data(2560*2+1:end)] **********************************************************************信道参数

%% ----------最大似然法，确定粗定时偏差 --------------------------
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

% ---------------极大似然作图----------------------------------------------
figure(1)
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;
% Max_i=400;
r=r(Max_i-30:Max_i-30+12*2560);

%% -------------------------- FFT变换 ---------------------------
N=IFFT_bin_length;
Ng=GI;
Y8=r; % 加延时，使截取点落入循环前缀之内。
r=Y8;
r_ofdmin=reshape(r(1:12*2560),2560,12).';
r_nocp=r_ofdmin(:,GI+1:end);
r_Sym=fft(r_nocp,2048,2);
%figure(2)
scatterplot(r_Sym(1:1705));
title('存在定时误差的FFT窗口取样符号');

%% ------------------细定时估计------------------------------------
% 找出分散导频位置
P=12;
Y_1=r_Sym(1,1:end);
Y_5=r_Sym(5,1:end);
% 找出导频位置点，逐点相关。
for i=1:P*4;
cc=0;
for p=0:142
cc=cc+Y_1(12*p+i)*conj(Y_5(12*p+i));
end
c(i)=abs(cc);
end
% 另一种找出导频点的方法，更好且简单，确定四种导频方式其中的一种。
for k1=0:3
dd=0;
for p=1:142
dd=dd+Y_1(12*p+1+3*k1)*conj(Y_5(12*p+1+3*k1));
end
d(k1+1)=abs(dd);
end

[cMax Max_ip]=max(d(1:4));
figure(3)
subplot(2,1,1)
plot(c);
subplot(2,1,2)
stem(d);
% 确定细定时偏差，假定定时偏差在（-85，+85）之间。
Y_SP=r_Sym(1,:); %取第1个OFDM符号。
sum=0;K1=-N/(2*pi*P);
SP_X=ScPilotX(Max_ip,:); % 取出模式Max_ip指定的分散导频
SP_X=[SP_X, zeros(1,2*P)]; %补点，防止下标溢出。
SP_start=3*(Max_ip-1)+1;
for i=0:142
sum=sum+conj(Y_SP(i*P+SP_start))*(SP_X(i*P+SP_start))...
*conj(conj(Y_SP((i+1)*P+SP_start))*(SP_X((i+1)*P+SP_start)));
%sum=sum+(Y_SP(i*P+SP_start))*conj(Y_SP((i+1)*P+SP_start));
end
err_avg=K1*angle(sum)/1;
Time_err=round(err_avg);
%Time_err=200; %保证取样点正确，为星座图便于理解的例外设定

% 修正定时误差产生的旋转量
k=0:2048-1;
X_modify=r_Sym(1:12,:);
% Time_err=-16;
for i1=1:12
X_modify1(i1,:)=X_modify(i1,:).*exp(j*2*pi*(Time_err)*k/N);
end
%figure(4)
scatterplot(X_modify1(Max_ip,:));
title('经定时误差修正的符号');
% ---------------------------------------------------------

%% -------------------- 信道估计 --------------------------
% 简化方式，只能对第一个OFDM符号做估计，且认为是分散导频模式0.
% 只保留分散导频和连续导频
r_chestimation=X_modify1(:,1:1705);
if (Max_ip==1)
First_ip=1;
end
if (Max_ip==2)
First_ip=4;
end
if (Max_ip==3)
First_ip=3;
end
if (Max_ip==4)
First_ip=2;
end
if (Max_ip>4)
First_ip=1;
end
% First_ip=4;
r_chestimation=X_modify1(First_ip:First_ip+8-1,:);
X_modify2=r_chestimation;

r_chestimation(:,TPS)=0;
for m=1:8
for k=1:1705
if (abs(Data_index(m,k))>0.5)
r_chestimation(m,k)=0;
end
end
end
r_chestimation_sum=zeros(1,1705); % *****************************************************************
X3_SPCP_sum=zeros(1,1705);
X3_SPCP12=[X3_SPCP;X3_SPCP;X3_SPCP];
CP_pilot12=[CP_pilot(1:4,:); CP_pilot(1:4,:); CP_pilot(1:4,:)];
X3_SPCP12(:,CP)=CP_pilot12(:,CP); % 发射端导频，数据载波为0；
% 将四个符号的导频合并到一个符号中
%for m=1:4
% r_chestimation_sum=r_chestimation_sum+r_chestimation(m,1:1705);
% X3_SPCP_sum=X3_SPCP_sum+X3_SPCP(m,1:1705);
%end

% 连续导频用第一个符号的连续导频置换。修改为Max_ip置换。

%CP_pilot
for i=1:8
for k=1:1705;
if(abs(X3_SPCP12(i,k))>0.001)
Hp(i,k)=r_chestimation(i,k)./X3_SPCP12(i,k);
end
end
end
Hp_CP(:,CP)=Hp(:,CP); %保留连续导频的信道估计待用

% 插值
% 列方向（时间）插值
i=1;
for n=0:568-1
while (i<=5)
Hp(i+1,3*n+1)=(Hp(i,3*n+1)+3*Hp(i+4,3*n+1))/4;
Hp(i+2,3*n+1)=(Hp(i,3*n+1)+Hp(i+4,3*n+1))/2;
Hp(i+3,3*n+1)=(3*Hp(i,3*n+1)+Hp(i+4,3*n+1))/4;
break
end
i=i+1;
if i==5
i=1;
end
end
% 第一列连续导频再替换回去。
Hp1=Hp;
Hp1(:,CP)=Hp_CP(:,CP);
% 行方向（频率）插值，对第5个OFDM符号的信道估计
for n=0:568-1
Hp1(5,3*n+2)=2*Hp1(5,3*n+1)/3+Hp1(5,3*n+4)/3;
Hp1(5,3*n+3)=Hp1(5,3*n+1)/3+2*Hp1(5,3*n+4)/3;
end
% Hp1=Hp;

%for k=0:568-1;
% Hp(1,3*k+2)=(Hp(1,3*k+1)*2/3)+(Hp(1,3*(k+1)+1))*1/3;
% Hp(1,3*k+3)=(Hp(1,3*k+1)*1/3)+(Hp(1,3*(k+1)+1))*2/3;
%end
% 信道估计修正接收的符号，第5个符号
% X_modify2(5,1:1705)=X_modify1(5,1:1705)./Hp1(5,1:1705);
X_modify3(5,1:1705)=X_modify2(5,1:1705)./Hp1(5,1:1705);
% figure(5)
scatterplot(X_modify3(5,1:1705));
title('经信道估计修正的符号');
xlim([-5 5]); ylim([-5 5]);
% --------------统计超过幅度的符号点--------------------------------
recod=zeros(2,1705);sum1=0;
for k=1:1705
if (abs(real(X_modify2(1,k)))>7)||(abs(imag(X_modify2(1,k)))>7)
recod(1,k)=X_modify2(1,k);
recod(2,k)=k;
sum1=sum1+1;
end
end

%% -------------------------------------------------------------

% 从前1705个子载波中提取1512个数据符号，去除导频和TPS。
k2=1;
for k=1:1705
if (abs(Data_index(5,k))>0.001)
S_data(k2)=X_modify3(5,k);
k2=k2+1;
end
end
%figure(6)
scatterplot(S_data);
title('抽取的数据符号');
xlim([-5 5]); ylim([-5 5]);

%% -------------解调和viterbi解码 --------------------------------
hDemod = modem.qamdemod('M', ModulateIndex, 'PhaseOffset', 0, ...
'SymbolOrder', 'Gray', 'OutputType', 'Bit');
% scatterplot(msg_rx_int);
msg_demod = demodulate(hDemod, S_data.');
msg_dec = vitdec(msg_demod, trellis, tblen, 'cont', 'hard');
[nChnlErrs BERChnl] = biterr(msg_enc(1:end/4), msg_demod);
[nCodErrs BERCoded] = biterr(msg_orig(1:end/4-tblen), msg_dec(1+tblen:end));
Time_err
Max_ip
Max_i
nChnlErrs
nCodErrs

系统误码率仿真结果如下：

基于DVB-T标准,COFDM调制系统的利用导频信号进行符号粗同步_COFDM调制系统