模拟幅度调制相干解调系统抗噪声性能仿真分析

1、引言

1.1 研究目的

        \space\space\space\space\space\space\space        在模拟通信系统中产生的信号（如人声等），其信号的频率往往很低，在20-2000Hz左右。因此模拟信号的波长过长，不利于信号传输。通过信号调制技术，可以将有用信号乘以载波，由此得到高频信号，以此增强信道的传输效率。由此诞生了信号调制解调与信道抗噪性能分析等技术。

1.2 研究方法

        \space\space\space\space\space\space\space        在本文中，我们将采用理论计算与MATLAB仿真的研究方法建立信道模型，通过快速傅里叶变换与傅里叶反变换同时分析信号在时域和频域的状态。利用随机过程与高斯分布产生加性高斯白噪声，并设计不同种类的FIR滤波器进行噪声抑制，其中包括切比雪夫滤波器、FIR低通滤波器等。

1.3 主要内容

此次实验的主要内容将以如下三个模块展开：
1、抑制载波双边带调制DSB-SC
2、常规调幅调制AM
3、抑制载波单边带调制SSC
就三种不同信道的调制解调、滤波技术、抗噪性能分析等指标进行分析，包括理论分析与仿真结果的对比，理论与实际误差的产生原因分析等。

2、系统模型

2.1 DSB-SC

图2.1-1 DSB-SC信道模型

抑制载波双边带调制DSB-SC的时域调制模型为
$$s(t)=A_{c}m(t)cos2\pi f_{c}t\text{\hspace{1em}(2-1)}$$
其中携带信息的调制信号$m(t)$用单音信号模拟，其信号频率为$f_m$
$$m(t)=cos2\pi f_{m}t\text{\hspace{1em}(2-2)}$$

因此，通过傅里叶变换，得到$s(t)、m(t)$的频域表达式
$$\begin{gathered} S(f)=\frac{1}{4}[\delta (f-110)+\delta (f-90) \\ +\delta (f+110)+\delta (f+90)]\text{\hspace{1em}(2-3)} \end{gathered}$$
$$M(f)=\frac{1}{4}[\delta (f-10)+\delta (f+10)]\text{\hspace{1em}(2-4)}$$
其中，$f_m=10Hz、f_c=100Hz$

已调信号$s(t)$通过叠加信道高斯白噪声后通过带通滤波器BPF，该带通滤波器的时域和频率可以分别用抽样信号与双边矩形信号来表示

图2.1-2 BPF时域表示

图2.1-3 BPF频域表示

通过BPF后的信道可以滤除有用信号在信道中被叠加的高斯白噪声，再通过解调器还原原始的信号波形。已知解调器由乘法器与低通滤波器组成，其中乘法器模型表示为：
$$s_d(t)=s(t)cos2\pi f_{c}t\text{\hspace{1em}(2-5)}$$
最后通过低通滤波器LPF得到原始信号$m_o(t)$，低通滤波器的时域和频域表示如下：

图2.1-4 LPF时域表示

图2.1-5 LPF频域表示

2.2 AM

常规调幅(AM)的调制模型为

图2.2-1 AM调制信道模型

常规调幅(AM)信号的时域表达式为
$$s_{AM}(t)=A_c[1+a\cdot m_n(t)]cos2\pi f_{c}t\text{\hspace{1em}(2-6)}$$
其中，$m_n(t)$携带信号，$a$为调幅系数，$A_c(t)cos2\pi f_{c}t$为载波信号。
常规AM调制信号时域已调信号如下图所示

图2.2-2 常规调幅AM信号时域表达

为了不产生过调制现象，需要满足$|m(t)|\le 1$，即$\overline{|m(t){|}^2}\le 1$，经验得到当调制系数${\eta }_{AM}\le 0.5$，通过傅里叶变换得到已调信号的频域模型$S_{AM}(f)$为
$$\begin{gathered} S_{AM}(f)=\frac{A_c}{2}[\delta (f-f_c)+M(f-f_c)] \\ +\frac{A_c}{2}[\delta (f+f_c)+M(f+f_c)]\text{\hspace{1em}(2-7)} \end{gathered}$$
$m(t)$与AM信号的频谱如下图所示：

图2.2-3 AM信号频域图

与DSB-SC相同，AM调制得到的信号频谱属于双边带信号，因此同样可以通过一个带通滤波器BPF滤除噪声，以达到提高信噪比的目的，其时域信道模型如下图所示：

图2.2-4 AM已调信号通过BPF滤除部分噪声

常规AM调制的最大优点是接收机简单，可以采用包络检波法。因此，其信道解调模型如下：

图2.2-5 AM信号的接收与解调

2.3 SSC

单边带调制SSB的调制模型，主要是通过➀滤波法➁ 相移法(希尔伯特变换)，得到已调信号的上边带或下边带信号，从而达到减小已调信号带宽、进一步抑制噪声功率的目的。其时域信道调制模型为：

图2.3-1 滤波法

图2.3-2 相移法

以下主要就第一种方法——滤波法，进行理论分析：

由上图单边带信道的调制模型，可以得到单边带调制的公式为
$$s_{SSB}(t)=A_{c}m(t)cos2\pi f_{c}t\ast h_{SSB}(t)\text{\hspace{1em}(2-8)}$$
其中，$h_{SSB}(t)$为上边带或下边带滤波器，进行傅里叶变换，可以得到
$$S_{SSB}(f)=\frac{A_c}{2}[M(f-f_c)+M(f+f_c)]H_{SSB}(f)\text{\hspace{1em}(2-9)}$$
因此得到频域的滤波原理图为：

图2.3-3 滤波法得到单边带原理

其中，(a)为下边带滤波，(b)为上边带滤波

3、抗噪声性能理论分析

3.1 抗噪性能通用模型

图3.1-1 抗噪通用模型

在上图中，$m(t)$是带宽为$W$，功率为$P_m$的模拟基带信号，$s(t)$是带宽为$B$，功率为$P_R$的已调信号，$n(t)$是功率密度谱为$\frac{N_0}{2}$，带宽为$B_T$的加性高斯白噪声，理想带通滤波器的带宽为B，滤波器的输出信号为：
$$s_n(t)=s(t)+n(t)\text{\hspace{1em}(3-1)}$$
抗噪理论分析，主要针对已通过BPF的输入信噪比$SN{R}_{in}$、接收机输出信号的输出信噪比$SN{R}_{out}$以及解调增益和系统增益进行衡量。
其中，我们将信噪比总体增益定义为：
$$G=\frac{SN{R}_{out}}{SN{R}_{in}}\text{\hspace{1em}(3-2)}$$
输入信噪比与输出信噪比分别定义为
$$SN{R}_{in}=\frac{P_M}{N_o{B}_T}\text{\hspace{1em}(3-3)}$$
$$SN{R}_{out}=\frac{k^2\overline{m^2(t)}}{P_n}\text{\hspace{1em}(3-4)}$$
由于载波函数$c(t)=A_{c}cos(2\pi f_{c}t)$可能具有不同的幅值$A_c$，因此式中定义$k$反映不同调制方式对消息信号的影响。

3.2 DSB-SC抗噪性能

图3.2-1 DSB-SC模型

对于已调的DSB信号$s(t)=A_{c}m(t)cos(2\pi f_{c}t)$，其信号的平均功率为
$$P_{DSB}=\frac{1}{2}{A}_c^2\overline{m^2(t)}\text{\hspace{1em}(3-5)}$$
其中，可另$P_m=\overline{m^2(t)}$，由此可以计算输入信噪比
$$SN{R}_{in}=\frac{\frac{1}{2}{A}_c^2{P}_m}{2N_{o}B}\text{\hspace{1em}(3-6)}$$
其中可以得到输入信号功率与输入噪声功率分别为
$$P_{in}=\frac{1}{2}{A}_c^2{P}_m\text{\hspace{1em}(3-7)}$$
$$N_{in}=2N_{o}B\text{\hspace{1em}(3-8)}$$
因此根据接收机的特点可以得到输出信号功率与输出噪声功率分别为
$$P_{out}=\frac{1}{2}{P}_{in}\text{\hspace{1em}(3-9)}$$
$$N_{out}=\frac{1}{4}{N}_{in}\text{\hspace{1em}(3-10)}$$
因此输出信噪比与输入信噪比的关系为
$$SN{R}_{out}=2\cdot SN{R}_{in}\text{\hspace{1em}(3-11)}$$

3.3 常规AM抗噪性能

图3.3-1 常规调幅AM抗噪模型

对于已调的AM信号$s_{AM}(t)=A_c[1+a\cdot m_n(t)]cos(2\pi f_{c}t)$，其信号的平均功率为
$$\begin{gathered} P_{AM}=\overline{A_c^2[1+m^2(t){]}^{2}co{s}^{2}2\pi f_{c}t} \\ =\frac{1}{2}{A}_c^2+\frac{1}{2}{A}_c^2\overline{m^2(t)}+A_c^2\overline{m(t)}\text{\hspace{1em}(3-12)} \end{gathered}$$
由于$m(t)$为纯交流信号，则$\overline{m(t)}=0$，因此AM信号的功率为
$$P_{AM}=P_c+P_m=\frac{1}{2}{A}_c^2+\frac{1}{2}{A}_c^2\overline{m^2(t)}\text{\hspace{1em}(3-13)}$$
由于AM与DSB在频域具有同样的带宽，因此对于AM的输入噪声功率，可以将其等同于DSB的输入噪声功率
$$S_{in}=2N_{o}B\text{\hspace{1em}(3-14)}$$
由此得到输入信噪比为
$$SN{R}_{in}=\frac{P_{AM}}{S_{in}}\text{\hspace{1em}(3-15)}$$
因此，若是以相干解调的方式对AM信号进行解调，可以得到AM调制的输出信噪比为
$$SN{R}_{out}=\frac{\frac{1}{2}{P}_{AM}}{\frac{1}{4}{S}_{in}}=2\cdot SN{R}_{in}\text{\hspace{1em}(3-16)}$$

3.4 SSC抗噪性能

图3.1-1 单边带SSC抗噪模型

对于单边带调制，其已调信号为
$$s_{SSB}(t)=A_{c}m(t)cos2\pi f_{c}t\pm A_c\hat{m}(t)sin2\pi f_{c}t\text{\hspace{1em}(3-17)}$$
不同于DSB与AM调制，SSC的频域带宽只有DSB-SC的一半，因此通过带通滤波器后，噪声功率谱密度的带宽也为DSB-SC的一半，可以得到
$$SN{R}_{in}=\frac{A_c^2\overline{m^2(t)}}{4N_{o}B}\text{\hspace{1em}(3-18)}$$
通过BPF、解调与低通滤波后，得到输出信号的信号功率与噪声功率为
$$P_{out}=\frac{P_{in}}{4}=\frac{A_c^2\overline{m^2(t)}}{4}\text{\hspace{1em}(3-19)}$$
$$S_{out}=\frac{S_{in}}{4}=N_{o}B\text{\hspace{1em}(3-20)}$$
因此
$$SN{R}_{out}=\frac{\frac{A_c^2\overline{m^2(t)}}{4}}{\frac{4N_{o}B}{4}}=SN{R}_{in}\text{\hspace{1em}(3-21)}$$

4、仿真实现与仿真结果

4.1 DSB-SC

4.1.1 仿真模型

图4.1-1 DSB-SC信道模型

4.1.2 MATLAB实现

%----------DSB调制与噪声抑制-------------
%@Author:Fangy
%@Date:2020.11.6
%--------------------------------------

clear;close all;

%--------------序列定义---------------
T_start = 0;                %开始时间
T_stop = 1;                 %截止时间
fm = 10;                    %消息信号频率
fc = 100;                   %载波频率
T = T_stop - T_start;       %仿真持续时间
T_sample = 1/1000;          %采样间隔
f_sample = 1/T_sample;      %采样速率
N_sample = T/T_sample;      %采样点数
n = 0:N_sample;             %矩阵序列
%-------------------------------------

%调制
m_t = cos(2*pi*fm*n*T_sample);          %定义调制信号
s_t = m_t.*cos(2*pi*fc*n*T_sample);     %用载波对调制信号进行DSB双边带调制

N0 = 4e-6;                              %噪声功率谱密度的强度
n_t = wgn(1,length(n),N0/2);            %产生加性高斯白噪声

SNR = 10;                                %信噪比SNR单位dB
sn_t = awgn(s_t,SNR,'measured');         %向s_t中加噪

%计算输入信号功率、噪声功率、信噪比(未经过BPF)
Pin(1,1) = sum(abs(s_t).^2)/length(s_t);         %输入信号功率
Sin(1,1) = sum(abs(sn_t-s_t).^2)/length(sn_t);   %噪声功率
SNR_in(1,1) = 10*log10(Pin(1,1)/Sin(1,1));       %信噪比

%-------------绘图1：时域信号--------------
figure;
subplot(311);plot(n*T_sample,m_t);
xlabel("Time t");ylabel("Ampl |m(t)|");title("m(t)");
subplot(312);plot(n*T_sample,s_t);
xlabel("Time t");ylabel("Ampl |s(t)|");title("s(t)");
subplot(313);plot(n*T_sample,sn_t);
xlabel("Time t");ylabel("Ampl |s(t)+n(t)|");title("s(t)叠加高斯白噪声");
%------------------------------------------

%傅里叶变换 频域频谱表达
f_res = f_sample/N_sample;  %频率分辨率
f_max = f_res*N_sample/2;   %最大频率

m_f = abs(fft(m_t))/f_sample;                                         %m(t)傅里叶变换得到M(f)幅度谱并归一化
% m_f_rearrange = [m_f(N_sample/2+1:N_sample-1),m_f(1:N_sample/2)];   %前后调转幅度谱顺序(原理实现)
m_f_rearrange = fftshift(m_f);                                        %前后调转幅度谱顺序(简单调用)

s_f = abs(fft(s_t))/f_sample;       %s(t)傅里叶变换得到S(f)幅度谱并归一化
s_f_rearrange = fftshift(s_f);      %前后调转幅度谱顺序(简单调用)

sn_f = abs(fft(sn_t))/f_sample;     %sn(t)傅里叶变换得到SN(f)幅度谱并归一化
sn_f_rearrange = fftshift(sn_f);    %前后调转幅度谱顺序(简单调用)

%-------------绘图2：频域信号--------------
figure;
subplot(311);plot((-N_sample/2+1:N_sample/2-1)*f_res,m_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |M(f)|");title("M(f)");
subplot(312);plot((-N_sample/2+1:N_sample/2-1)*f_res,s_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |S(f)|");title("S(f)");
subplot(313);plot((-N_sample/2+1:N_sample/2-1)*f_res,sn_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |S(f)+N(f)|");title("S(f)加噪频谱");
%------------------------------------------

%----------------带通滤波器-----------------
%通过filterDesigner设计：
%ResponseType选Bandpass
%fs=1000;Fstop1=84;Fpass1=90;Fpass2=110;Fstop2=116;
%Design Method选择FIR->Window
%点击DesignFilter后导出数据命名为Num_BPF
load Num_BPF_DSB;
%-------------------------------------------

%有噪信号通过BPF
sn_BPF_t = conv(sn_t,Num_BPF);              %时域卷积相当于频域相乘
%---------------抵消fdatool产生的相移---------------
count = 0;
sn_BPF_t_rearrange = zeros(1,1001);
for i=length(sn_BPF_t)/2-500:1:length(sn_BPF_t)/2+500
    count = count + 1;
    sn_BPF_t_rearrange(1,count) = sn_BPF_t(1,i);
end
%--------------------------------------------------

n_BPF_t = s_t - sn_BPF_t_rearrange;         %通过BPF后剩下的噪声

%计算输入信号功率、噪声功率、信噪比(已通过BPF)
Pin(2,1) = sum(abs(s_t).^2)/length(s_t);                %输入信号功率
Sin(2,1) = sum(abs(n_BPF_t).^2)/length(n_BPF_t);        %噪声功率
SNR_in(2,1) = 10*log10(Pin(2,1)/Sin(2,1));              %信噪比

%-------------绘图3：通过带通滤波器--------------
figure;
subplot(311);plot(n*T_sample,sn_BPF_t_rearrange);
xlabel("Time t");ylabel("Ampl |sn_{BPF}(t)|");title("sn(t)通过BPF");
subplot(312);plot(n*T_sample,s_t);
xlabel("Time t");ylabel("Ampl |s(t)|");title("无噪s(t)");
subplot(313);plot(n*T_sample,n_BPF_t);
xlabel("Time t");ylabel("Ampl |n_{BPF}(t)|");title("通过BPF后剩余的噪声");
%---------------------------------------------

%解调——时域
sd_t = s_t .* cos(2*pi*fc*n*T_sample);                      %无噪解调
sd_n_t = sn_BPF_t_rearrange .* cos(2*pi*fc*n*T_sample);     %有噪解调

%解调——频域
sd_f = abs(fft(sd_t))/f_sample;                             %无噪解调
sd_f_rearrange = fftshift(sd_f);                            %无噪解调
sd_n_f = abs(fft(sd_n_t))/f_sample;                         %有噪解调
sd_n_f_rearrange = fftshift(sd_n_f);                        %有噪解调

%通过低通滤波器（无噪）
%----------------低通滤波器-----------------
%通过filterDesigner设计：
%ResponseType选择Lowpass
%Fs=1000;Fpass=10;Fstop=16;
%Design Method选择FIR->Window
%点击DesignFilter后导出数据命名为Num_LPF
load Num_LPF_DSB
%-------------------------------------------
mo_t = conv(sd_t,Num_LPF);
%---------------抵消fdatool产生的相移---------------
count = 0;
mo_t_rearrange = zeros(1,1001);
for i=length(mo_t)/2-500:1:length(mo_t)/2+500
    count = count + 1;
    mo_t_rearrange(1,count) = mo_t(1,i);
end
%--------------------------------------------------

%通过低通滤波器（有噪）
mo_n_t = conv(sd_n_t,Num_LPF);
%---------------抵消fdatool产生的相移---------------
count = 0;
mo_n_t_rearrange = zeros(1,1001);
for i=length(mo_n_t)/2-500:1:length(mo_n_t)/2+500
    count = count + 1;
    mo_n_t_rearrange(1,count) = mo_n_t(1,i);
end
%--------------------------------------------------

%-------------绘图4：解调-时域--------------
figure;
subplot(411);plot(n*T_sample,sd_t);
xlabel("Time t");ylabel("Ampl |s_{d}(t)|");title("解调s_{d}(t)未通过LPF(无噪)");
subplot(412);plot(n*T_sample,sd_n_t);
xlabel("Time t");ylabel("Ampl |sn_{d}(t)|");title("解调s_{d}(t)未通过LPF(有噪)");
subplot(413);plot(n*T_sample,mo_t_rearrange);
xlabel("Time t");ylabel("Ampl |m_{o}(t)|");title("通过LPF后输出信号(无噪)");
subplot(414);plot(n*T_sample,mo_n_t_rearrange);
xlabel("Time t");ylabel("Ampl |m_{no}(t)|");title("通过LPF后输出信号(有噪)");
%---------------------------------------------

%傅里叶变换频域观察输出频谱
mo_f = abs(fft(mo_t_rearrange))/f_sample;       %频域解调(无噪)
mo_f_rearrange = fftshift(mo_f);                %频域解调(无噪)
mo_n_f = abs(fft(mo_n_t_rearrange))/f_sample;   %频域解调(有噪)
mo_n_f_rearrange = fftshift(mo_n_f);            %频域解调(有噪)

%-------------绘图5：解调-频域--------------
figure;
subplot(411);plot((-N_sample/2+1:N_sample/2-1)*f_res,sd_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |s_{d}(f)|");title("解调s_{f}(t)未通过LPF(无噪)");
subplot(412);plot((-N_sample/2+1:N_sample/2-1)*f_res,sd_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |sn_{d}(f)|");title("解调s_{d}(t)未通过LPF(有噪)");
subplot(413);plot((-N_sample/2+1:N_sample/2-1)*f_res,mo_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |m_{o}(f)|");title("通过LPF输出m_{o}(f)(无噪)");
subplot(414);plot((-N_sample/2+1:N_sample/2-1)*f_res,mo_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |m_{no}(f)|");title("通过LPF输出m_{no}(f)(有噪)");
%---------------------------------------------

no_t = mo_t_rearrange-mo_n_t_rearrange;
%计算输出信号功率、噪声功率、信噪比(已通过LPF)
Pin(3,1) = sum(abs(mo_t_rearrange).^2)/length(mo_t_rearrange);          %输出信号功率
Sin(3,1) = sum(abs(no_t).^2)/length(no_t);                              %噪声功率
SNR_in(3,1) = 10*log10(Pin(3,1)/Sin(3,1));                              %信噪比

disp("输入信噪比："+SNR_in(2,1));
disp("输出信噪比："+SNR_in(3,1));

4.1.3 仿真结果

图4.1-2 消息信号，调制信号，加入噪声信号的时域

图4.1-3 消息信号与调制后信号，加入噪声信号的傅里叶变换

图4.1-4 有噪与无噪的DSB信号通过BPF、通过BPF后剩余噪声

图4.1-5 无噪与有噪的解调信号、输出信号（时域）

图4.1-6 无噪与有噪的解调信号、输出信号(频域)

4.1.4 仿真结果与理论结果对比

图4.1-7 DSB-SC理论结果与仿真结果对比

产生误差的原因：
1、由MATLAB产生的加性高斯白噪声并不完全符合随机分布，因此双边功率谱密度不完全遵守幅值均为$N_0/2$
2、信号通过滤波器时(如BPF、LPF)，在MATLAB中用卷积抽样函数得以实现，因此产生一定相移，即时通过将相位还原，滤波后的信号与原始信号在相位上仍然存在一定的误差
3、仿真通过离散的抽样点实现模拟信号的仿真，因此仿真信号并不是真正意义上的连续信号，因此产生误差，若是提高时域的采样速率，误差会有一定程度的减小，但是对应的频域采样频率增大，并不利于频域分析，因此选择一个合适的采样频率很重要，这个采样频率不能过大或过小。

4.2 AM

4.2.1 仿真模型

图4.2-1 常规AM抗噪信道模型

4.2.2 MATLAB实现

%----------AM调制与噪声抑制-------------
%@Author:Fangy
%@Date:2020.11.7
%--------------------------------------

clear;close all;

%--------------序列定义---------------
T_start = 0;                %开始时间
T_stop = 1;                 %截止时间
fm = 10;                    %消息信号频率
fc = 100;                   %载波频率
T = T_stop - T_start;       %仿真持续时间
T_sample = 1/1000;          %采样间隔
f_sample = 1/T_sample;      %采样速率
N_sample = T/T_sample;      %采样点数
n = 0:N_sample;             %矩阵序列
%-------------------------------------

%调制
m_t = cos(2*pi*fm*n*T_sample);          %定义调制信号
a = 0.5;                                %调制系数
c_t = cos(2*pi*fc*n*T_sample);          %载波信号
Sam_t = (1+a*m_t).*c_t;                 %AM调制

SNR = 1;                                    %信噪比SNR单位dB
Sam_n_t = awgn(Sam_t,SNR,'measured');        %向s_t中加噪

%计算输入信号功率、噪声功率、信噪比(未经过BPF)
Pin(1,1) = sum(abs(Sam_t).^2)/length(Sam_t);                   %输入信号功率
Sin(1,1) = sum(abs(Sam_n_t-Sam_t).^2)/length(Sam_n_t-Sam_t);   %噪声功率
SNR_in(1,1) = 10*log10(Pin(1,1)/Sin(1,1));                     %信噪比

%-------------绘图1：时域信号--------------
figure;
subplot(311);plot(n*T_sample,m_t);
xlabel("Time t");ylabel("Ampl |m(t)|");title("m(t)");
subplot(312);plot(n*T_sample,Sam_t);
xlabel("Time t");ylabel("Ampl |s_{AM}(t)|");title("AM调制信号s_{AM}(t)");
subplot(313);plot(n*T_sample,Sam_n_t);
xlabel("Time t");ylabel("Ampl |s_{AM}(t)+n(t)|");title("s_{AM}(t)叠加高斯白噪声");
%------------------------------------------

%傅里叶变换 频域频谱表达
f_res = f_sample/N_sample;  %频率分辨率
f_max = f_res*N_sample/2;   %最大频率

m_f = abs(fft(m_t))/f_sample;                                         %m(t)傅里叶变换得到M(f)幅度谱并归一化
% m_f_rearrange = [m_f(N_sample/2+1:N_sample-1),m_f(1:N_sample/2)];   %前后调转幅度谱顺序(原理实现)
m_f_rearrange = fftshift(m_f);                                        %前后调转幅度谱顺序(简单调用)

Sam_f = abs(fft(Sam_t))/f_sample;         %s(t)傅里叶变换得到S(f)幅度谱并归一化
Sam_f_rearrange = fftshift(Sam_f);        %前后调转幅度谱顺序(简单调用)

Sam_n_f = abs(fft(Sam_n_t))/f_sample;     %sn(t)傅里叶变换得到SN(f)幅度谱并归一化
Sam_n_f_rearrange = fftshift(Sam_n_f);    %前后调转幅度谱顺序(简单调用)

%-------------绘图2：频域信号--------------
figure;
subplot(311);plot((-N_sample/2+1:N_sample/2-1)*f_res,m_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |M(f)|");title("M(f)");
subplot(312);plot((-N_sample/2+1:N_sample/2-1)*f_res,Sam_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |S_{AM}(f)|");title("无噪S_{AM}(f)");
subplot(313);plot((-N_sample/2+1:N_sample/2-1)*f_res,Sam_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |SN_{AM}(f)|");title("有噪SN_{AM}(f)");
%-----------------------------------------

%----------------带通滤波器-----------------
%通过filterDesigner设计：
%ResponseType选Bandpass
%fs=1000;Fstop1=84;Fpass1=90;Fpass2=110;Fstop2=116;
%Design Method选择FIR->Window
%点击DesignFilter后导出数据命名为Num_BPF
load Num_BPF_DSB;
%-------------------------------------------

%有噪信号通过BPF
Sam_n_BPF_t = conv(Sam_n_t,Num_BPF);              %时域卷积相当于频域相乘
%---------------抵消fdatool产生的相移---------------
count = 0;
Sam_n_BPF_t_rearrange = zeros(1,1001);
for i=length(Sam_n_BPF_t)/2-500:1:length(Sam_n_BPF_t)/2+500
    count = count + 1;
    Sam_n_BPF_t_rearrange(1,count) = Sam_n_BPF_t(1,i);
end
%--------------------------------------------------

n_BPF_t = Sam_t - Sam_n_BPF_t_rearrange;         %通过BPF后剩下的噪声

%计算输入信号功率、噪声功率、信噪比(已通过BPF)
Pin(2,1) = sum(abs(Sam_t).^2)/length(Sam_t);            %输入信号功率
Sin(2,1) = sum(abs(n_BPF_t).^2)/length(n_BPF_t);        %噪声功率
SNR_in(2,1) = 10*log10(Pin(2,1)/Sin(2,1));              %信噪比

%-------------绘图3：通过带通滤波器--------------
figure;
subplot(311);plot(n*T_sample,Sam_n_BPF_t_rearrange);
xlabel("Time t");ylabel("Ampl |Sam_{n BPF}(t)|");title("有噪Sam_{n BPF}(t)通过BPF");
subplot(312);plot(n*T_sample,Sam_t);
xlabel("Time t");ylabel("Ampl |Sam(t)|");title("无噪Sam(t)");
subplot(313);plot(n*T_sample,n_BPF_t);
xlabel("Time t");ylabel("Ampl |n_{BPF}(t)|");title("通过BPF后剩余的噪声");
%---------------------------------------------

%解调——时域
Sam_d_t = Sam_t .* cos(2*pi*fc*n*T_sample);                      %无噪解调
Sam_d_n_t = Sam_n_BPF_t_rearrange .* cos(2*pi*fc*n*T_sample);    %有噪解调

%解调——频域
Sam_d_f = abs(fft(Sam_d_t))/f_sample;                             %无噪解调
Sam_d_f_rearrange = fftshift(Sam_d_f);                            %无噪解调
Sam_d_n_f = abs(fft(Sam_d_n_t))/f_sample;                         %有噪解调
Sam_d_n_f_rearrange = fftshift(Sam_d_n_f);                        %有噪解调

%通过低通滤波器（无噪）
%----------------低通滤波器-----------------
%通过filterDesigner设计：
%ResponseType选择Lowpass
%Fs=1000;Fpass=10;Fstop=16;
%Design Method选择FIR->Window
%点击DesignFilter后导出数据命名为Num_LPF
load Num_LPF_DSB
%-------------------------------------------
mo_t = conv(Sam_d_t,Num_LPF);
%---------------抵消fdatool产生的相移---------------
count = 0;
mo_t_rearrange = zeros(1,1001);
for i=length(mo_t)/2-500:1:length(mo_t)/2+500
    count = count + 1;
    mo_t_rearrange(1,count) = mo_t(1,i);
end
%--------------------------------------------------

%通过低通滤波器（有噪）
mo_n_t = conv(Sam_d_n_t,Num_LPF);
%---------------抵消fdatool产生的相移---------------
count = 0;
mo_n_t_rearrange = zeros(1,1001);
for i=length(mo_n_t)/2-500:1:length(mo_n_t)/2+500
    count = count + 1;
    mo_n_t_rearrange(1,count) = mo_n_t(1,i);
end
%--------------------------------------------------

%-------------绘图4：解调-时域--------------
figure;
subplot(411);plot(n*T_sample,Sam_d_t);
xlabel("Time t");ylabel("Ampl |Sam_{d}(t)|");title("解调Sam_{d}(t)未通过LPF(无噪)");
subplot(412);plot(n*T_sample,Sam_d_n_t);
xlabel("Time t");ylabel("Ampl |Sam_{n d}(t)|");title("解调Sam_{n d}(t)未通过LPF(有噪)");
subplot(413);plot(n*T_sample,mo_t_rearrange);
xlabel("Time t");ylabel("Ampl |m_{o}(t)|");title("通过LPF后输出信号(无噪)");
subplot(414);plot(n*T_sample,mo_n_t_rearrange);
xlabel("Time t");ylabel("Ampl |m_{n o}(t)|");title("通过LPF后输出信号(有噪)");
%---------------------------------------------

%傅里叶变换频域观察输出频谱
mo_f = abs(fft(mo_t_rearrange))/f_sample;       %频域解调(无噪)
mo_f_rearrange = fftshift(mo_f);                %频域解调(无噪)
mo_n_f = abs(fft(mo_n_t_rearrange))/f_sample;   %频域解调(有噪)
mo_n_f_rearrange = fftshift(mo_n_f);            %频域解调(有噪)

%-------------绘图5：解调-频域--------------
figure;
subplot(411);plot((-N_sample/2+1:N_sample/2-1)*f_res,Sam_d_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |Sam_{d}(f)|");title("解调Sam_{d}(f)未通过LPF(无噪)");
subplot(412);plot((-N_sample/2+1:N_sample/2-1)*f_res,Sam_d_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |Sam_{n d}(f)|");title("解调Sam_{n d}(f)未通过LPF(有噪)");
subplot(413);plot((-N_sample/2+1:N_sample/2-1)*f_res,mo_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |m_{o}(f)|");title("通过LPF输出m_{o}(f)(无噪)");
subplot(414);plot((-N_sample/2+1:N_sample/2-1)*f_res,mo_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |m_{no}(f)|");title("通过LPF输出m_{no}(f)(有噪)");
%---------------------------------------------

no_t = mo_t_rearrange-mo_n_t_rearrange;
%计算输出信号功率、噪声功率、信噪比(已通过LPF)
Pin(3,1) = sum(abs(mo_t_rearrange).^2)/length(mo_t_rearrange);          %输出信号功率
Sin(3,1) = sum(abs(no_t).^2)/length(no_t);                              %噪声功率
SNR_in(3,1) = 10*log10(Pin(3,1)/Sin(3,1));                              %信噪比

disp("输入信噪比："+SNR_in(2,1));
disp("输出信噪比："+SNR_in(3,1));

4.2.3 仿真结果

图4.2-2 调制信号、AM调制信号、AM调制加噪信号（时域）

图4.2-3 调制信号、AM调制信号、AM调制加噪信号（频域）

图4.2-4 通过BPF的有噪AM信号、无噪AM信号、剩余噪声（时域）

图4.2-5 解调后有噪与无噪AM信号、通过LPF后有噪和无噪输出信号

图4.2-5 解调后有噪与无噪AM信号、通过LPF后有噪和无噪输出信号（频域）

4.2.4 仿真结果与理论结果对比

图4.2-6 常规AM理论结果与仿真结果对比

产生误差的原因：
1、由MATLAB产生的加性高斯白噪声并不完全符合随机分布，因此双边功率谱密度不完全遵守幅值均为$N_0/2$
2、信号通过滤波器时(如BPF、LPF)，在MATLAB中用卷积抽样函数得以实现，因此产生一定相移，即时通过将相位还原，滤波后的信号与原始信号在相位上仍然存在一定的误差
3、仿真通过离散的抽样点实现模拟信号的仿真，因此仿真信号并不是真正意义上的连续信号，因此产生误差，若是提高时域的采样速率，误差会有一定程度的减小，但是对应的频域采样频率增大，并不利于频域分析，因此选择一个合适的采样频率很重要，这个采样频率不能过大或过小。

4.3 SSC

4.3.1 仿真模型

图4.3-1 单边带SSC信道模型

4.3.2 MATLAB实现

%----------SSC调制与噪声抑制-------------
%@Author:Fangy
%@Date:2020.11.8
%--------------------------------------

clear;close all;

%--------------序列定义---------------
T_start = 0;                %开始时间
T_stop = 1;                 %截止时间
fm = 10;                    %消息信号频率
fc = 100;                   %载波频率
T = T_stop - T_start;       %仿真持续时间
T_sample = 1/1000;          %采样间隔
f_sample = 1/T_sample;      %采样速率
N_sample = T/T_sample;      %采样点数
n = 0:N_sample;             %矩阵序列
%-------------------------------------

%调制
m_t = cos(2*pi*fm*n*T_sample);          %定义调制信号
s_t = m_t.*cos(2*pi*fc*n*T_sample);     %用载波对调制信号进行DSB双边带调制

SNR = 5;                               %加入噪声的信噪比
sn_t = awgn(s_t,SNR);                  %向已调信号加入高斯白噪声

%--------------希尔伯特滤波器--------------
%说明：此希尔伯特滤波器并不是严格定义的希尔伯特滤波器
%      而是用带宽为fm的带通滤波器做了同等效果处理
%通过filterDesigner设计：
%ResponseType选Bandpass
%fs=1000;Fstop1=84;Fpass1=90;Fpass2=100;Fstop2=106;
%Design Method选择FIR->Window
%点击DesignFilter后导出数据命名为Num_BPF_SSC
load Num_BPF_SSC                        
%-----------------------------------------

ss_t = conv(Num_BPF_SSC,s_t);           %获得下边带信号（无噪）
ss_n_t = conv(Num_BPF_SSC,sn_t);        %获得下边带信号（有噪）
%---------------抵消fdatool产生的相移---------------
count = 0;
ss_n_t_rearrange = zeros(1,1001);
for i=length(ss_n_t)/2-500:1:length(ss_n_t)/2+500
    count = count + 1;
    ss_n_t_rearrange(1,count) = ss_n_t(1,i);
end
%--------------------------------------------------
ss_n_t = ss_n_t_rearrange;
%---------------抵消fdatool产生的相移---------------
count = 0;
ss_t_rearrange = zeros(1,1001);
for i=length(ss_t)/2-500:1:length(ss_t)/2+500
    count = count + 1;
    ss_t_rearrange(1,count) = ss_t(1,i);
end
%--------------------------------------------------
ss_t = ss_t_rearrange;

%计算输入信号功率、噪声功率、信噪比(未经过BPF)
Pin(1,1) = sum(abs(ss_t).^2)/length(ss_t);              %输入信号功率
Sin(1,1) = sum(abs(ss_n_t-ss_t).^2)/length(ss_n_t);     %噪声功率
SNR_in(1,1) = 10*log10(Pin(1,1)/Sin(1,1));              %信噪比

%-------------绘图1：时域信号--------------
figure;
subplot(311);plot(n*T_sample,m_t);
xlabel("Time t");ylabel("Ampl |m(t)|");title("m(t)");
subplot(312);plot(n*T_sample,ss_t);
xlabel("Time t");ylabel("Ampl |ss(t)|");title("单边带(下边带)已调信号ss(t)");
subplot(313);plot(n*T_sample,ss_n_t);
xlabel("Time t");ylabel("Ampl |s(t)+n(t)|");title("ss(t)叠加高斯白噪声");
%------------------------------------------

%傅里叶变换 频域频谱表达
f_res = f_sample/N_sample;  %频率分辨率
f_max = f_res*N_sample/2;   %最大频率

m_f = abs(fft(m_t))/f_sample;                                         %m(t)傅里叶变换得到M(f)幅度谱并归一化
% m_f_rearrange = [m_f(N_sample/2+1:N_sample-1),m_f(1:N_sample/2)];   %前后调转幅度谱顺序(原理实现)
m_f_rearrange = fftshift(m_f);                                        %前后调转幅度谱顺序(简单调用)

ss_f = abs(fft(ss_t))/f_sample;       %s(t)傅里叶变换得到S(f)幅度谱并归一化
ss_f_rearrange = fftshift(ss_f);      %前后调转幅度谱顺序(简单调用)

ss_n_f = abs(fft(ss_n_t))/f_sample;     %sn(t)傅里叶变换得到SN(f)幅度谱并归一化
ss_n_f_rearrange = fftshift(ss_n_f);    %前后调转幅度谱顺序(简单调用)

%-------------绘图2：频域信号--------------
figure;
subplot(311);plot((-N_sample/2+1:N_sample/2-1)*f_res,m_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |M(f)|");title("M(f)");
subplot(312);plot((-N_sample/2+1:N_sample/2-1)*f_res,ss_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |S(f)|");title("S(f)");
subplot(313);plot((-N_sample/2+1:N_sample/2-1)*f_res,ss_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |S(f)+N(f)|");title("S(f)加噪频谱");
%------------------------------------------

%解调——时域
ssd_t = ss_t .* cos(2*pi*fc*n*T_sample);                    %无噪解调
ssd_n_t = ss_n_t .* cos(2*pi*fc*n*T_sample);                %有噪解调

%解调——频域
ssd_f = abs(fft(ssd_t))/f_sample;                            %无噪解调
ssd_f_rearrange = fftshift(ssd_f);                           %无噪解调
ssd_n_f = abs(fft(ssd_n_t))/f_sample;                        %有噪解调
ssd_n_f_rearrange = fftshift(ssd_n_f);                       %有噪解调

%通过低通滤波器（无噪）
%----------------低通滤波器-----------------
%通过filterDesigner设计：
%ResponseType选择Lowpass
%Fs=1000;Fpass=10;Fstop=16;
%Design Method选择FIR->Window
%点击DesignFilter后导出数据命名为Num_LPF
load Num_LPF_DSB
%-------------------------------------------
mo_t = conv(ssd_t,Num_LPF);
%---------------抵消fdatool产生的相移---------------
count = 0;
mo_t_rearrange = zeros(1,1001);
for i=length(mo_t)/2-500:1:length(mo_t)/2+500
    count = count + 1;
    mo_t_rearrange(1,count) = mo_t(1,i);
end
%--------------------------------------------------

%通过低通滤波器（有噪）
mo_n_t = conv(ssd_n_t,Num_LPF);
%---------------抵消fdatool产生的相移---------------
count = 0;
mo_n_t_rearrange = zeros(1,1001);
for i=length(mo_n_t)/2-500:1:length(mo_n_t)/2+500
    count = count + 1;
    mo_n_t_rearrange(1,count) = mo_n_t(1,i);
end
%--------------------------------------------------

%-------------绘图3：解调-时域--------------
figure;
subplot(411);plot(n*T_sample,ssd_t);
xlabel("Time t");ylabel("Ampl |ss_{d}(t)|");title("解调ss_{d}(t)未通过LPF(无噪)");
subplot(412);plot(n*T_sample,ssd_n_t);
xlabel("Time t");ylabel("Ampl |ssn_{d}(t)|");title("解调ss_{d}(t)未通过LPF(有噪)");
subplot(413);plot(n*T_sample,mo_t_rearrange);
xlabel("Time t");ylabel("Ampl |m_{o}(t)|");title("通过LPF后输出信号(无噪)");
subplot(414);plot(n*T_sample,mo_n_t_rearrange);
xlabel("Time t");ylabel("Ampl |m_{no}(t)|");title("通过LPF后输出信号(有噪)");
%---------------------------------------------

%傅里叶变换频域观察输出频谱
mo_f = abs(fft(mo_t_rearrange))/f_sample;       %频域解调(无噪)
mo_f_rearrange = fftshift(mo_f);                %频域解调(无噪)
mo_n_f = abs(fft(mo_n_t_rearrange))/f_sample;   %频域解调(有噪)
mo_n_f_rearrange = fftshift(mo_n_f);            %频域解调(有噪)

%-------------绘图4：解调-频域--------------
figure;
subplot(411);plot((-N_sample/2+1:N_sample/2-1)*f_res,ssd_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |ss_{d}(f)|");title("解调ss_{f}(t)未通过LPF(无噪)");
subplot(412);plot((-N_sample/2+1:N_sample/2-1)*f_res,ssd_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |ssn_{d}(f)|");title("解调ss_{d}(t)未通过LPF(有噪)");
subplot(413);plot((-N_sample/2+1:N_sample/2-1)*f_res,mo_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |m_{o}(f)|");title("通过LPF输出m_{o}(f)(无噪)");
subplot(414);plot((-N_sample/2+1:N_sample/2-1)*f_res,mo_n_f_rearrange(1:N_sample-1));
xlabel("frequency f");ylabel("Ampl |m_{no}(f)|");title("通过LPF输出m_{no}(f)(有噪)");
%---------------------------------------------

no_t = mo_t_rearrange-mo_n_t_rearrange;
%计算输出信号功率、噪声功率、信噪比(已通过LPF)
Pin(3,1) = sum(abs(mo_t_rearrange).^2)/length(mo_t_rearrange);          %输出信号功率
Sin(3,1) = sum(abs(no_t).^2)/length(no_t);                              %噪声功率
SNR_in(3,1) = 10*log10(Pin(3,1)/Sin(3,1));                              %信噪比

disp("输入信噪比："+SNR_in(1,1));
disp("输出信噪比："+SNR_in(3,1));

4.3.3 仿真结果

图4.3-2 调制信号、SSC调制信号、SSC调制加噪信号（时域）

图4.3-3 调制信号、SSC调制信号、SSC调制加噪信号（频域）

图4.3-4 解调后有噪与无噪SSC信号、通过LPF后有噪和无噪输出信号

图4.3-5 解调后有噪与无噪SSC信号、通过LPF后有噪和无噪输出信号(频域)

4.3.4 仿真结果与理论结果对比

图4.3-6 SSC理论结果与仿真结果对比

产生误差的原因：
1、由MATLAB产生的加性高斯白噪声并不完全符合随机分布，因此双边功率谱密度不完全遵守幅值均为$N_0/2$
2、信号通过滤波器时(如BPF、LPF)，在MATLAB中用卷积抽样函数得以实现，因此产生一定相移，即时通过将相位还原，滤波后的信号与原始信号在相位上仍然存在一定的误差
3、仿真通过离散的抽样点实现模拟信号的仿真，因此仿真信号并不是真正意义上的连续信号，因此产生误差，若是提高时域的采样速率，误差会有一定程度的减小，但是对应的频域采样频率增大，并不利于频域分析，因此选择一个合适的采样频率很重要，这个采样频率不能过大或过小。

5、小结

本次仿真实验对信号的调制解调、噪声的抑制以及滤波器的设计三个方面给出理论设计以及仿真设计。
在对滤波器的设计中，采用MATLAB内置filterDesigner功能对带通滤波器BPF与低通滤波器LPF做出设计，由于设计得到的滤波器的采样点数与原始信号不一致，因此原始时域信号与滤波器函数进行卷积后的信号长度并不是完全遵守二倍信号函数长度，因此需要对滤波后的信后做出相位调整，从而在后续比对剩余噪声时可以极大减少误差。
在噪声的抑制方面，采用高斯加性白噪声，由于MATLAB所产生的随机噪声信号的随机性并不是完全符合随机分布，因此在进行噪声比对时会产生一定误差。

6、参考文献

[1] 现代通信原理6.1 常规调幅调制(AM)与抑制载波双边带(DSB-SC)调制
[2] 仿真作业3：噪声通过DSB-SC解调器
[3] 现代通信原理A.2：FIR低通滤波器设计
[4] 现代通信原理6.2：单边带(SSB)调制
[5] 现代通信原理仿真报告总结
[6] 模拟幅度调制相干解调系统抗噪声性能仿真分析