📄 zf_mmse_bj.m
字号:
%## 基于国家数字电视地面广播（DTTB）标准,单载波调制系统的信道估计的原始模型 ##%
%################## 单载波应用模式：帧头模式2，帧体模式C=1 ###################%
%############ 仿真条件：PN595帧头，16QAM调制，多径衰落信道模型 #############%
%######################## 比较ZF和MMSE对多径的抑制作用 #####################%
clear all;
close all;
len_head=595; %帧头长度（帧头模式2）
len_carriers=3744; %有效信息符号长度
len_sysmes=36; %系统信息符号长度
len_body=len_carriers+len_sysmes; %帧体长度
len_frame=len_head+len_body; %信号帧长度
symbols=8;
fd=2000; %采样间隔（2000Hz）
bits=4; %帧体每符号上的比特数
fs=7.56*1e+6; %基带采样率（7.56MHz）
ts=1/fs; %基带符号周期（1/7.56us）
%----------------------------------TRANSMISSION---------------------------%
%--------------------------------------------------------------------------
%% 发送数据调制 %%
N=symbols*len_carriers*bits; %发送有效数据包含比特数
rand('state',0); %确定随机发生器种子
sendbits=round(rand(1,N)); %产生数据位
%% 符号映射 %%
data_qam=fun_GB_QAM(sendbits,bits); %调用4QAM调制函数
trans_data=reshape(data_qam,len_carriers,symbols).'; %串并变换，每一行代表一个信息符号
%--------------------------------------------------------------------------
%% 系统信息调制 %%
sysmes=zeros(1,len_sysmes); %36位系统信息
sysmes(1:4)=[0 0 0 0]; %4位系统信息：子载波个数C=1模式
bits_sysmes='00010001101000010100011110111100'; %32位系统信息：16QAM以及LDPC码率1，交织模式1
symmes(5:36)=bits_sysmes-48;
sysmes_bits=zeros(1,2*len_sysmes); %进行I,Q路数据相同的4QAM调制
sysmes_bits(1:2:end)=sysmes;
sysmes_bits(2:2:end)=sysmes;
sysmes_qam=fun_GB_QAM(sysmes_bits,2); %进行I,Q路数据相同的4QAM调制
trans_sysmes=repmat(sysmes_qam,symbols,1); %表示symbols个符号的系统信息的矩阵
%--------------------------------------------------------------------------
%% 帧体数据处理 %%
trans_body=[trans_sysmes trans_data];
frame_day=[zeros(symbols,len_head) trans_body];
len=size(frame_day);
frame_day_serial=reshape(frame_day.',1,len(1)*len(2)); %并串变换
% frame_day_serial(1:2:end)=frame_day_serial(1:2:end)-1; %偶数符号，实部-1（第一个符号看作第0个符号，为偶数）
% frame_day_serial(2:2:end)=frame_day_serial(2:2:end)+1; %奇数符号，实部+1（第二个符号看作第1个符号，为奇数）
frame_day_parallel=reshape(frame_day_serial.',len(2),len(1)).'; %串并变换
body_process=frame_day_parallel(:,len_head+1:end); %处理完成后的帧体数据
%--------------------------------------------------------------------------
%% 已知帧头信号 %%
head_PN595=zeros(symbols,len_head); %PN595序列生成
g=bin2dec('10010000001'); %伪随机二进制序列生成多项式:1+x3+x10
state=bin2dec('0000000001'); %伪随机二进制序列生成寄存器初始状态
N=2^10-1; %生成二进制序列长度
m=mgen(g,state,N); %完成帧头0到+1，1到-1的映射
PN595=m(1:len_head); %帧头数据取m序列的前595个数值
for i=1:len_head %调用mgen函数，生成1023位二进制随机序列
if PN595(i)==0
PN595(i)=-4.5*(1+1j);
else
PN595(i)=4.5*(1+1j);
end
end
head_PN595=repmat(PN595,symbols,1); %组成symbols个帧头
body_power=mean(mean(abs(body_process).^2,2))/(4.5^2*2); %帧体数据平均功率
head_power=body_power; %帧头数据的平均功率和帧体相同
head_PN595=head_PN595.*sqrt(head_power); %帧头发送数据
%--------------------------------------------------------------------------
%% 组帧 %%
trans_signal=[head_PN595 body_process]; %将每个符号的帧头和帧体组成一个完整信号帧
%--------------------------------------------------------------------------
%% 并串变换 %%
trans_frame=reshape(trans_signal.',1,symbols*len_frame); %并串变换，变换为串行数据
SNR=15; %信噪比
frame_passchannel = awgn(trans_frame,SNR,'measured');
trans_frame=frame_passchannel;
%-----------------------------------PASS CHANNEL--------------------------%
%--------------------------------------------------------------------------
%% 信道模型:带多普勒频移的多径衰落信道 %%
X7=trans_frame;
r=3; %多径数
a=[0.2 0.3 0.4]; %多径的幅度
d=[30 50 80]; %各径的延迟
channel1=zeros(size(X7));
channel1(1+d(1):length(X7))=a(1)*X7(1:length(X7)-d(1));
channel2=zeros(size(X7));
channel2(1+d(2):length(X7))=a(2)*X7(1:length(X7)-d(2));
channel3=zeros(size(X7));
channel3(1+d(3):length(X7))=a(3)*X7(1:length(X7)-d(3));
Tx_data=X7+channel1+channel2+channel3;
%--------------------------------------------------------------------------
%% 加高斯白噪声 %%
SNR=15; %信噪比
frame_passchannel = awgn(Tx_data,SNR,'measured');
%------------------------------------RECEPTION----------------------------%
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
%% 串并变换 %%
receiv_frame=reshape(frame_passchannel,len_frame,symbols); %并串变换，将变换为串行数据
receiv_frame=receiv_frame.';
%--------------------------------------------------------------------------
%% 去帧头PN序列 %%
receiv_head=receiv_frame(:,1:len_head); %接收到的帧头序列
receiv_body=receiv_frame(:,len_head+1:end); %接收到的帧体数据
%--------------------------------------------------------------------------
%% FFT %%
head_FFT=fft(receiv_head,len_head,2);
body_FFT=fft(receiv_body,len_body,2);
%--------------------------------------------------------------------------
%% 信道估计 %%
PN595_FFT=fft(head_PN595,len_head,2);
H_head=head_FFT./PN595_FFT; %由帧头的PN序列估计信道
%--------------------------------------------------------------------------
%% ZF频域均衡 %%
h_head=ifft(H_head,len_head,2);
H_body=fft(h_head,len_body,2); %所得信道估计
W_zf=1./H_body;
body_equa_zf=body_FFT.*W_zf; %迫零均衡（ZF）
body_zf=ifft(body_equa_zf,len_body,2);
%--------------------------------------------------------------------------
%% MMSE频域均衡 %%
H_mmse_fm=abs(H_body).^2+1/10^(SNR/10);
W_mmse=conj(H_body)./H_mmse_fm;
body_equa_mmse=body_FFT.*W_mmse; %最小均方误差均衡（MMSE）
body_mmse=ifft(body_equa_mmse,len_body,2);
%--------------------------------------------------------------------------
%% 自相关图 %%
figure(1); %发射信号自相关图
subplot(221);
Nc=len_frame;
cy=xcorr(trans_frame(596:4375),Nc); %生成长度为2*Nc+1的自相关序列
tt=(-Nc*ts:ts:Nc*ts)*1e6;
plot(tt,20*log10(abs(cy)/abs(cy(Nc+1))));
axis([-20 20 -40 0]);
zoom on;
xlabel('时间（us)');
ylabel('幅度（dB）');
title('发射数据自相关');
subplot(222); %接收信号自相关图
Nc=len_frame;
cyy=xcorr(receiv_body(1,:),Nc); %生成长度为2*Nc+1的自相关序列
tt=(-Nc*ts:ts:Nc*ts)*1e6;
plot(tt,20*log10(abs(cyy)/abs(cyy(Nc+1))));
axis([-20 20 -40 0]);
zoom on;
xlabel('时间（us)');
ylabel('幅度（dB）');
title('通过信道后数据自相关');
subplot(223); %ZF均衡后数据自相关图
CY=xcorr(body_zf(1,:),Nc); %生成长度为2*Nc+1的自相关序列
tt=(-Nc*ts:ts:Nc*ts)*1e6;
plot(tt,20*log10(abs(CY)/abs(CY(Nc+1))));
axis([-20 20 -40 0]);
zoom on;
xlabel('时间（us)');
ylabel('幅度（dB）');
title('ZF均衡后数据自相关');
subplot(224); %MMSE均衡后数据自相关图
CYY=xcorr(body_mmse(1,:),Nc); %生成长度为2*Nc+1的自相关序列
tt=(-Nc*ts:ts:Nc*ts)*1e6;
plot(tt,20*log10(abs(CYY)/abs(CYY(Nc+1)))); %均衡后数据自相关图
axis([-20 20 -40 0]);
zoom on;
xlabel('时间（us)');
ylabel('幅度（dB）');
title('MMSE均衡后数据自相关');
% %--------------------------------------------------------------------------
% %% 计算副峰抑制值 %%
% a=20*log10(abs(cyy)/abs(cyy(Nc+1)));
% b=20*log10(abs(CY)/abs(CY(Nc+1)));
% c=20*log10(abs(CYY)/abs(CYY(Nc+1)));
%
% zf1=a(30+Nc+1)-b(30+Nc+1);
% zf2=a(50+Nc+1)-b(50+Nc+1);
% zf3=a(80+Nc+1)-b(80+Nc+1);
%
% mmse1=a(30+Nc+1)-c(30+Nc+1);
% mmse2=a(50+Nc+1)-c(50+Nc+1);
% mmse3=a(80+Nc+1)-c(80+Nc+1);
%--------------------------------------------------------------------------
%% MMSE均衡自相关图 %%
figure(3);
subplot(211);
cyy=xcorr(receiv_body(1,:),Nc); %生成长度为2*Nc+1的自相关序列
tt=(-Nc*ts:ts:Nc*ts)*1e6;
plot(tt,20*log10(abs(cyy)/abs(cyy(Nc+1)))); %接收数据自相关图
axis([-20 20 -40 0]);
zoom on;
xlabel('时间（us)');
ylabel('幅度（dB）');
title('通过信道后数据自相关');
subplot(212);
CYY=xcorr(body_mmse(1,:),Nc); %生成长度为2*Nc+1的自相关序列
tt=(-Nc*ts:ts:Nc*ts)*1e6;
plot(tt,20*log10(abs(CYY)/abs(CYY(Nc+1)))); %均衡后数据自相关图
axis([-20 20 -40 0]);
zoom on;
xlabel('时间（us)');
ylabel('幅度（dB）');
title('MMSE均衡后数据自相关');
%--------------------------------------------------------------------------
%% MMSE均衡互相关图 %%
figure(5);
Nc=595;
x7=Tx_data(596:4375);
CYY=xcorr(x7,Nc); %生成长度为2*Nc+1的自相关序列
tt=(-Nc*ts:ts:Nc*ts)*1e6;
plot(tt,20*log10(abs(CYY)/max(abs(CYY))));
hold on
cyy2=xcorr(x7,body_mmse(1,:),Nc);
plot(tt,20*log10(abs(cyy2)/max(abs(cyy2))),'r'); %接收数据自相关图
xlabel('时间（us)');
ylabel('幅度（dB）');
hold on;
legend('均衡前数据自相关','均衡后数据互相关');
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -