数字通信系统误码性能仿真分析

数字通信系统误码性能仿真分析

一、概述

从《通信原理》课本中我们知道，系统的误码率定义为错误的码元个数除以传输的总码元个数。显然可以通过对于设计出来的系统，进行大量的试验测试，统计结果来求出其误码率，但是这不仅不利于我们设计新的通信系统，而且会花费大量时间和精力。因此我们希望能够通过理论分析的方法来获得误码率的结果（差错概率）。这个结果实际上是码元错误的概率，而上述测试得到的结果会无限接近这个概率。因而，对各类系统的误码性能进行理论分析，是非常重要且活跃的研究领域。对于数字传输系统，要想分析其误码性能，首先要了解对于不同的数字传输系统的模型结构才能进行分析。

数字传输系统的模型主要是：在信息发送端，信源发送信号，通过信源编码、信道编码以及数字调制器，之后通过信道传输，在信息接收端，通过数字解调器、信道译码以及信源译码，通过信宿接收。这里信道我们主要研究高斯白噪声信道（$AWGN$)，从带宽来看，$AWGN$信道主要分为：带宽无限的$AWGN$信道和带宽有限的$AWGN$信道，而其中的接收滤波器一般都用匹配滤波器，从通信信道的传递函数来看，$AWGN$信道主要分为：基带信道和带通信道，其分别进行的是基带传输和频带传输，下面我们就以下对四种情况进行分析：

1. 信道带宽无限时的单极性基带传输

2. 信道带宽受限时的双极性基带传输

3. 信道带宽受限时的$QPSK$传输

4. 信道带宽受限时的$16QAM$传输

本文通过对系统模型分析，理论推导不同情况下的误码率，然后再利用$MATLAB$进行仿真，比较仿真值与理论值的区别，并进行分析。

二、系统模型

2.1 数字基带传输系统

2.1.1 数字系统基带传输模型

图2-1 数字系统基带传输模型

2.1.2 信道带宽无限时的单极性基带传输

对于单极性信号，我们用二进制编码信号进行分析，其离散幅度值为 { a n } = { A , 0 } \{a_n\}=\{A,0\} {an​}={A,0}，由于不归零码可以看成归零码的一个特例（当$\tau =T_s$时），所以归零码和不归零码可以看成是不同信息速率的区别，因此这里可以用不归零码进行分析，其发送滤波器冲激响应为：
$$g_T(t)=\{\begin{array}{l}1，0\le t\le T_s\\ 0，otherwise\end{array}\text{\hspace{1em}(2.1.1)}$$
这样发射机的发射信号为：
$$s_i(t)=\{\begin{array}{l}{s}_1(t)=A，b_i{=}^{\prime }{1}^{\prime }\\ s_2(t)=0，b_i{=}^{\prime }{0}^{\prime }\end{array}\text{\hspace{1em}(2.1.2)}$$
所谓带宽无限$AWGN$信道，是指信道模块的频率传递函数$C(f)$在整个频率轴上都是一个恒定的值，如图$2-3(a)$所示，发送信号$s_i(t)$通过$C(f)$时不会发生形状上的变化。进一步，信号在进入接收机之前，叠加加性高斯白噪声$n_W(t)$，因此可以得到接收信号为：
$$r(t)=s(t)+n_w(t)\text{\hspace{1em}(2.1.3)}$$
为了达到更佳的误码率，这里的接收滤波器采用匹配滤波器，让其与输入信号$s(t)$进行匹配，得到匹配滤波器的冲激响应为：
$$g_R(t)=Ks(T-t)\text{\hspace{1em}(2.1.4)}$$
其中$T=T_s$，令$K=1$，得到匹配滤波器的冲击响应如下式子：
$$g_R(t)=s(T_s-t)\text{\hspace{1em}(2.1.5)}$$

图2-2 带宽无限传输模型

图2-3 频带无限与频带有限信道的频率传递函数示意图<\center>

2.1.3 信道带宽受限时的双极性基带传输

对于双极性信号，如上分析仍采用二进制编码信号以及归零码进行分析，其离散幅度值为 { a n } = { A , − A } \{a_n\}=\{A,-A\} {an​}={A,−A}，其发送滤波器冲激响应为：
$$g_T(t)=\{\begin{array}{l}\pm 1，0\le t\le T_s\\ 0，otherwise\end{array}\text{\hspace{1em}(2.1.6)}$$
这样发射机的发射信号为：
$$s_i(t)=\{\begin{array}{l}{s}_1(t)=A，b_i{=}^{\prime }{1}^{\prime }\\ s_2(t)=-A，b_i{=}^{\prime }{0}^{\prime }\end{array}\text{\hspace{1em}(2.1.7)}$$
所谓带宽有限的$AWGN$信道，指的是信道模块的频率传递函数$C(f)$类似于理想低通滤波器，设其带宽为$B$，如上图$2-3(b)$。显然，对于这类信道，发送信号$s(t)$通过$C(f)$时，其波形及频谱会发生改变，得到的接收信号为：
$$r(t)=s(t)\ast c(t)+n_w(t)\text{\hspace{1em}(2.1.8)}$$
因此，这里我们将发送滤波器、信道、接受滤波器看成一个成型滤波器进行分析，则成型滤波器的冲激响应为：
$$h(t)=g_t(t)c(t)g_R(t)\text{\hspace{1em}(2.1.9)}$$

$$H(f)=G_T(f)C(f)G_R(f)\text{\hspace{1em}(2.1.10)}$$

这样由于信道频带有限，使得输出的信号频带有限，这样信号的时域就是无限长的了，就可能产生码间干扰，因此为了避免码间干扰，分析知，需要$h(t)$满足如下条件：
$$h(n{T}_s)=\{\begin{array}{l}1，n=0\\ 0，n=1\end{array}\text{\hspace{1em}(2.1.11)}$$
而由奈奎斯特第一准则得，满足式子$(2.1.11)$的充分必要条件是$h(t)$的傅里叶变换$H(f)$必须满足：
$$\sum _{m=-\infty }^{+\infty }H(f+\frac{m}{T_s})=T_s\text{\hspace{1em}(2.1.12)}$$
经过分析知，要满足这个条件，$h(t)$需要为升余弦信号，即：
$$h(t)=sa(\frac{\pi t}{T_s})\frac{cos(\frac{\alpha \pi t}{T_s})}{1-\frac{4a^2{t}^2}{T_s}}\text{\hspace{1em}(2.1.13)}$$

H ( f ) = { T s ， 0 ≤ ∣ f ∣ ≤ 1 − α 2 T s T s 2 { 1 + c o s [ π T s α ( ∣ f ∣ − 1 − α 2 T s ) ] } ， 1 − α 2 T s ≤ ∣ f ∣ ≤ 1 + α 2 T s 0 ， ∣ f ∣ > 1 + α 2 T s (2.1.14) H(f)=\begin{cases} T_s，0\leq|f|\leq\frac{1-\alpha}{2T_s}\\ \frac{T_s}{2}\{1+cos[\frac{\pi T_s}{\alpha}(|f|-\frac{1-\alpha}{2T_s})]\}，\frac{1-\alpha}{2T_s}\leq|f|\leq\frac{1+\alpha}{2T_s}\\ 0，|f|>\frac{1+\alpha}{2T_s} \tag{2.1.14} \end{cases} H(f)=⎩⎪⎨⎪⎧​Ts​，0≤∣f∣≤2Ts​1−α​2Ts​​{1+cos[απTs​​(∣f∣−2Ts​1−α​)]}，2Ts​1−α​≤∣f∣≤2Ts​1+α​0，∣f∣>2Ts​1+α​​(2.1.14)

其中$\alpha$称为滚动因子，其取值为$0\le \alpha \le 1$。

这里假定信道是理想带限信道，其频率幅度函数为：
$$|C(f)|=Rect(\frac{f}{2B})\text{\hspace{1em}(2.1.15)}$$
由为了让误码率尽量小，$G_R(f)$选用匹配滤波器，与$G_T(f)$匹配，同时，由于信道是带限的，若发送滤波器的冲激响应仍采用方波，那么其频谱$Sa$函数带宽无限，通过信道会损失信息，所以发送信号的带宽不应该超过信道带宽$B$，同时应充分利用信道，故假设$g_T(t)$的带宽为$BW=B$。

推导得到如下式子：
$$|G_R(f)|=|G_T(f)|=\sqrt{H_{rcos}(f)}\text{\hspace{1em}(2.1.16)}$$
又由式子$$(2.1.11)、(2.1.12)$$知，$$h(t)$$为升余弦信号，因此得出发送滤波器的单位冲激响应为：
$$g_T(t)=F^{-1}[\sqrt{H_{rcos}(f)}]\text{\hspace{1em}(2.1.17)}$$

2.2 数字频带传输系统

2.2.1数字系统频带传输模型

数字频带传输系统其实就是在数字基带传输系统的基础上加上了一个带通调制和解调的环节，以此将要传输的信号从基带调制到频带。

图2-4 数字系统频带传输模型

2.2.2信道带宽受限时的$QPSK$传输

$QPSK$名叫四相移相键控，又名四进制移相键控，该信号的正弦载波有$4$个可能的离散相位状态，每个载波相位携带$2$个二进制符号，下面我们先假设信道宽度为无限宽，采用矩形脉冲成形滤波器进行推导，之后再类比到信道带宽受限的情况去。

$QPSK$的信号表示式为：
$$s_i(t)=Acos({\omega }_{c}t+{\theta }_i)i=1,2,3,40\le t\le T_s\text{\hspace{1em}(2.2.1)}$$
其中$T_s$为四进制符号间隔，${\theta }_i(i=1,2,3,4)$为正弦载波的相位，有$4$种可能状态。

这里可以用矢量图表示$QPSK$信号，当${\theta }_i=(i-1)\frac{\pi }{2}$时，如图$2-5(a)$所示，当${\theta }_i=(2i-1)\frac{\pi }{4}$时，如图$2-5(b)$所示。

图2-5 QPSK信号矢量图

我们着重讨论图$2-5(b)$这种情况的$QPSK$信号的产生和解调。

式子$(2.2.1)$可以写成：
$$\begin{gathered} s_i(t)=Acos({\omega }_{c}t+{\theta }_i) \\ =A(cos{\theta }_{i}cos{\omega }_{c}t-sin{\theta }_{i}sin{\omega }_{c}t)0\le t\le T_s\text{\hspace{1em}(2.2.2)} \end{gathered}$$
其中，${\theta }_i$为$\frac{\pi }{4}、\frac{3\pi }{4}、\frac{5\pi }{4}、\frac{7\pi }{4}$，则：
$$cos{\theta }_i=\pm \frac{1}{\sqrt{2}}；sin{\theta }_i=\pm \frac{1}{\sqrt{2}}\text{\hspace{1em}(2.2.3)}$$
于是，式子$(2.2.2)$可以写成：
$$s_i(t)=\frac{A}{\sqrt{2}}[I(t)cos{\omega }_{c}t-Q(t)sin{\omega }_{c}t]0\le t\le T_s\text{\hspace{1em}(2.2.4)}$$
这样可以得到如图$2-6$所示的正交调制框图。

图2-6 产生QPSK信号的正交调制原理图

图$2-6$中的两条支路分别称为同相支路和正交支路。

通过原理图和上面的推导也可以看出$QPSK$信号可以看成是两个$2PSK$信号的叠加，所以在解调时，可对两路信号分别进行$2PSK$的解调，然后进行串并变换，进而得到所传输的数据，同时为了达到更佳的误码率，这里的接收滤波器也采用匹配滤波器，如图$2-7$所示。

图2-7 QPSK匹配滤波器最佳接收

由于信道频带有限，使得输出的信号频带有限，这样信号的时域就是无限长的了，就可能产生码间干扰，通过前面$2.1.3$中的讨论知，为了避免码间干扰，需要满足发送滤波器和匹配滤波器都是根升余弦滤波器，即满足式子$(2.1.16)、(2.1.17)$，同样假设信道传输过程中的噪声为加性高斯白噪声$n_w(t)$，得到信道带宽受限时的$QPSK$传输系统如图$2-8$所示

图2-8 信道带宽受限时的QPSK传输系统

2.2.3信道带宽受限时的$16QAM$传输

$QAM$名叫正交幅度调制，是由两个正交载波多电平振幅键控信号叠加而成的，他与$MPSK$的不同之处在于两个支路的多电平幅度序列是相互独立的，$MQAM$的信号表达式为：
$$s_{QAM}(t)=a_{i_c}{g}_T(t)cos{\omega }_{c}t-a_{i_s}{g}_T(t)sin{\omega }_{c}ti=1,2,\dots ,M0\le t\le T_s\text{\hspace{1em}(2.2.6)}$$
式中， { a i c } \{a_{i_c}\} {aic​​}及 { a i s } \{a_{i_s}\} {ais​​}是一组离散电平的集合，$g_T(t)$是基带成形滤波器的冲激响应

$MQAM$信号波形也可以表示为两个归一化正交基函数的线性组合，即：
$$s_i(t)=s_{i_1}{f}_1(t)+s_{i_2}{f}_2(t)i=1,2,\dots ,M0\le t\le T_s\text{\hspace{1em}(2.2.7)}$$
其中，两个归一化正交基函数为：
$$f_1(t)=\sqrt{\frac{2}{E_g}}{g}_T(t)cos{\omega }_{c}ti=1,2,\dots ,M\text{\hspace{1em}(2.2.8)}$$

$$f_2(t)=-\sqrt{\frac{2}{E_g}}{g}_T(t)sin{\omega }_{c}ti=1,2,\dots ,M\text{\hspace{1em}(2.2.9)}$$

系数
$$s_{i_1}={\int }_0^{T_s}{s}_i(t)f_1(t)dt=a_{i_c}\sqrt{\frac{E_g}{2}}i=1,2,\dots ,M\text{\hspace{1em}(2.2.10)}$$

$$s_{i_2}={\int }_0^{T_s}{s}_i(t)f_2(t)dt=a_{i_s}\sqrt{\frac{E_g}{2}}i=1,2,\dots ,M\text{\hspace{1em}(2.2.11)}$$

因此$MQAM$信号波形的二维矢量表示为：
$$\overrightarrow{s_i}=[s_{i_1},s_{i_2}]=[a_{i_c}\sqrt{\frac{E_g}{2}},a_{i_s}\sqrt{\frac{E_g}{2}}]i=1,2,\dots ,M\text{\hspace{1em}(2.2.12)}$$
其中$E_g$为脉冲$g_T(t)$的能量。

又由于矩形星座的$MQAM$信号尽管不是最优的星座结构，但在满足给定的最小欧氏距离条件下，其所需要平均发送的功率仅比最优$MQAM$的星座结构的信号平均发送功率稍大，而矩形星座的$MQAM$信号的产生及解调在实际中比较容易实现，因此这里采用矩形星座的$MQAM$信号。

而对于$M=2^K$，且$K$为偶数的矩形星座的$MQAM$信号，可等效为同相及正交支路的$\sqrt{M}$进制$ASK$信号之和，每个支路具有$\sqrt{M}=2^{\frac{K}{2}}$个信号电平。

这样，对于矩形星座的$16QAM$信号，就可以等效为同相及正交支路的四进制$ASK$信号之和，每个支路具有个$4$信号电平，其大小为：
$$a_{i_c},a_{i_s}=2i-1-4i=1,2,3,4\text{\hspace{1em}(2.2.13)}$$
则矩形星座的$16QAM$信号的信号空间图和产生框图如图$2-9$和图$2-10$所示

图2-9 矩形16QAM信号空间图

图2-10 矩形星座16QAM信号的产生框图

由于信道带宽受限，因此选择其发送成形滤波器为根升余弦滤波器，其窄脉冲信号为不归零脉冲，在加性高斯白噪声信道条件下，其最佳接收框图如图$2-11$所示。

图2-11 矩形星座16QAM信号的最佳接收框图

三、误码性能分析

3.1 数字基带传输系统

3.1.1 信道带宽无限时的单极性基带传输

​ 对于单极性信号，信号通过$AWGN$信道到达接收滤波器前的信号如式子$(2.1.3)$为$r(t)=s(t)+n_w(t)$，其中$s(t)=s_i(t)$满足式子$(2.1.2)$，$n_w(t)$为高斯白噪声。又根据式子$(2.1.4)$，这里令匹配滤波器与$s_1(t)$相匹配，得到：
$$y(t)={\int }_0^{t}r(\tau )s_1(T-t+\tau )d\tau \text{\hspace{1em}(3.1.1)}$$
​ 在$s(t)$波形结束的时候进行取值，即$t_0=T_s$，代入上式得：
$$y\equiv y(T_s)={\int }_0^{T_s}[s(\tau )+n_w(\tau )]s_1(\tau )d\tau ={\int }_0^{T_s}s(\tau )s_1(\tau )d\tau +{\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau \text{\hspace{1em}(3.1.2)}$$

• 当输入信号为$s(t)=s_1(t)$（发$“1”$）时：
  $$y\equiv y(T_s)={\int }_0^{T_s}[s_1(\tau ){]}^{2}d\tau +{\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau \text{\hspace{1em}(3.1.3)}$$

• 当输入信号为$s(t)=s_2(t)$（发$“1”$）时：
  $$y\equiv y(T_s)={\int }_0^{T_s}{s}_2(\tau )s_1(\tau )d\tau +{\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau \text{\hspace{1em}(3.1.4)}$$

令$Z={\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau$推导得其满足$Z\sim N(0,\frac{N_0}{2}{E}_{b1})$，其中$N_0$为高斯白噪声信号$n_w(\tau )$的半功率谱密度，$E_{b1}$为信号$s_1(t)$的能量。

而式子${\int }_0^{T_s}[s_1(\tau ){]}^{2}d\tau =E_{b1}$，又由于$s_2(t)=0$，因此可以得出${\int }_0^{T_s}{s}_2(\tau )s_1(\tau )d\tau =0$，就得到了判决变量$y$的概率分布为：

• 发"1"：

$$y\sim N(E_{b1},\frac{N_0}{2}{E}_{b1})\text{\hspace{1em}(3.1.5)}$$

• 发"0"：

$$y\sim N(0,\frac{N_0}{2}{E}_{b1})\text{\hspace{1em}(3.1.6)}$$

令${\sigma }_n^2=\frac{N_0}{2}{E}_{b1}$，我们可以得到发$“1”$以及发$“0”$时，判决变量的条件概率密度函数分别为：
$$p(y|"0")=\frac{1}{\sqrt{2\pi {\sigma }_n^2}}{e}^{-\frac{y^2}{2{\sigma }_n^2}}\text{\hspace{1em}(3.1.7)}$$

$$p(y|"1")=\frac{1}{\sqrt{2\pi {\sigma }_n^2}}{e}^{-\frac{(y-E_{b1}{)}^2}{2{\sigma }_n^2}}\text{\hspace{1em}(3.1.8)}$$

其示意图如图$3-1$（图中$A=E_{b1}$）：

图3-1 单极性判决变量条件概率密度函数示意图

假设信源发$“1”$和发$“0”$的概率是一样的，都等于$0.5$，这样由图$3-1$很容易看出，最佳门限应为：
$$V_T=\frac{A}{2}=\frac{E_{b1}}{2}\text{\hspace{1em}(3.1.9)}$$
​ 显然，从图$3-1$来看，阴影部分的面积，$P_r(y\le V_T|“1”)$，也就是发$“1”$的错误译码的概率$P_{e1}$，它可以进行如下计算：
$$\begin{array}{rl}{P}_{e1}& =P_r(y\le V_T|"1")\\ & ={\int }_{-\infty }^{V_T}p(y|"1")dy\\ & ={\int }_{-\infty }^{V_T}\frac{1}{\sqrt{2\pi {\sigma }_n^2}}{e}^{-\frac{(y-E_{b1}{)}^2}{2{\sigma }_n^2}}dy\\ & ={\int }_{-\infty }^{V_T}\frac{1}{\sqrt{2\pi {\sigma }_n^2}}{e}^{-\frac{(y-E_{b1}/{\sigma }_n{)}^2}{2}}d(\frac{y-E_{b1}}{{\sigma }_n})\\ & ={\int }_{-\infty }^{\frac{V_T-E_{b1}}{{\sigma }_n}}\frac{1}{\sqrt{2\pi }}{e}^{\frac{-z^2}{2}}dz\\ & ={\int }_{\frac{E_{b1}-V_T}{{\sigma }_n}}^{+\infty }\frac{1}{\sqrt{2\pi }}{e}^{\frac{-z^2}{2}}dz\\ & =Q(\frac{E_{b1}-V_T}{{\sigma }_n})\end{array}\text{\hspace{1em}(3.1.10)}$$
同时可以得出：
$$P_{e0}=Q(\frac{V_T}{{\sigma }_n})\text{\hspace{1em}(3.1.11)}$$
​ 将公式$(3.1.9)$和公式${\sigma }_n^2=\frac{N_0}{2}{E}_{b1}$代入得：
$$P_{e1}=Q(\sqrt{\frac{E_{b1}}{2N_0}})\text{\hspace{1em}(3.1.12)}$$

$$P_{e0}=Q(\sqrt{\frac{E_{b1}}{2N_0}})\text{\hspace{1em}(3.1.13)}$$

计算得误码率为：
$$P_e=P_0{P}_{e0}+P_1{P}_{e1}=0.5\times Q(\sqrt{\frac{E_{b1}}{2N_0}})+0.5\times Q(\sqrt{\frac{E_{b1}}{2N_0}})=Q(\sqrt{\frac{E_{b1}}{2N_0}})\text{\hspace{1em}(3.1.14)}$$
又对于单极性码，计算符号的平均能量为$$E_b=0.5\times E_{b1}$$，代入得：
$$P_e=Q(\sqrt{\frac{E_b}{N_0}})\text{\hspace{1em}(3.1.15)}$$

3.1.2 信道带宽受限时的双极性基带传输

​ 对于双极性信号，由于我们假设输入信号的带宽等于信道的带宽，因此信号通过$AWGN$信道到达接收滤波器前的信号如式子$(2.1.8)$为$r(t)=s(t)+n_w(t)$，其中$s(t)=s_i(t)$满足式子$(2.1.7)$，$n_w(t)$为高斯白噪声。又这里的接收滤波器也采用的是匹配滤波器，所以令匹配滤波器与$s_1(t)$相匹配，得到：
$$y(t)={\int }_0^{t}r(\tau )s_1(T-t+\tau )d\tau \text{\hspace{1em}(3.1.16)}$$
​ 在$s(t)$波形结束的时候进行取值，即$t_0=T_s$，代入上式得：
$$y\equiv y(T_s)={\int }_0^{T_s}[s(\tau )+n_w(\tau )]s_1(\tau )d\tau ={\int }_0^{T_s}s(\tau )s_1(\tau )d\tau +{\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau \text{\hspace{1em}(3.1.17)}$$

• 当输入信号为$s(t)=s_1(t)$（发$“1”$）时：
  $$y\equiv y(T_s)={\int }_0^{T_s}[s_1(\tau ){]}^{2}d\tau +{\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau \text{\hspace{1em}(3.1.18)}$$

• 当输入信号为$s(t)=s_2(t)$（发$“1”$）时：
  $$y\equiv y(T_s)={\int }_0^{T_s}{s}_2(\tau )s_1(\tau )d\tau +{\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau \text{\hspace{1em}(3.1.19)}$$

令$Z={\int }_0^{T_s}{n}_w(\tau )s_1(\tau )d\tau$推导得其满足$Z\sim N(0,\frac{N_0}{2}{E}_{b1})$，其中$N_0$为高斯白噪声信号$n_w(\tau )$的半功率谱密度，$E_{b1}$为信号$s_1(t)$的能量。

而式子${\int }_0^{T_s}[s_1(\tau ){]}^{2}d\tau =E_{b1}$，又由于$s_2(t)=-s_1(t)$，因此可以得出${\int }_0^{T_s}{s}_2(\tau )s_1(\tau )d\tau =-E_{b1}$，就得到了判决变量$y$的概率分布为：

• 发"1"：

$$y\sim N(E_{b1},\frac{N_0}{2}{E}_{b1})\text{\hspace{1em}(3.1.20)}$$

• 发"0"：

$$y\sim N(-E_{b1},\frac{N_0}{2}{E}_{b1})\text{\hspace{1em}(3.1.21)}$$

令${\sigma }_n^2=\frac{N_0}{2}{E}_{b1}$，我们可以得到发$“1”$以及发$“0”$时，判决变量的条件概率密度函数分别为：
$$p(y|"0")=\frac{1}{\sqrt{2\pi {\sigma }_n^2}}{e}^{-\frac{(y+E_{b1}{)}^2}{2{\sigma }_n^2}}\text{\hspace{1em}(3.1.22)}$$

$$p(y|"1")=\frac{1}{\sqrt{2\pi {\sigma }_n^2}}{e}^{-\frac{(y-E_{b1}{)}^2}{2{\sigma }_n^2}}\text{\hspace{1em}(3.1.23)}$$

其示意图如图$3-2$（图中$A=E_{b1}$）：

图3-2 双极性判决变量条件概率密度函数示意图

假设信源发$“1”$和发$“0”$的概率是一样的，都等于$0.5$，这样由图$3-1$很容易看出，最佳门限应为：
$$V_T=0\text{\hspace{1em}(3.1.24)}$$
​ 显然，从图$3-2$来看，阴影部分的面积，$P_r(y\le V_T|“1”)$，也就是发$“1”$的错误译码的概率$P_{e1}$，根据前面单极性码误码率的推导，将${\sigma }_n^2=\frac{N_0}{2}{E}_{b1}$代入得：
$$P_{e1}=Q(\sqrt{\frac{2E_{b1}}{N_0}})\text{\hspace{1em}(3.1.25)}$$
同理，可以得到：
$$P_{e0}=Q(\sqrt{\frac{2E_{b1}}{N_0}})\text{\hspace{1em}(3.1.26)}$$
计算出符号的平均能量为$E_b=E_{b1}$，推导出此时的误码率为：
$$P_e=P_0{P}_{e0}+P_1{P}_{e1}=0.5\times Q(\sqrt{\frac{2E_b}{N_0}})+0.5\times Q(\sqrt{\frac{2E_b}{N_0}})=Q(\sqrt{\frac{2E_b}{N_0}})\text{\hspace{1em}(3.1.27)}$$

3.2 数字频带传输系统

3.2.1 信道带宽受限时的$QPSK$传输

在给定二进制信息速率的条件下，$QPSK$的同相支路及正交支路的四进制符号速率是二进制信息速率的一半，即$T_s=2T_b$。在给定信号总发送功率的条件下，$QPSK$同相支路或正交支路的信号功率是总发送功率的一半，于是在加性高斯白噪声信道的条件下，得到两支路的平均错误概率为：
$$P_{eI}=P_{eQ}=Q(\sqrt{\frac{(\frac{A}{\sqrt{2}}{)}^2(2T_b)}{2N_0}})=Q(\sqrt{\frac{2E_b}{N_0}})\text{\hspace{1em}(3.2.1)}$$
其中$E_b={\int }_0^{T_s}[s_1(\tau ){]}^{2}d\tau$，是平均比特能量，$N_0$为加性噪声的单边功率谱密度。

当$QPSK$发送信源输出的二进制符号$"1"$和$"0"$等概率出现时，二进制码元经串并变换后在同相支路及正交支路也是等该分布的，所以在收端的同相及正交支路解调的输出经并串变换后的数据，其总的平均误比特率与$I$支路或$Q$支路的平均误判概率是相同的，即：
$$P_b=P_I{P}_{eI}+P_Q{P}_{eQ}\text{\hspace{1em}(3.2.2)}$$
其中$$P_I=P_Q=\frac{1}{2}$$，因此$$QPSK$$的平均误比特率为：
$$P_b=P_{eI}=P_{eQ}=Q(\sqrt{\frac{2E_b}{N_0}})\text{\hspace{1em}(3.2.3)}$$
这也说明，对于$QPSK$而言，在$QPSK$与$2PSK$的输入二进制信息速率相同，二者的发送功率相同，加性噪声的单边功率谱密度$N_0$相同的条件下，$QPSK$与$2PSK$的平均误比特率是相同的。
而式中，由于符号速率为原二进制序列速率的一半，因此每个符号携带两个二进制数据，因此每个二进制数据的平均能量为：
$$E_b=\frac{1}{2}{E}_s\text{\hspace{1em}(3.2.4)}$$
下面推导$E_s$计算公式：

每个符号$s_i(t)$的能量为：
E s = ∫ 0 T s s i 2 ( t ) d t = ∫ 0 T s { A 2 [ I ( t ) c o s ω c t − Q ( t ) s i n ω c t ] } 2 d t = ∫ 0 T s A 2 2 { [ I ( t ) c o s ω c t ] 2 + Q ( t ) I ( t ) c o s ω c t s i n ω c t + [ Q ( t ) s i n ω c t ] 2 } d t (3.2.5) \begin{aligned} E_s&=\int_{0}^{T_s}{s_i^2(t)}dt\\ &=\int_{0}^{T_s}{\{\frac{A}{\sqrt{2}}[I(t)cos\omega_ct-Q(t)sin\omega_ct]\}^2dt}\\ &=\int_{0}^{T_s}{\frac{A^2}{2}\{[I(t)cos\omega_ct]^2+Q(t)I(t)cos\omega_ctsin\omega_ct+[Q(t)sin\omega_ct]^2\}}dt\\ \end{aligned} \tag{3.2.5} Es​​=∫0Ts​​si2​(t)dt=∫0Ts​​{2 ​A​[I(t)cosωc​t−Q(t)sinωc​t]}2dt=∫0Ts​​2A2​{[I(t)cosωc​t]2+Q(t)I(t)cosωc​tsinωc​t+[Q(t)sinωc​t]2}dt​(3.2.5)
根据正交性：
$${\int }_0^{T_s}Q(t)I(t)cos{\omega }_{c}tsin{\omega }_{c}tdt=0\text{\hspace{1em}(3.2.6)}$$
因此得：
E s = ∫ 0 T s A 2 2 { [ I ( t ) c o s ω c t ] 2 + [ Q ( t ) s i n ω c t ] 2 } d t = A 2 4 ∫ 0 T s I 2 ( t ) + Q 2 ( t ) d t = A 2 2 ∫ 0 T s I 2 ( t ) d t (3.2.7) \begin{aligned} E_s&=\int_{0}^{T_s}{\frac{A^2}{2}\{[I(t)cos\omega_ct]^2+[Q(t)sin\omega_ct]^2\}}dt\\ &=\frac{A^2}{4}\int_{0}^{T_s}{I^2(t)+Q^2(t)}dt\\ &=\frac{A^2}{2}\int_{0}^{T_s}{I^2(t)}dt\\ \end{aligned} \tag{3.2.7} Es​​=∫0Ts​​2A2​{[I(t)cosωc​t]2+[Q(t)sinωc​t]2}dt=4A2​∫0Ts​​I2(t)+Q2(t)dt=2A2​∫0Ts​​I2(t)dt​(3.2.7)
其中${\int }_0^{T_s}{I}^2(t)dt={\int }_0^{T_s}{Q}^2(t)dt$为发送滤波器单位冲激信号的能量。这样结合式子$3.2.4$得到：
$$E_b=\frac{1}{2}{E}_s=\frac{A^2}{4}{\int }_0^{T_s}{I}^2(t)dt\text{\hspace{1em}(3.2.8)}$$

3.2.2 信道带宽受限时的$16QAM$传输

矩形星座$MQAM$的最小欧氏距离为：
$$d_{min}=\sqrt{2E_g}=\sqrt{\frac{6E_{b}lo{g}_2^M}{M-1}}\text{\hspace{1em}(3.2.9)}$$
其中，$E_g$为脉冲$g_T(t)$的能量，$E_b$表示符号能量。

此时的判决门限为：
$$V_T=d_{min}\text{\hspace{1em}(3.2.10)}$$
矩形星座$MQAM$的最佳接收误符率与$MASK$的一样，取决于数字基带$MPAM$的误符率。$MQAM$的正确判决符号的概率为：
$$P_c=(1-P_{\sqrt{M}}{)}^2\text{\hspace{1em}(3.2.11)}$$
式中，$P_{\sqrt{M}}$表示同相或正交支路$\sqrt{M}$进制$ASK$的误符率，该$\sqrt{M}$进制$ASK$的平均功率是$MQAM$信号总的平均功率$P_s$的一半，即：
$$\begin{array}{rl}{P}_{\sqrt{M}}& =2(1-\frac{1}{\sqrt{M}})Q(\sqrt{\frac{3P_s{T}_s}{(M-1)N_0}})\\ & =2(1-\frac{1}{\sqrt{M}})Q(\sqrt{\frac{3E_s}{(M-1)N_0}})\\ & =2(1-\frac{1}{\sqrt{M}})Q(\sqrt{\frac{d_{min}^2}{2N_0}})\end{array}\text{\hspace{1em}(3.2.12)}$$

而式子中$E_s$表示平均能量，即：
$$E_i={\int }_0^{T_s}{s}_i^2(t)dti=1,2,\dots ,M\text{\hspace{1em}(3.2.13)}$$

$$E_s=\frac{1}{M}\sum _{i=1}^M{E}_i\text{\hspace{1em}(3.2.14)}$$

从而得，$MQAM$的误符率为：
$$P_M=1-P_c=1-(1-P_{\sqrt{M}}{)}^2=2P_{\sqrt{M}}-P_{\sqrt{M}}^2\text{\hspace{1em}(3.2.15)}$$

$$P_M\approx 2P_{\sqrt{M}}=4(1-\frac{1}{\sqrt{M}})Q(\sqrt{\frac{3E_s}{(M-1)N_0}})\text{\hspace{1em}(3.2.16)}$$

当采用格雷码，且$E_b/N_0$较大时，$MQAM$的误比特率近似等于误符率除以$lo{g}_2^M$，即：
$$P_b\approx \frac{P_M}{lo{g}_2^M}\text{\hspace{1em}(3.2.17)}$$
这样，$16QAM$的误符率和误比特率就为：
$$P_{16}\approx 4(1-\frac{1}{4})Q(\sqrt{\frac{3E_s}{(16-1)N_0}})=3Q(\sqrt{\frac{E_s}{5N_0}})\text{\hspace{1em}(3.2.18)}$$

$$P_b\approx \frac{P_{16}}{lo{g}_2^{16}}\text{\hspace{1em}(3.2.19)}$$

四、仿真模型与仿真过程

4.1 数字基带传输系统

4.1.1 信道带宽无限时的单极性基带传输

取比特传输速率$r_b=100bit/s$，持续时间为：$10s$，符号传输周期为$T_s=10ms$，$A=1$，采样频率为$f_{s}ample=1000Hz$，采用二进制编码且为单极性归零码，同时通过计算与仿真结果比较，选取噪声单边功率谱密度为$N_0={10}^{-3}$时，得到的仿真效果比较好，代码如下：

1. 系统参数设置

   %------------------
   %系统参数设置
   %-----------------
   T_start=0;%开始时间
   T_stop=10;%截止时间
   T=T_stop-T_start;%仿真持续时间
   T_sample=1/1000;%采样间隔
   f_sample=1/T_sample; % 采样速率
   N_sample=T/T_sample;% 采样点数
   r_b=100;%比特传输速率
   T_b=1/r_b;%比特传输周期
   r_s=r_b;%符号传输速率
   T_s=T_b;%符号传输周期
   Bits=2;%是二进制码
   NumMa=T*r_s;%在仿真时间内总的传输符号数
   NumBits=NumMa*log2(2);%在仿真时间内总的传输的bit数（由于本仿真为二进制，因此=NumMa，下面有的地方有部分混用）
   A=1;%设置gt(t)的幅度值为A
   
   

2. 发射端

   %------------------
   %发射端
   %-----------------
   %randi()函数生成均匀分布的伪随机整数，范围为imin--imax，如果没指定imin，则默认为1
   b1=randi([0,1],NumMa,1);%b1是列向量，输入序列
   %由于实际抽样得到的序列比需要的符号数多，因此需要进行处理
   b2=zeros(f_sample/r_s,NumMa);%f_sample/r_s：1s抽样符号数是1s所需传输符号数的几倍
   b2(1,:)=b1;%b2矩阵第一列等于b1向量
   b=reshape(b2,1,f_sample/r_s*NumMa);%转换一维向量{bn}
   %通过窄脉冲调制后的序列
   a=A*b;%通过二进制编码得到{an}，也可以看成是通过窄脉冲调制得到的一堆冲击响应序列
   g_T = ones(1,T_s/T_sample);%发送滤波器
   s = conv(a,g_T);%信号通过发送滤波器，发送端发送信号
   

3. 无限带宽$AWGN$信道

   %--------------------------------
   %无限带宽AWGN信道
   %--------------------------------
   rng(7,'twister');%随机因子
   N_0=10^(-3);%单边功率谱密度
   noise_w=wgn(1,length(s),N_0/2*f_sample,'linear');%产生白噪声
   r=s+noise_w;%接收端接收信号
   

4. 接收端

   %------------------
   %接收端
   %------------------
   V_T=0.5*A*A*T_s/T_sample;%判决门限
   %采用匹配滤波器
   g_R = g_T;%接收滤波器与s1(t)匹配
   yt = conv(r,g_R);
   %抽样
   y = zeros(1,NumMa);%采样得到的序列
   for i = 1:NumMa
       k = T_s/T_sample;%时域采样间隔
       x = i * k;
       y(i)=yt(x);
   end
   %判决
   b0=zeros(1,NumMa);%输出序列
   for i=1:NumMa
       if(y(i)<V_T)
           b0(i)=0;
       else
           b0(i)=1;
       end
   end
   

5. 画图

   %------------------
   %画图
   %------------------   
   figure(1)
   %------输入序列，通过窄脉冲调制后的序列，发送端发送信号,接收端接收信号，采样得到的序列，输出序列------
   %b{n}，输入序列
   subplot(3,2,1);
   nb1=0:length(b1)-1;
   stem(nb1,b1,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b{n}');
   title('输入序列')
   
   %a{n}，通过窄脉冲调制后的序列
   subplot(3,2,3)
   na=0:length(a)-1;
   stem(na,a,'*');
   axis([0,100,-0.2,1.2]);
   xlabel('n');
   ylabel('a{n}');
   title('通过窄脉冲调制后的序列')
   
   %s(t)，发送端发送信号
   subplot(3,2,5)
   ns=0:length(s)-1;
   plot(ns,s);
   axis([0,100,-0.2,1.2]);
   xlabel('t');
   ylabel('s(t)');
   title('发送端发送信号')
   
   %r(t)，接收端接收信号
   subplot(3,2,2)
   nr=0:length(r)-1;
   plot(nr,r);
   axis([0,100,-2,2]);
   xlabel('t');
   ylabel('r(t)');
   title('接收端接收信号')
   
   %y，采样得到的序列
   subplot(3,2,4)
   ny=0:length(y)-1;
   stem(ny,y,'*');
   axis([0,10,-20,20]);
   xlabel('n');
   ylabel('y');
   title('采样得到的序列')
   
   %b0{n}，输出序列
   subplot(3,2,6)
   nb0=0:length(b0)-1;
   stem(nb0,b0,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b^~{n}');
   title('输出序列')
   

6. 误码率

   %------------------
   %误码率
   %------------------ 
   BER_shiji=length(find(b1~= b0'))/NumBits%实际误比特率
   Es1 = A*A*T_s;
   Eb=0.5*Es1;
   BER_lilun=qfunc(sqrt(Eb/(N_0)))%理论误比特率
   

4.1.2 信道带宽受限时的双极性基带传输

取比特传输速率$r_b=100bit/s$，持续时间为：$10s$，符号传输周期为$T_s=10ms$，$A=1$，采样频率为$1000Hz$，采用二进制编码且为双极性归零码，同时通过计算与仿真结果比较，选取噪声单边功率谱密度为$N_0={8\times 10}^{-5}$时，得到的仿真效果比较好，代码如下：

1. 系统参数设置

   %------------------
   %系统参数设置
   %-----------------
   T_start=0;%开始时间
   T_stop=10;%截止时间
   T=T_stop-T_start;%仿真持续时间
   T_sample=1/1000;%采样间隔
   f_sample=1/T_sample; % 采样速率
   N_sample=T/T_sample;% 采样点数
   r_b=100;%比特传输速率
   T_b=1/r_b;%比特传输周期
   r_s=r_b;%符号传输速率
   T_s=T_b;%符号传输周期
   Bits=2;%是二进制码
   NumMa=T*r_s;%在仿真时间内总的传输符号数
   NumBits=NumMa*log2(2);%在仿真时间内总的传输的bit数（由于本仿真为二进制，因此=NumMa，下面有的地方有部分混用）
   A=1;%设置gt(t)的幅度值为A
   NumCoff=40%number of coefficients of RRC(阶数)(对根升余弦滤波器的冲激函数进行40点取样)
   alpha=0.25;%df=alpha*rs=25Hz
   

2. 发射端

   %------------------
   %发射端
   %-----------------
   %randi()函数生成均匀分布的伪随机整数，范围为imin--imax，如果没指定imin，则默认为1
   b1_1=randi([0,1],NumMa,1);%b1是列向量,输入序列
   %b1 = sign(b1_1-0.5);%转换成-1和1（双极性码）
   b1 = 2*b1_1-1;%转换成-1和1（双极性码）
   %由于实际抽样得到的序列比需要的符号数多，因此需要进行处理(内插)(窄脉冲滤波器)
   b2=zeros(f_sample/r_s,NumMa);%f_sample/r_s：1s抽样符号数是1s所需传输符号数的几倍
   b2(1,:)=b1;%b2矩阵第一列等于b1向量
   b=reshape(b2,1,f_sample/r_s*NumMa);%转换一维向量{bn}，冲击串
   %通过窄脉冲调制后的序列
   a=A*b;%通过二进制编码得到{an}，也可以看成是通过窄脉冲调制得到的一堆冲击响应序列
   %根升余弦滤波器 阶数，截止频率，过度带宽，采样频率
   g_T=firrcos(NumCoff,r_s/2,alpha*r_s,f_sample);%发送滤波器
   s = conv(a,g_T);%信号通过发送滤波器，发送端发送信号
   

3. 无限带宽$AWGN$信道

   %--------------------------------
   %有限带宽AWGN信道
   %--------------------------------
   rng(7,'twister');%随机因子
   N_0=8*10^(-5);%单边功率谱密度
   noise_w=wgn(1,length(s),N_0/2*f_sample,'linear');%产生白噪声
   r=s+noise_w;%接收端接收信号
   

4. 接收端

   %------------------
   %接收端
   %------------------
   V_T=0;%判决门限
   %采用匹配滤波器
   g_R = g_T;%接收滤波器与s1(t)匹配
   yt = conv(r,g_R);
   %抽样和判决
   sample1=zeros(f_sample/r_s,NumBits);
   sample1(1,:)=ones(1,NumBits);
   sample2=reshape(sample1,1,f_sample/r_s*NumBits);
   sample3=zeros(1,length(yt));
   sample3(NumCoff+1:NumCoff+f_sample/r_s*NumBits)=sample2;
   y=yt.*sample3;%采样得到的序列
   y(:,all(y==0,1))=[];%去掉多余的零
   b0=(sign(y)+1)*0.5;%判决，输出序列
   

5. 画图

   %------------------
   %画图
   %------------------   
   figure(1)
   %------输入序列，通过窄脉冲调制后的序列，发送端发送信号,接收端接收信号，采样得到的序列，输出序列------
   %b{n}，输入序列
   subplot(3,2,1);
   nb1_1=0:length(b1_1)-1;
   stem(nb1_1,b1_1,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b{n}');
   title('输入序列')
   
   %a{n}，通过窄脉冲调制后的序列
   subplot(3,2,3)
   na=0:length(a)-1;
   stem(na,a,'*');
   axis([0,100,-1.2,1.2]);
   xlabel('n');
   ylabel('a{n}');
   title('通过窄脉冲调制后的序列')
   
   %s(t)，发送端发送信号
   subplot(3,2,5)
   ns=0:length(s)-1;
   plot(ns,s);
   axis([0,100,-0.2,0.2]);
   xlabel('t');
   ylabel('s(t)');
   title('发送端发送信号')
   
   %r(t)，接收端接收信号
   subplot(3,2,2)
   nr=0:length(r)-1;
   plot(nr,r);
   axis([0,100,-0.5,0.5]);
   xlabel('t');
   ylabel('r(t)');
   title('接收端接收信号')
   
   %y，采样得到的序列
   subplot(3,2,4)
   ny=0:length(y)-1;
   stem(ny,y,'*');
   axis([0,10,-0.2,0.2]);
   xlabel('n');
   ylabel('y');
   title('采样得到的序列')
   
   %b0{n}，输出序列
   subplot(3,2,6)
   nb0=0:length(b0)-1;
   stem(nb0,b0,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b^~{n}');
   title('输出序列')
   

6. 误码率

   %------------------
   %误码率
   %------------------ 
   BER_shiji=length(find(b1_1~= b0'))/NumBits%实际误比特率
   E_NTs = 0:NumCoff;
   E_Ts = E_NTs*T_sample;
   Es1 = trapz(E_Ts,g_T.^2);
   BER_lilun = qfunc(sqrt((2*Es1)/N_0))%理论误比特率
   

4.1.3 误比特率与信噪比曲线

定义信噪比为：
$$SNR=\frac{E_b}{N_0}\text{\hspace{1em}(4.1.1)}$$
这样结合上面的式子$(3.1.15)$和式子$(3.1.27)$得，此时两个情况的误比特率分别为：
$$P_{e1}=Q(\sqrt{SNR})\text{\hspace{1em}(4.1.2)}$$

$$P_{e2}=Q(\sqrt{2SNR})\text{\hspace{1em}(4.1.3)}$$

将上面的两个代码改成函数，并主代码进行循环调用即可得到基带传输的理论与实际的误比特率与信噪比曲线。

代码如下：

%基带传输SNR
sn=0.1:0.1:30;%信噪比序列
pe1_1=[];%用于存实际误码率
pe1_2=[];%用于存理论误码率
pe2_1=[];%用于存实际误码率
pe2_2=[];%用于存理论误码率
for i=1:length(sn)
    [pe1_shiji,pe1_lilun]=test4_4_1(sn(i));
    pe1_1=[pe1_1,pe1_shiji];
    pe1_2=[pe1_2,pe1_lilun];
    [pe2_shiji,pe2_lilun]=test4_4_3(sn(i));
    pe2_1=[pe2_1,pe2_shiji];
    pe2_2=[pe2_2,pe2_lilun];
end
snr=10*log10(sn);%将信噪比转换为dB
%pe1=qfunc(sqrt(sn));
%pe2=qfunc(sqrt(2.*sn));
figure(2);
semilogy(snr,pe1_2);
hold on;
semilogy(snr,pe1_1,'*');
hold on;
semilogy(snr,pe2_2);
hold on;
semilogy(snr,pe2_1,'*');
hold on;
axis([-10,10,5*10^-3,0.5])
legend('单极性理论误比特率','单极性实际误比特率','双极性理论误比特率','双极性实际误比特率');
title('平均误比特率曲线');
xlabel('SNR(dB)')
ylabel('误比特率')

4.2 数字频带传输系统

4.2.1 信道带宽受限时的$QPSK$传输

取比特传输速率$r_b=100bit/s$，持续时间为：$10s$，符号传输周期为$T_s=10ms$，$A=1$，采样频率为$2000Hz$，采用四进制编码，同时通过计算与仿真结果比较，选取噪声单边功率谱密度为$N_0=2\times 10^{-6}$时，得到的仿真效果比较好，代码如下：

1. 系统参数设置

   %------------------
   %系统参数设置
   %-----------------
   T_start=0;%开始时间
   T_stop=10;%截止时间
   T=T_stop-T_start;%仿真持续时间
   T_sample=1/2000;%采样间隔
   f_sample=1/T_sample; %采样速率
   N_sample=T/T_sample;%采样点数
   r_b=100;%比特传输速率
   T_b=1/r_b;%比特传输周期
   r_s=r_b;%符号传输速率
   T_s=T_b;%符号传输周期
   NumMa=T*r_s;%在仿真时间内总的传输符号数
   NumBits=NumMa*log2(2);%在仿真时间内总的传输的bit数
   A=1;%设置gt(t)的幅度值为A
   NumCoff=40%number of coefficients of RRC(阶数)(对根升余弦滤波器的冲激函数进行40点取样)
   alpha=0.25;%df=alpha*rs=25Hz
   

2. 发射端

   %------------------
   %发射端
   %-----------------
   %randi()函数生成均匀分布的伪随机整数，范围为imin--imax，如果没指定imin，则默认为1
   b1_1=randi([0,1],NumMa,1);%b1是列向量，输入序列
   %b1 = sign(b1_1-0.5);%转换成-1和1（双极性码）
   b1 = 2*b1_1-1;%转换成-1和1（双极性码）
   %-----串并变换-----
   bI0=[];%串并变换后得到的偶序列
   bQ0=[];%串并变换后得到的奇序列
   for i=1:NumMa
       if(mod(i,2)==1)
           bQ0=[bQ0,b1(i)];
       else
           bI0=[bI0,b1(i)];
       end
   end
   r_s=r_s/2;
   bQ=upsample(bQ0,f_sample/r_s);%内插，相当于进行窄脉冲调制
   bI=upsample(bI0,f_sample/r_s);%内插，相当于进行窄脉冲调制
   bQ=A*bQ;%可以看成是通过窄脉冲调制得到的一堆冲击响应序列
   bI=A*bI;%可以看成是通过窄脉冲调制得到的一堆冲击响应序列
   %根升余弦滤波器 阶数，截止频率，过渡带宽，采样频率
   g_T=firrcos(NumCoff,r_s/2,alpha*r_s,f_sample);%发送滤波器
   sI = conv(bI,g_T);%信号通过发送滤波器
   sQ = conv(bQ,g_T);%信号通过发送滤波器
   t=0:f_sample*T-1+40;%NumCoff=40
   t=t.*T_sample;
   c_t1=cos(2*pi*400*t);%取载波频率为400Hz(31.25<400<1000)满足抽样定理fc+31.25<=fmax<=2000/2=1000
   c_t2=-sin(2*pi*400*t);%取载波频率为400Hz(31.25<400<1000)满足抽样定理fc+31.25<=fmax<=2000/2=1000
   s1=sI.*c_t1;%频谱搬移，I支路发送信号
   s2=sQ.*c_t2;%频谱搬移，Q支路发送信号
   s=s1+s2;%发送端发送信号
   

3. 无限带宽$AWGN$信道

   %--------------------------------
   %有限带宽AWGN信道
   %--------------------------------
   N_0=2*10^(-6);%单边功率谱密度
   rng(6,'twister');%加了一个随机数种子
   noise_w=wgn(1,length(s),N_0/2*f_sample,'linear');%产生白噪声
   r=s+noise_w;%接收端接收信号
   

4. 接收端

   %------------------
   %接收端
   %------------------
   rQ=r.*c_t2; 
   rI=r.*c_t1;
   V_T=0;%判决门限
   %采用匹配滤波器
   g_R = g_T;%接收滤波器与s1(t)匹配
   yQt=conv(rQ,g_R);
   yIt=conv(rI,g_R);
   %-----抽样-----
   NumMa=NumMa/2;%符号数减半
   sample1=zeros(f_sample/r_s,NumMa);
   sample1(1,:)=ones(1,NumMa);
   sample2=reshape(sample1,1,f_sample/r_s*NumMa);
   sample4=sample2;%抽样函数
   %由于g_T的中心不再原点，因此应该进行一定的平移
   sample3=zeros(1,length(yQt));
   sample3(NumCoff+1:NumCoff+f_sample/r_s*NumMa)=sample2;
   y1=yQt.*sample3;%采样得到的Q支路序列
   sample3=zeros(1,length(yIt));
   sample3(NumCoff+1:NumCoff+f_sample/r_s*NumMa)=sample4;
   y2=yIt.*sample3;%采样得到的I支路序列
   %-----判决-----
   y1(:,all(y1==0,1))=[];
   b1=(sign(y1)+1)*0.5;%判决得到的Q支路序列
   y2(:,all(y2==0,1))=[];
   b2=(sign(y2)+1)*0.5;%判决得到的I支路序列
   %-----并串变换-----
   b0=[];%输出序列
   for i=1:length(y2)
       b0=[b0,b1(i)];
       b0=[b0,b2(i)];
   end
   

5. 画图

   %------------------
   %画图
   %------------------
   figure(1);
   %------输入序列，串并变换后得到的偶序列，串并变换后得到的奇序列------
   %b{n}，输入序列
   subplot(3,1,1);
   nb1=1:length(b1_1);
   stem(nb1,b1_1,'*');
   axis([1,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b{n}');
   title('输入序列');
   
   %bI0，串并变换后得到的偶序列
   subplot(3,1,2);
   nbI0=1:length(bI0);
   stem(nbI0,bI0,'*');
   axis([1,10,-1.2,1.2]);
   xlabel('k');
   ylabel('b2n[k]');
   title('串并变换后得到的偶序列');
   
   %bQ0，串并变换后得到的奇序列
   subplot(3,1,3);
   nbQ0=1:length(bQ0);
   stem(nbQ0,bQ0,'*');
   axis([1,10,-1.2,1.2]);
   xlabel('k');
   ylabel('b2n-1[k]');
   title('串并变换后得到的奇序列');
   
   figure(2);
   %------I支路发送信号,Q支路发送信号，发送端发送信号------
   %s1(t)，I支路发送信号
   subplot(3,1,1);
   ns1=0:length(s1)-1;
   plot(ns1,s1);
   axis([0,600,-0.035,0.035]);
   xlabel('t');
   ylabel('sI(t)');
   title('I支路发送信号');
   
   %s2(t)，Q支路发送信号
   subplot(3,1,2);
   ns2=0:length(s2)-1;
   plot(ns2,s2);
   axis([0,600,-0.035,0.035]);
   xlabel('t');
   ylabel('sQ(t)');
   title('Q支路发送信号');
   
   %s(t)，发送端发送信号
   subplot(3,1,3);
   ns=0:length(s)-1;
   plot(ns,s);
   axis([0,600,-0.035,0.035]);
   xlabel('t');
   ylabel('s(t)');
   title('发送端发送信号');
   
   figure(3);
   %------接收端接收信号，采样得到的I支路序列，采样得到的Q支路序列------
   %r(t)，接收端接收信号
   subplot(3,1,1);
   nr=0:length(r)-1;
   plot(nr,r);
   axis([0,3000,-0.2,0.2]);
   xlabel('t');
   ylabel('r(t)');
   title('接收端接收信号');
   
   %y2，采样得到的I支路序列
   subplot(3,1,2);
   ny2=0:length(y2)-1;
   stem(ny2,y2,'*');
   axis([0,10,-0.05,0.05]);
   xlabel('n');
   ylabel('yI');
   title('采样得到的I支路序列');
   
   %y1，采样得到的Q支路序列
   subplot(3,1,3);
   ny1=0:length(y1)-1;
   stem(ny1,y1,'*');
   axis([0,10,-0.05,0.05]);
   xlabel('n');
   ylabel('yQ');
   title('采样得到的Q支路序列');
   
   figure(4);
   %------判决得到的I支路序列，判决得到的Q支路序列，输出序列------
   %b2{n}，判决得到的I支路序列
   subplot(3,1,1);
   nb2=1:length(b2);
   stem(nb2,b2,'*');
   axis([1,10,-0.2,1.2]);
   xlabel('n');
   ylabel('bI^~{n}');
   title('判决得到的I支路序列');
   
   %b1{n}，判决得到的Q支路序列
   subplot(3,1,2)
   nb1=1:length(b1);
   stem(nb1,b1,'*');
   axis([1,10,-0.2,1.2]);
   xlabel('n');
   ylabel('bQ^~{n}');
   title('判决得到的Q支路序列');
   
   %b0{n}，输出序列
   subplot(3,1,3)
   nb0=1:length(b0);
   stem(nb0,b0,'*');
   axis([1,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b^~{n}');
   title('输出序列');
   

6. 误码率

   %------------------
   %误码率
   %------------------ 
   BER_shiji=length(find(b1_1~= b0'))/NumBits%实际误比特率
   E_NTs = 0:NumCoff;
   E_Ts = E_NTs*T_sample;
   Es1 = ((A^2)*trapz(E_Ts,g_T.^2));
   Eb=Es1/2;%Eb=0.5Es1，一个符号携带两个二进制
   BER_lilun = qfunc(sqrt((2*Eb)/N_0))%理论误比特率
   

4.2.2信道带宽受限时的$16QAM$传输

取比特传输速率$r_b=100bit/s$，持续时间为：$10s$，符号传输周期为$T_s=10ms$，$A=1$，采样频率为$2000Hz$，采用十六进制编码，同时通过计算与仿真结果比较，选取噪声单边功率谱密度为$N_0=2\times 10^{-4}$时，得到的仿真效果比较好，代码如下：

1. 系统参数设置

   %------------------
   %系统参数设置
   %-----------------
   T_start=0;%开始时间
   T_stop=10;%截止时间
   T=T_stop-T_start;%仿真持续时间
   T_sample=1/2000;%采样间隔
   f_sample=1/T_sample; % 采样速率
   N_sample=T/T_sample;% 采样点数
   r_b=100;%信息传输速率
   T_b=1/r_b;%信息传输周期
   r_s=r_b;%符号传输速率
   T_s=T_b;%符号传输周期
   NumMa=T*r_s;%在仿真时间内总的传输符号数
   NumBits=NumMa*log2(2);%在仿真时间内总的传输的bit数
   A=1;%设置gt(t)的幅度值为A
   NumCoff=40%number of coefficients of RRC(阶数)(对根升余弦滤波器的冲激函数进行40点取样)
   alpha=0.25;%df=alpha*rs=25Hz
   M=16;%16路
   

2. 发射端

   %------------------
   %发射端
   %-----------------
   %randi()函数生成均匀分布的伪随机整数，范围为imin--imax，如果没指定imin，则默认为1
   b1_1=randi([0,1],NumMa,1);%b1是列向量，输入序列
   %b1 = sign(b1_1-0.5);%转换成-1和1（双极性码）
   %-----串并变换-----
   bI0=[];%输入偶数序列
   bQ0=[];%输入奇数序列
   for i=1:NumMa
       if(mod(i,2)==1)
           bQ0=[bQ0,b1_1(i)];
       else
           bI0=[bI0,b1_1(i)];
       end
   end
   %-----二进制变四进制-----
   %00:-3
   %01:-1
   %11:+1
   %10:+3
   bI1=[];%四进制的同相支路序列
   bQ1=[];%四进制的正交支路序列
   for i=1:2:length(bI0)
       if(bI0(i)==0&&bI0(i+1)==0)
          bI1=[bI1,-3]; 
       elseif(bI0(i)==0&&bI0(i+1)==1)
          bI1=[bI1,-1]; 
       elseif(bI0(i)==1&&bI0(i+1)==1)
          bI1=[bI1,+1];
       else
          bI1=[bI1,+3];
       end
   end
   for i=1:2:length(bQ0)
       if(bQ0(i)==0&&bQ0(i+1)==0)
          bQ1=[bQ1,-3]; 
       elseif(bQ0(i)==0&&bQ0(i+1)==1)
          bQ1=[bQ1,-1]; 
       elseif(bQ0(i)==1&&bQ0(i+1)==1)
          bQ1=[bQ1,+1];
       else
          bQ1=[bQ1,+3];
       end
   end
   bQ=upsample(bQ1,(f_sample/r_s)*4);%内插，相当于进行窄脉冲调制
   bI=upsample(bI1,(f_sample/r_s)*4);%内插，相当于进行窄脉冲调制
   bQ=bQ;%可以看成是通过窄脉冲调制得到的一堆冲击响应序列
   bI=bI;%可以看成是通过窄脉冲调制得到的一堆冲击响应序列
   %-----成形滤波器-----
   r_s=r_b/4;
   T_s=T_b*4;
   %新的根升余弦滤波器 阿尔法 截断的符号范围 单个符号范围的采样个数
   g_T = rcosdesign(alpha,6,T_s/T_sample,'sqrt');
   sI = conv(bI,g_T);%信号通过发送滤波器
   sQ = conv(bQ,g_T);%信号通过发送滤波器
   %-----频谱搬移-----
   t=0:length(sI)-1;
   t=t.*T_sample;
   c_t1=cos(2*pi*400*t);%取载波频率为400Hz(15.625<400<1000)满足抽样定理fc+15.625<=fmax<=2000/2=1000
   t=0:length(sQ)-1;
   t=t.*T_sample;
   c_t2=-sin(2*pi*400*t);%取载波频率为400Hz(15.625<400<1000)满足抽样定理fc+15.625<=fmax<=2000/2=1000
   s1=sI.*c_t1;%频谱搬移，同相支路信号
   s2=sQ.*c_t2;%频谱搬移，正交支路信号
   s=s1+s2;%16QAM信号
   

3. 无限带宽$AWGN$信道

   %--------------------------------
   %有限带宽AWGN信道
   %--------------------------------
   N_0=2*10^(-4);%单边功率谱密度
   rng(3,'twister');%加了一个随机数种子
   noise_w=wgn(1,length(s),N_0/2*f_sample,'linear');%产生白噪声
   r=s+noise_w;%接收机输入信号
   

4. 接收端

   %------------------
   %接收端
   %------------------
   %-----计算得到Eg-----
   E_NTs = 0:length(g_T)-1;
   E_Ts = E_NTs*T_sample;
   Eg = trapz(E_Ts,g_T.^2);
   %-----得到f1,f2-----
   t=0:length(g_T)-1;
   t=t.*T_sample;
   c_t1=cos(2*pi*400*t);%取载波频率为400Hz(15.625<400<1000)满足抽样定理fc+15.625<=fmax<=2000/2=1000
   t=0:length(g_T)-1;
   t=t.*T_sample;
   c_t2=-sin(2*pi*400*t);%取载波频率为400Hz(15.625<400<1000)满足抽样定理fc+15.625<=fmax<=2000/2=1000
   
   f1=sqrt(2/Eg).*g_T.*c_t1;
   f2=sqrt(2/Eg).*g_T.*c_t2;
   
   %相关接收机（匹配滤波器）
   %由于选取的fc是fs整数倍，coswct关于T/2偶对称，sinwct关于T/2奇对称
   yQt=conv(r,-f2);%卷积的时候需要整一个轴对称，因此乘以-f2
   yIt=conv(r,f1);
   %-----抽样-----
   NumMa1=NumMa/4;%由于rs=rb/4，因此总的符号数要减小
   sample1=zeros(f_sample/r_s,NumMa1);
   sample1(1,:)=ones(1,NumMa1);
   sample2=reshape(sample1,1,f_sample/r_s*NumMa1);%抽样函数
   sample4=sample2;%抽样函数
   %由于g_T的中心不再原点，因此应该进行一定的平移
   %g_T引入了暂态响应，因此sample3相当于是忽略前面一部分进行抽样
   sample3=zeros(1,length(yQt));
   sample3(length(f2):length(f2)+f_sample/r_s*NumMa1-1)=sample2;
   y2=yQt.*sample3;
   sample3=zeros(1,length(yIt));
   sample3(length(f1):length(f1)+f_sample/r_s*NumMa1-1)=sample4;
   y1=yIt.*sample3;
   
   %-----判决-----
   %去掉多余的零
   y1(:,all(y1==0,1))=[];%同相支路采样得到的序列
   y2(:,all(y2==0,1))=[];%正交支路采样得到的序列
   %受卷积影响，乘f_sample，dmin=sqrt(2Eg)
   %119和120行匹配滤波器卷积中本来的卷积运算是积分，但序列的卷积是求和
   %而积分比求和多考虑一个最小间隔T_sample
   %因此匹配滤波器卷积少乘了一个T_sample
   %所以就需要判决门限V_T除以一个T_sample，也就是乘上一个f_sample
   V_T = sqrt(2*Eg)*f_sample%判决门限
   yb1 = [1,length(y1)];%I，同相支路判决得到的序列
   yb2 = [1,length(y2)];%Q，正交支路判决得到的序列
   %00:-3
   %01:-1
   %11:+1
   %10:+3
   %-----四进制变二进制-----
   %采取另一种调用数组的方法
   for i=1:length(y1)
       if(y1(i)<-V_T)%00
           yb1(2*i-1)=0;
           yb1(2*i)=0;
       elseif(y1(i)>=-V_T&&y1(i)<0)%01
           yb1(2*i-1)=0;
           yb1(2*i)=1;
       elseif(y1(i)>=0&&y1(i)<V_T)%11
           yb1(2*i-1)=1;
           yb1(2*i)=1;
       else%10
           yb1(2*i-1)=1;
           yb1(2*i)=0;
       end
   end
   for i=1:length(y2)
       if(y2(i)<-V_T)%00
           yb2(2*i-1)=0;
           yb2(2*i)=0;
       elseif(y2(i)>=-V_T&&y2(i)<0)%01
           yb2(2*i-1)=0;
           yb2(2*i)=1;
       elseif(y2(i)>=0&&y2(i)<V_T)%11
           yb2(2*i-1)=1;
           yb2(2*i)=1;
       else%10
           yb2(2*i-1)=1;
           yb2(2*i)=0;
       end
   end
   %-----并串变换-----
   b0=[1,length(yb1)+length(yb2)];%输出序列
   for i=1:length(yb1)
       b0(2*i-1)=yb2(i);
       b0(2*i)=yb1(i);
   end
   
   

5. 画图

   %------------------
   %画图
   %------------------
   figure(1);
   %------输入序列，四进制的同相支路序列，四进制的同相支路序列------
   %b{n}输入序列
   subplot(3,1,1);
   nb1=0:length(b1_1)-1;
   stem(nb1,b1_1,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b{n}');
   title('输入序列')
   
   %bI1{n}四进制的同相支路序列
   subplot(3,1,2);
   nbI1=0:length(bI1)-1;
   stem(nbI1,bI1,'*');
   axis([0,10,-3.5,3.5]);
   xlabel('n');
   ylabel('bI1[n]');
   title('四进制的同相支路序列')
   
   %bI1{n}四进制的正交支路序列
   subplot(3,1,3);
   nbQ1=0:length(bQ1)-1;
   stem(nbI1,bQ1,'*');
   axis([0,10,-3.5,3.5]);
   xlabel('n');
   ylabel('bQ1[n]');
   title('四进制的正交支路序列')
   
   figure(2);
   %------同相支路信号，正交支路信号，16QAM信号------
   %s1(t)同相支路信号
   subplot(3,1,1);
   ns1=0:length(s1)-1;
   plot(ns1,s1);
   axis([0,800,-0.6,0.6]);
   xlabel('t');
   ylabel('s1(t)');
   title('同相支路信号')
   
   %s2(t)正交支路信号
   subplot(3,1,2);
   ns2=0:length(s2)-1;
   plot(ns2,s2);
   axis([0,800,-0.6,0.6]);
   xlabel('t');
   ylabel('s2(t)');
   title('正交支路信号')
   
   %s(t)16QAM信号
   subplot(3,1,3);
   ns=0:length(s)-1;
   plot(ns,s);
   axis([0,800,-0.6,0.6]);
   xlabel('t');
   ylabel('s(t)');
   title('16QAM信号')
   
   figure(3);
   %------接收机接收信号，同相支路采样得到的序列，正交支路采样得到的序列------
   %r(t)接收机接收信号
   subplot(3,1,1);
   nr=0:length(r)-1;
   plot(nr,r);
   axis([0,5000,-1.2,1.2]);
   xlabel('t');
   ylabel('r(t)');
   title('接收机接收信号')
   
   %y1{n}同相支路采样得到的序列
   subplot(3,1,2);
   ny1=0:length(y1)-1;
   stem(ny1,y1,'*');
   axis([0,10,-125,125]);
   xlabel('n');
   ylabel('y1{n}');
   title('同相支路采样得到的序列')
   
   %y2{n}正交支路采样得到的序列
   subplot(3,1,3);
   ny2=0:length(y2)-1;
   stem(ny2,y2,'*');
   axis([0,10,-125,125]);
   xlabel('n');
   ylabel('y2{n}');
   title('正交支路采样得到的序列')
   
   figure(4);
   %------同相支路判决得到的序列，正交支路判决得到的序列，输出序列------
   %yb1{n}同相支路判决得到的序列
   subplot(3,1,1);
   nyb1=0:length(yb1)-1;
   stem(nyb1,yb1,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('yb1{n}');
   title('同相支路判决得到的序列')
   
   %yb2{n}正交支路判决得到的序列
   subplot(3,1,2);
   nyb2=0:length(yb2)-1;
   stem(nyb2,yb2,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('yb2{n}');
   title('正交支路判决得到的序列')
   
   %b0{n}输出序列
   subplot(3,1,3);
   nb0=0:length(b0)-1;
   stem(nb0,b0,'*');
   axis([0,10,-0.2,1.2]);
   xlabel('n');
   ylabel('b0{n}');
   title('输出序列')
   

6. 误码率

   %------------------
   %误码率
   %------------------ 
   BER_shiji=length(find(b1_1~= b0'))/NumBits%实际误比特率
   dmin=sqrt(2*Eg);%最小欧氏距离
   PsqrtM=2*(1-1/sqrt(M))*qfunc(sqrt((dmin^2)/(2*N_0)));%P根号下M
   PM=2*PsqrtM - (PsqrtM)^2;
   BER_lilun = PM/(log2(M))%理论误比特率（估算公式，认为每个码元只错一个bit）
   

4.2.3 误比特率与信噪比曲线

定义信噪比为：
$$SNR=\frac{E_b}{N_0}\text{\hspace{1em}(4.2.1)}$$
这样结合上面的式子$(3.2.3)$、式子$(3.2.11)$、式子$(3.2.14)$和式子$(3.2.16)$可以计算得到这两个情况的理论误比特率，将如上的两个代码改成函数，并主代码进行循环调用即可得到频带传输的理论与实际的误比特率与信噪比曲线。

代码如下：

sn=0.1:0.1:30;%信噪比序列
pe1_1=[];%用于存实际误码率
pe1_2=[];%用于存理论误码率
pe2_1=[];%用于存实际误码率
pe2_2=[];%用于存理论误码率
for i=1:length(sn)
    [pe1_shiji,pe1_lilun]=test4_6_3(sn(i));
    pe1_1=[pe1_1,pe1_shiji];
    pe1_2=[pe1_2,pe1_lilun];
    [pe2_shiji,pe2_lilun]=test4_6_2(sn(i));
    pe2_1=[pe2_1,pe2_shiji];
    pe2_2=[pe2_2,pe2_lilun];
end
snr=10*log10(sn);%将信噪比转换为dB
figure(3);
semilogy(snr,pe1_2);
hold on;
semilogy(snr,pe1_1,'*');
hold on;
semilogy(snr,pe2_2);
hold on;
semilogy(snr,pe2_1,'*');
hold on;
axis([-10,10,1*10^-2,0.5])
legend('QPSK理论误比特率','QPSK实际误比特率','16QAM理论误比特率','16QAM实际误比特率');
title('平均误比特率曲线');
xlabel('SNR(dB)')
ylabel('误比特率')

五、仿真结果分析

5.1数字基带传输系统

5.1.1 信道带宽无限时的单极性基带传输

图5-1 信道带宽无限时的单极性基带传输仿真波形图

图 5-2 信道带宽无限时的单极性基带传输误比特率理论和实际图

分析：从图$5-1$可以看出，尽管加入高斯白噪声信号之后，总的信号波形变得杂乱无章，但是经过检测得到的码和发射机发出的码是一样的，说明这个误码率还是比较小的。对于图$5-2$，其中$BE{R}_{s}hiji$表示实际（也就是仿真）得到的误比特率，而$BE{R}_{l}ilun$表示理论计算（公式$3.1.15$）得到的理论误比特率，两者十分接近，其误差为：
$$相对误差=|\frac{测量值-理论值}{理论值}|=|\frac{0.0130-0.0127}{0.0127}|\approx 2.36$$

$$绝对误差=|测量值-理论值|=|0.0130-0.0127|=3\times {10}^{-4}$$

5.1.2 信道带宽有限时的双极性基带传输

图5-3 信道带宽有限时的双极性基带传输仿真波形图

图 5-4 信道带宽有限时的双极性基带传输误比特率理论和实际图

分析：从图$5-3$可以看出，尽管加入高斯白噪声信号之后，总的信号波形变得杂乱无章，但是经过检测得到的码和发射机发出的码是一样的，说明这个误码率还是比较小的。对于图$5-4$，其中$BE{R}_{s}hiji$表示实际（也就是仿真）得到的误比特率，而$BE{R}_{l}ilun$表示理论计算（公式$3.1.27$）得到的理论误比特率，两者十分接近，其误差为：
$$相对误差=|\frac{测量值-理论值}{理论值}|=|\frac{0.0660-0.0639}{0.0639}|\approx 3.29$$

$$绝对误差=|测量值-理论值|=|0.0660-0.0639|=2.1\times {10}^{-3}$$

5.1.3 误比特率与信噪比曲线

图5-5 理论的误比特率与信噪比仿真曲线图

图5-6 理论和实际的误比特率与信噪比仿真曲线图

由该曲线可以看出：

1. 对于信道带宽无限时的单极性基带传输系统，由图$5-6$可以看出理论计算的误比特率与实际仿真得到的误比特率在信噪比大的地方和小的地方都近似相等，十分接近，可以看出，仿真结果与理论基本一样，同时，随着信噪比的增大，误比特率呈非线性减小。
2. 对于信道带宽有限时的双极性基带传输系统，由图$5-6$可以看出在信噪比较低的时候，其理论计算的误比特率与实际仿真得到的误比特率比较接近，这是由于此时的信噪比比较低，导致这时候的误差比较大，因此相对来看理论值与实际值之间的较小的差距不明显；但是随着信噪比的增大，误比特率理论值与实际值就产生了比较明显的误差，通过仿真发现是由于仿真时取的传输信号的幅值较小，导致当信噪比较小时，实际仿真的误差较大，但是信噪比在一定范围内，算出来的理论误比特率较小，这就导致中间一段的理论与实际误比特率值有一定的误差；但是通过仿真（图中没有呈现）可知，随着信噪比的增大，实际的误比特率迅速减小至零，而理论计算的信噪比也会达到一个很小的值，这个时候两者的误差又变得很小了。
3. 同时通过图$5-5$可以看到，对于带宽有限时的双极性基带传输系统，当信噪比为$0.7$时，此时的误比特率约为$10^{-3}$，与理论带进去计算得到的值基本一样。
4. 整体来看，两曲线都是随着信噪比的增大，误比特率呈非线性减小，但是在信噪比相等的时候，双极性码的误比特率始终要小于单极性码的误比特率，这与理论分析的结论是一样的。

5.2数字频带传输系统

5.2.1 信道带宽受限时的$QPSK$传输

图5-7 信道带宽有限时的QPSK频带传输仿真波形图1

图5-8 信道带宽有限时的QPSK频带传输仿真波形图2

图5-9 信道带宽有限时的QPSK频带传输仿真波形图3

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-doFKS5xZ-1612168580263)(D:\swjtu syf\大三上\通信原理\实验\仿真实验最终\有限带宽QPSK仿真波形4.jpg)]

图5-10 信道带宽有限时的QPSK频带传输仿真波形图4

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Tj5migPy-1612168580264)(C:\Users\pc\AppData\Roaming\Typora\typora-user-images\image-20210131215226895.png)]

图 5-11 信道带宽有限时的QPSK频带传输误比特率理论和实际图

分析：图$5-7$至图$5-10$中选取噪声单边功率谱密度为$N_0=2\times 10^{-6}$比较小，得到的波形和传输之后解调出来的信号都比较理想，可以看出，尽管加入高斯白噪声信号之后，总的信号波形变得杂乱无章，但是经过检测得到的码和发射机发出的码是一样的，说明这个误比特率还是比较小的。图5-11中$BE{R}_{s}hiji$表示实际（也就是仿真）得到的误比特率，而$BE{R}_{l}ilun$表示理论计算（公式$3.2.3$）得到的理论误比特率，两者十分接近，其误差为：
$$相对误差=|\frac{测量值-理论值}{理论值}|=|\frac{0.0150-0.0142}{0.0142}|\approx 5.63$$

$$绝对误差=|测量值-理论值|=|0.0150-0.0142|=8\times 10^{-4}$$

5.2.2 信道带宽受限时的$16QAM$传输

图5-12 信道带宽有限时的16QAM频带传输仿真波形图1

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-rME0FZZw-1612168580265)(D:\swjtu syf\大三上\通信原理\实验\仿真实验最终\有限带宽16QAM仿真波形2.jpg)]

图5-13 信道带宽有限时的16QAM频带传输仿真波形图2

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-WvZaNrZY-1612168580266)(D:\swjtu syf\大三上\通信原理\实验\仿真实验最终\有限带宽16QAM仿真波形3.jpg)]

图5-14 信道带宽有限时的16QAM频带传输仿真波形图3

图5-15 信道带宽有限时的16QAM频带传输仿真波形图4

图 5-16 信道带宽有限时的16QAM频带传输误比特率理论和实际图

图 5-17 信道带宽有限时的16QAM频带传输实际脉冲能量和判决门限图

分析：图$5-12$至图$5-15$中选取噪声单边功率谱密度为$N_0=2\times 10^{-4}$，得到的波形和传输之后解调出来的信号都比较理想，可以看出，尽管加入高斯白噪声信号之后，总的信号波形变得杂乱无章，但是经过检测得到的码和发射机发出的码是一样的，说明这个误比特率还是比较小的。图$5-16$中$BE{R}_{s}hiji$表示实际（也就是仿真）得到的误比特率，而$BE{R}_{l}ilun$表示理论计算（公式$3.2.3$）得到的理论误比特率，两者十分接近，其误差为：
$$相对误差=|\frac{测量值-理论值}{理论值}|=|\frac{0.0420-0.0409}{0.0409}|\approx 2.69$$

$$绝对误差=|测量值-理论值|=|0.0420-0.0409|=1.1\times 10^{-3}$$

图$5-17$中，$E_g$表示脉冲$g_T(t)$的能量，$V_T$表示的是判决门限，这里注意的是，$V_T$的计算不仅根据公式$(3.2.9)$和公式$(3.2.10)$，还要乘上$f_s$，分析这是由于，卷积运本来应该是积分，但仿真中进行的是序列的卷积，也就是求和，而积分比求和多考虑一个最小时间间隔$T_s$，这样匹配滤波器卷积后得到的序列相当于少乘了一个最小时间间隔$T_s$，因此，计算判决门限时，就需要除以这个最小时间间隔$T_s$，就相当于乘以$f_s$。根据仿真结果可以看出，上面的分析、公式的改进是正确的。

5.2.3 误比特率与信噪比曲线

图5-18 理论和实际的误比特率与信噪比仿真曲线图

由该曲线可以看出：

1. 对于道带宽有限时的$QPSK$频带传输系统，由图$5-17$可以看出理论计算的误比特率与实际仿真得到的误比特率在信噪比大的地方和小的地方都近似相等，十分接近，可以看出，仿真结果与理论基本一样，同时，随着信噪比的增大，误比特率呈非线性减小。
2. 对于道带宽有限时的$16QAM$频带传输系统，由图5-18可以看出在信噪比较高的时候，其理论计算的误比特率与实际仿真得到的误比特率比较接近，而在信噪比比较低的时候，两者出现了很大的差别。这是由于，$MQAM$的计算误比特率的公式（公式$(3.2.16)$）是一个近似的公式，其成立条件是使用格雷码而且一定要信噪比较大的时候，因为这个时候可以近似认为每个码元出错只能到相邻的两个码元上，这样错误的总比特数就近似等于错误的码元数，从而得到公式$(3.2.16)$。仿真结果也正印证了这个条件的必要性，与理论知识相符合。
3. 整体来看，两曲线都是随着信噪比的增大，误比特率呈非线性减小，同时以仿真结果来看，当信噪比相同的时候（对于$16QAM$，只看信噪比比较大的时候），$16QAM$的误比特率要略大于$QPSK$的误比特率。

六、结论

由实际仿真结果看出，四种传输系统在噪声为高斯白噪声的时候的误比特率是相对较小的，都还能比较准确地传输信号。

基带传输和频带传输系统各有优劣，基带传输系统相对比较简单，便于仿真、实现，而频带传输系统能够充分利用频带资源，也为无线电的远程传输提供了可能。

带宽有限和带宽无限时候采用的成形滤波器也不一样，带宽无限的时候，为了方便在时域观察波形，仿真的时候经常采用矩形波；而带宽有限的时候，由于频带有限了，那么时域就变成类似于$S_a$函数的样子，这个时候就需要考虑码间干扰，因此采用根升余弦滤波器，在避免码间干扰的同时，可以一定程度上改进频带利用了（改进$\alpha$）。

由仿真得到的误比特率与理论的误比特率还是有一定的差别，分析原因可能是仿真的时候是以离散的点来代替连续的信号，噪声信号、$s(t)$的波形都有一定的失真，从而导致产生误差，可以通过增加抽样的点数使其能更好地模拟连续信号，以此来减小误差。

由仿真和理论分析的单双极性平均误比特率曲线可以看出，在信噪比相等的时候，双极性码的误比特率要小于单极性码的误比特率，这是由于，双极性码本身传输的符号为$+1$和$-1$，比单极性码的$0$和$+1$之间的差距更大，因此受到噪声影响更小，从而传输后得到的抽样信号$y$的两种情况的均值距离更大，使得平均误码率更小。

在仿真和理论分析的道带宽有限时的$16QAM$频带传输系统的时候也印证了理论推导出来的其误比特率公式的前提条件（格雷码、信噪比较大）的正确性和必要性。

本次实验也能够反映出理论一定要与实践相结合，在实践中检验理论，理解其中的原理，不断改进，才能更加深刻地理解理论。

七、参考文献

• [1] Tanghonghanhaoli.现代通信原理10.1：带宽无限信道下采用低通滤波器（LPF）接收时的误码性能分析[OL].2020-11-15]
• [2]Tanghonghanhaoli.现代通信原理9.2：数字基带传输系统模型[OL]. 2019-11-06