超级干货：BPSK/QPSK数字调制系统误码率MATLAB仿真

本文旨在通过简单实例来对基于MATLAB的数字调制解调系统仿真进行一个较为全面的介绍，并加深对一些基础知识的理解。且有详细解释大多数人在进行数字调制MATLAB仿真时遇到的大部分问题。

数字调制是把数字信号即二进制比特或一定长度的比特串变换为连续信号波形，进而在连续信道中传输的过程。其向信道编译码器提供了一条透明的离散编码信道。大多数情况下，我们将调制过程分成以下几个子过程：

1. 比特（串）到星座图中信号的映射（mapping）。（离散信号） 2. 离散信号（可以为一维实信号或二维复信号）向波形（waveform）的转换。（连续基带信号） 3. 基带波形与载波相乘变换到射频。（可能先到中频再上变频到射频）

所以一般情况下，不进行特殊说明，调制指的是过程1，即映射。映射包含了两层含义：

后者也称为标识（labeling）方式，当星座图确定时，映射指的也就是标识，不过以标识称之更加准确。 M M M阶调制存在 M ! M! M! 种可能的标识方案，标识方案对系统的性能有很大的影响。最常见的标识方式为 Gray 标识，其特点是星座图中邻近的两个信号，其标识仅有一位不同，这样当误符号率一定时，误比特率（BER）最小。

上述星座图并不是唯一的，也可以有任意的旋转角度。

对于BPSK调制，比特到信号的映射关系 s = μ ( b ) s=\mu(b) s=μ(b)可以表示为： s = 1 − 2 b , b = 0 , 1 s=1-2b, b=0, 1 s=1−2b,b=0,1 对于QPSK，映射 s = μ ( b 1 , b 0 ) s=\mu (b_1, b_0) s=μ(b1​,b0​)可以表示为： s = ( 1 − 2 b 0 ) + j ( 1 − 2 b 1 ) 2 s=\dfrac{(1-2b_0)+j(1-2b_1)}{\sqrt{2}} s=2 ​(1−2b0​)+j(1−2b1​)​ 利用上面两个映射关系，就可以实现自己调制，而不用使用MATLAB中的函数或工具箱中的函数（或Object）。

以下代码完成调制功能：

注意：BPSK调制尽量不要使用MATLAB中自带函数。 如pskmod函数和comm.BPSKModulator，你可以自行尝试一下，调制的结果是复数形式，虽然虚部为0。但是当你使用awgn函数加噪声时，该函数会给两个维度（实部和虚部）都加上噪声，而不是只给一个维度上加，这必然导致实际的噪声功率为设想的一半。换句话说，仿真结果相比理论情况有3dB的增益，不过这是错误的。

E b E_b Eb​表示信息比特能量， E s E_s Es​表示符号能量

在一个星座图中，每一个符号 s s s都有其符号能量 E s = ∣ s ∣ 2 E_s=|s|^2 Es​=∣s∣2,该星座图的平均符号能量为： E ‾ s = ∑ m = 0 M − 1 ∣ s m ∣ 2 \overline E_s=\sum_{m=0}^{M-1}|s_m|^2 Es​=m=0∑M−1​∣sm​∣2

N o N_o No​表示高斯白噪声的单边功率谱密度

S N R SNR SNR表示每一个采样时刻信号功率与噪声功率的比值 S N R = σ s 2 σ n 2 SNR=\dfrac{\sigma_s^2}{\sigma_n^2} SNR=σn2​σs2​​

注意此处为每个符号所携带的信息比特数，因为在存在信道编码的场景下，符号所携带的还有冗余校验比特。对于采用 M M M阶调制，以及码率为 R R R的信道编码的情况，上面公式中的 i n f o r m a t i o n _ b i t s _ p e r _ s y m b o l = log ⁡ 2 M ∗ R information\_bits\_per\_symbol=\log_2^M*R information_bits_per_symbol=log2M​∗R

有了以上3点， E b N o , E s N o , S N R \dfrac{E_b}{N_o}, \dfrac{E_s}{N_o}, SNR No​Eb​​,No​Es​​,SNR之间的转化就变得很简单。

对于通信系统性能的仿真，我们均是选择 E b N o \frac{E_b}{N_o} No​Eb​​，可以将其视为归一化的信噪比。通信的目的是传输有效信息，在信道条件给定的情况下，不管系统采用何种调制方式、何种纠错码，我们总是希望平均为每一比特所花费的能量是最少的。所以我们总是选择 E b N o \frac{E_b}{N_o} No​Eb​​ 来进行仿真和比较。

上面提到我们更喜欢使用 E b N o \frac{E_b}{N_o} No​Eb​​，那么这里就以 E b N o \frac{E_b}{N_o} No​Eb​​为例进行说明。

使用awgn函数

r = awgn(s, snr, 'measured');

上面一行代码实现给发送端的信号 s s s加上信噪比 S N R = s n r SNR=snr SNR=snr的AWGN噪声。详情可参考MATLAB文档。那么现在的问题就是如何将 E b N o \frac{E_b}{N_o} No​Eb​​转化为 S N R SNR SNR。下面分别以BPSK与QPSK为例进行分析：

可以看出在将 E s N o \frac{E_s}{N_o} No​Es​​转化为 S N R SNR SNR时，QPSK与BPSK时不同的。这在前面已经解释了：

先进行如下符号约定：

L ( s l k ) = ln ⁡ p ( s l k = 0 ∣ r l ) p ( s l k = 1 ∣ r l ) = ln ⁡ p ( s l k = 0 , r l ) p ( s l k = 1 , r l ) = ln ⁡ ∑ s l ∈ χ 0 k p ( r l ∣ s l ) ⋅ p ( s l ) ∑ s l ∈ χ 1 k p ( r l ∣ s l ) ⋅ p ( s l ) = ln ⁡ ∑ s l ∈ χ 0 k [ p ( r l ∣ s l ) ⋅ ∏ k = 0 K − 1 p ( s l k ) ] ∑ s l ∈ χ 1 k [ p ( r l ∣ s l ) ⋅ ∏ k = 0 K − 1 p ( s l k ) ] \begin{aligned} L(s_l^k) &=\ln\dfrac{p\left(s_l^k=0\mid r_l \right)}{p\left(s_l^k=1\mid r_l \right)} \\ &=\ln\dfrac{p\left(s_l^k=0, r_l \right)}{p\left(s_l^k=1, r_l \right)} \\ &=\ln\dfrac{\sum\limits_{s_l\in \chi_0^k}p\left( r_l\mid s_l\right)\cdot p\left(s_l\right)}{\sum\limits_{s_l\in \chi_1^k}p\left( r_l\mid s_l\right)\cdot p\left(s_l\right)}\\ &=\ln\dfrac{\sum\limits_{s_l\in \chi_0^k} \left[ p\left( r_l\mid s_l\right)\cdot \prod\limits_{k=0}^{K-1}p\left(s_l^k\right) \right]}{\sum\limits_{s_l\in \chi_1^k} \left[ p\left( r_l\mid s_l\right)\cdot \prod\limits_{k=0}^{K-1}p\left(s_l^k\right) \right]} \end{aligned} L(slk​)​=lnp(slk​=1∣rl​)p(slk​=0∣rl​)​=lnp(slk​=1,rl​)p(slk​=0,rl​)​=lnsl​∈χ1k​∑​p(rl​∣sl​)⋅p(sl​)sl​∈χ0k​∑​p(rl​∣sl​)⋅p(sl​)​=lnsl​∈χ1k​∑​[p(rl​∣sl​)⋅k=0∏K−1​p(slk​)]sl​∈χ0k​∑​[p(rl​∣sl​)⋅k=0∏K−1​p(slk​)]​​ 由于 L a ( b ) = p ( b = 0 ) p ( b = 1 ) L_a(b)=\dfrac{p(b=0)}{p(b=1)} La​(b)=p(b=1)p(b=0)​，所以： p ( b = 0 ) = e L a ( b ) 1 + e L a ( b ) p ( b = 1 ) = 1 1 + e L a ( b ) \begin{aligned} p(b=0)=\dfrac{e^{L_a(b)}}{1+e^{L_a(b)}} \\ p(b=1)=\dfrac{1}{1+e^{L_a(b)}} \end{aligned} p(b=0)=1+eLa​(b)eLa​(b)​p(b=1)=1+eLa​(b)1​​ 综上可得： L ( s l k ) = ln ⁡ ∑ s l ∈ χ 0 k [ p ( r l ∣ s l ) e x p ( ∑ k = 0 K − 1 s l k ‾ ⋅ L a ( s l k ) ) ] ∑ s l ∈ χ 1 k [ p ( r l ∣ s l ) e x p ( ∑ k = 0 K − 1 s l k ‾ ⋅ L a ( s l k ) ) ] L(s_l^k) =\ln\dfrac{\sum\limits_{s_{l}\in\chi_{0}^{k}}\left[ p\left(r_l|s_l\right) {\rm exp}\left( \sum\limits_{k=0}^{K-1}\overline{s_l^k}\cdot L_a\left(s_l^k \right) \right) \right] }{\sum\limits_{s_l\in\chi_{1}^{k}}\left[ p\left(r_l|s_l\right) {\rm exp}\left( \sum\limits_{k=0}^{K-1}\overline{s_l^k}\cdot L_a\left(s_l^k \right) \right)\right]} L(slk​)=lnsl​∈χ1k​∑​[p(rl​∣sl​)exp(k=0∑K−1​slk​​⋅La​(slk​))]sl​∈χ0k​∑​[p(rl​∣sl​)exp(k=0∑K−1​slk​​⋅La​(slk​))]​ 之所以推导MAP解调器公式，是因为这个模块不仅可用于硬判决、软判决，也可以当作一个软输入软输出解调器用在迭代检测系统中。对于前两种情况，不存在先验信息，或认为 L a = 0 L_a=0 La​=0，即并不知道某个比特更倾向于为0或为1。上式可以简化为： L ( s l k ) = ln ⁡ ∑ s l ∈ χ 0 k p ( r l ∣ s l ) ∑ s l ∈ χ 1 k p ( r l ∣ s l ) L(s_l^k) =\ln\dfrac{\sum\limits_{s_{l}\in\chi_{0}^{k}}p\left(r_l|s_l\right) }{\sum\limits_{s_{l}\in\chi_{1}^{k}}p\left(r_l|s_l\right)} L(slk​)=lnsl​∈χ1k​∑​p(rl​∣sl​)sl​∈χ0k​∑​p(rl​∣sl​)​ 对于无信道编码的系统，直接进行硬判决得到信息比特： s l k ^ = 1 − s i g n ( L ( s l k ) ) 2 \hat{s_l^k}=\dfrac{1-sign(L(s_l^k))}{2} slk​^​=21−sign(L(slk​))​ 对于有信道编码的系统，可以进行硬判决，然后将硬判决结果传给译码器；也可以直接将似然比交予软输入译码器进行译码。前者因为丢失信息，往往有2dB左右的性能损失。

前文提到不建议用MATLAB中自带函数或Object进行BPSK调制，并说明了原因。对于解调，我们也完全可以用以下3行代码完成：