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

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

数字调制的概念

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

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

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

• 星座图
• 星座图中各信号点对应的比特串

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

BPSK与QPSK仿真实例

星座图与标识方式

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

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

调制模块代码

以下代码完成调制功能：

%% BPSK Modulation
frameLength = 10000;
b = randi([0, 1], frameLength, 1);
s_bpsk = 1 - 2 * b;

%% QPSK Modulation

% method 1
frameLength = 10000; % bits per data frame
b = randi([0, 1], frameLength, 1);
qpskMod = comm.QPSKModulator('BitInput', true, 'PhaseOffset', pi/4, 'SymbolMapping', 'Gray');
s_qpsk_1 = qpskMod(b);

% method 2
frameLength = 10000; % bits per data frame
M = 4; % Modulation Order
symLength = frameLength / M; % symbols per data frame
sym = randi([0, M-1], symLength, 1);
s_qpsk_2 = pskmod(sym, M, pi/4, 'gray');

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

添加AWGN噪声

$\frac{E_b}{N_o},\frac{E_s}{N_o},SNR$的辨析

基本概念

1. $E_b$表示信息比特能量， $E_s$表示符号能量

   在一个星座图中，每一个符号$s$都有其符号能量$E_s=|s{|}^2$,该星座图的平均符号能量为：
   $${\overline{E}}_s=\sum _{m=0}^{M-1}|s_m{|}^2$$

2. $N_o$表示高斯白噪声的单边功率谱密度

3. $SNR$表示每一个采样时刻信号功率与噪声功率的比值
   $$SNR=\frac{{\sigma }_s^2}{{\sigma }_n^2}$$

相互之间的关系

1. $E_b$与$E_s$之间的关系
   $$E_b=\frac{E_s}{information\_bits\_per\_symbol}$$

注意此处为每个符号所携带的信息比特数，因为在存在信道编码的场景下，符号所携带的还有冗余校验比特。对于采用$M$阶调制，以及码率为$R$的信道编码的情况，上面公式中的$information\_bits\_per\_symbol={\log }_2^M\ast R$

2. $N_o$与${\sigma }_n^2$的关系
   • 对于实信号，即一维调制（BPSK，PAM）
     $$N_o=2{\sigma }_n^2$$
   • 对于复信号，即IQ调制
     $$N_o={\sigma }_n^2$$
     $${\sigma }_{nI}^2={\sigma }_{nQ}^2=\frac{{\sigma }_n^2}{2}=\frac{N_o}{2}$$
     综上，可知对任意维度，都有
     $$N_o=2{\sigma }^2$$
     推导过程在附录中给出。
3. $E_s$与${\sigma }_s^2$的关系
   在不存在上下采样的情况下（大多数时候），显然
   $$E_s={\sigma }_s^2$$
   若每个符号采样$N$点，$N=\frac{T_{sym}}{T_s}$那么有：
   $$E_s=N{\sigma }_s^2$$

有了以上3点，$\frac{E_b}{N_o},\frac{E_s}{N_o},SNR$之间的转化就变得很简单。

$\frac{E_b}{N_o},\frac{E_s}{N_o},SNR$的选择

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

MATLAB中添加AWGN噪声

上面提到我们更喜欢使用$\frac{E_b}{N_o}$，那么这里就以$\frac{E_b}{N_o}$为例进行说明。

1. 使用awgn函数

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

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

	% BPSK
	frameLength = 1e6;
	EbNo_dB = 5;
	M = 2; % BPSK
	K = log2(M); % coded bits per symbol
	R = 1; %code rate, information bits per coded bit
	infoBitsPerSymbol = K * R;
	EsNo_dB = EbNo_dB + 10*log10(infoBitsPerSymbol);
	snr = EsNo_dB + 10*log10(2);
	b = randi([0,1],frameLength,1);
	s = 1-2*b;
	r = awgn(s, snr, 'measured');

% QPSK
EbNo_dB = 5;
frameLength = 1e6;
M = 4; % QPSK
K = log2(M); % coded bits per symbol
R = 1; %code rate, information bits per coded bit
infoBitsPerSymbol = K * R;
EsNo_dB = EbNo_dB + 10*log10(infoBitsPerSymbol);
snr = EsNo_dB;
b = randi([0,1],frameLength,1);
qpskMod = comm.QPSKModulator('BitInput', true, 'PhaseOffset', pi/4, 'SymbolMapping', 'Gray');
s = qpskMod(b);
r = awgn(s, snr, 'measured');

可以看出在将$\frac{E_s}{N_o}$转化为$SNR$时，QPSK与BPSK时不同的。这在前面已经解释了：

• QPSK
  $$SNR=\frac{{\sigma }_s^2}{{\sigma }_n^2}=\frac{E_s}{{\sigma }_{nI}^2+{\sigma }_{nQ}^2}=\frac{E_s}{N_o}$$
• BPSK
  $$SNR=\frac{{\sigma }_s^2}{{\sigma }_n^2}=\frac{E_s}{{\sigma }_{nI}^2}=\frac{2E_s}{N_o}$$

2. 自己手动添加

%BPSK
EbNo_dB = 5;
M = 2; % BPSK
K = log2(M); % coded bits per symbol
R = 1; %code rate, information bits per coded bit
infoBitsPerSymbol = K * R;
EsNo_dB = EbNo_dB + 10*log10(infoBitsPerSymbol);
EsNo_Lin = 10^(EsNo_dB/10);
Es = 1; % PSK or other constellation(e.g. QAM) with unit average symbol energy
No = Es / EsNo_Lin;
sigma = sqrt(No/2); % each dimension
b = randi([0,1],frameLength,1);
s = 1 - 2 * b;
r = s + sigma*randn(size(s));

% QPSK
EbNo_dB = 5;
M = 4; % QPSK
K = log2(M); % coded bits per symbol
R = 1; %code rate, information bits per coded bit
infoBitsPerSymbol = K * R;
EsNo_dB = EbNo_dB + 10*log10(infoBitsPerSymbol);
EsNo_Lin = 10^(EsNo_dB/10);
Es = 1; % PSK or other constellation(e.g. QAM) with unit average symbol energy
No = Es / EsNo_Lin;
sigma = sqrt(No/2); % each dimension
b = randi([0,1],frameLength,1);
qpskMod = comm.QPSKModulator('BitInput', true, 'PhaseOffset', pi/4, 'SymbolMapping', 'Gray');
s = qpskMod(b);
r = s + sigma*(randn(size(s)) + 1j * randn(size(s)));

最大后验（MAP）解调算法

先进行如下符号约定：

符号	含义
$s_l$	发送的第$l$个符号
$r_l$	接收到的第$l$个符号
${\chi }_b^k$	第$k$个比特为$b$的符号子集，$b=0,1$
$s_l^k$	发送的第$l$个符号中的第$k$比特
$\overline{s_l^k}$	$s_l^k$取反
$L_a(b)$	比特$b$的先验似然比，$L_a(b)=\frac{p(b=0)}{p(b=1)}$
$L(b)$	比特$b$的后验似然比，$L(b)=\frac{p(b=0\mid r)}{p(b=1\mid r)}$

$$\begin{array}{rl}L(s_l^k)& =\ln \frac{p\left(s_l^k=0\mid r_l\right)}{p\left(s_l^k=1\mid r_l\right)}\\ & =\ln \frac{p\left(s_l^k=0,r_l\right)}{p\left(s_l^k=1,r_l\right)}\\ & =\ln \frac{\sum _{s_l\in {\chi }_0^k}p\left(r_l\mid s_l\right)\cdot p\left(s_l\right)}{\sum _{s_l\in {\chi }_1^k}p\left(r_l\mid s_l\right)\cdot p\left(s_l\right)}\\ & =\ln \frac{\sum _{s_l\in {\chi }_0^k}\left[p\left(r_l\mid s_l\right)\cdot \prod _{k=0}^{K-1}p\left(s_l^k\right)\right]}{\sum _{s_l\in {\chi }_1^k}\left[p\left(r_l\mid s_l\right)\cdot \prod _{k=0}^{K-1}p\left(s_l^k\right)\right]}\end{array}$$
由于 $L_a(b)=\frac{p(b=0)}{p(b=1)}$，所以：
$$\begin{array}{r}p(b=0)=\frac{e^{L_a(b)}}{1+e^{L_a(b)}}\\ p(b=1)=\frac{1}{1+e^{L_a(b)}}\end{array}$$
综上可得：
$$L(s_l^k)=\ln \frac{\sum _{s_l\in {\chi }_0^k}\left[p\left(r_l|s_l\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\sum _{k=0}^{K-1}\overline{s_l^k}\cdot L_a\left(s_l^k\right)\right)\right]}{\sum _{s_l\in {\chi }_1^k}\left[p\left(r_l|s_l\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\sum _{k=0}^{K-1}\overline{s_l^k}\cdot L_a\left(s_l^k\right)\right)\right]}$$
之所以推导MAP解调器公式，是因为这个模块不仅可用于硬判决、软判决，也可以当作一个软输入软输出解调器用在迭代检测系统中。对于前两种情况，不存在先验信息，或认为 $L_a=0$，即并不知道某个比特更倾向于为0或为1。上式可以简化为：
$$L(s_l^k)=\ln \frac{\sum _{s_l\in {\chi }_0^k}p\left(r_l|s_l\right)}{\sum _{s_l\in {\chi }_1^k}p\left(r_l|s_l\right)}$$
对于无信道编码的系统，直接进行硬判决得到信息比特：
$$\widehat{s_l^k}=\frac{1-sign(L(s_l^k))}{2}$$
对于有信道编码的系统，可以进行硬判决，然后将硬判决结果传给译码器；也可以直接将似然比交予软输入译码器进行译码。前者因为丢失信息，往往有2dB左右的性能损失。

MATLAB中BPSK与QPSK的解调器

BPSK解调器

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

b_hat = r; % 接收到的信道输出
b_hat(b_hat>-0) = 0;
b_hat(b_hat<0) = 1;

QPSK解调器

% 硬判决
qpskDemodHard = comm.QPSKDemodulator('BitOutput', true,...
                        'PhaseOffset', pi/4,'SymbolMapping','Gray');
b_hat = qpskDemodHard(r);

% 软判决
qpskDemodSoft = comm.QPSKDemodulator('BitOutput', true,...
                        'PhaseOffset', pi/4,'SymbolMapping','Gray',...
                        'DecisionMethod', 'Log-likelihood ratio',...
                        'VarianceSource','Property',...
                        'Variance', noiseVar);
Lb = qpskDemodSoft(r);
% passed to soft-in decoder ...

仿真结果

此处就不画仿真结果了。QPSK与BPSK具有相同的误比特率：
$$P_b=\frac{1}{2}erfc(\sqrt{\frac{E_b}{N_o}})$$

补充

若$X\sim \mathcal{N}(0,1)$，则其累积分布函数为:
$$\Phi \left(x\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^x{e}^{\frac{t^2}{2}}dt$$
$Q$函数定义为标准正态分布的右尾函数：
$$Q(x)=\frac{1}{\sqrt{2\pi }}{\int }_x^{\infty }{e}^{-\frac{t^2}{2}}dt=1-\Phi (x)$$
误差函数 $erf(x)$ 的物理意义是服从均值为$0$, 方差为$\frac{1}{2}$正态分布的随机变量$Y$落在区间$(-x,x)$的概率：
$$erf(x)=\frac{2}{\sqrt{\pi }}{\int }_0^x{e}^{t^2}dt$$
互补误差函数$erfc(x)$的物理意义即$Y$落在上述区间外的概率：
$$erfc(x)=\frac{2}{\sqrt{\pi }}{\int }_x^{\infty }{e}^{t^2}dt$$
显然：
$$erfc(x)=1-erfc(x)$$

$$Q(x)=\frac{1}{2}erfc(\frac{x}{\sqrt{2}})$$

附录

假设匹配滤波器具有幅频响应$H(\omega )$，输入$r_i(t)=s_i(t)+n_i(t)$，输出为$r_o(t)=s_o(t)+n_o(t)$。其中：
$$s_o(t)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }S(\omega )H(\omega )e^{j\omega t}d\omega$$
输出信号$s_o(t)$的功率为：
$$P_s=|s_o(t){|}^2=|\frac{1}{2\pi }{\int }_{-\infty }^{\infty }S(\omega )H(\omega )e^{j\omega t}d\omega {|}^2$$
输出噪声的功率为：
$$P_n=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\frac{N_o}{2}|H(\omega ){|}^{2}d\omega$$
抽样时刻$t_0$时的输出信噪比为：
$$\begin{array}{rl}\frac{P_s}{P_n}& =\frac{|\frac{1}{2\pi }{\int }_{-\infty }^{\infty }S(\omega )H(\omega )e^{j\omega t_0}d\omega {|}^2}{\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\frac{N_o}{2}|H(\omega ){|}^{2}d\omega }\\ & \le \frac{\frac{1}{2\pi }{\int }_{-\infty }^{\infty }|H(\omega ){|}^{2}d\omega \cdot \frac{1}{2\pi }{\int }_{-\infty }^{\infty }|S(\omega ){|}^{2}d\omega }{\frac{N_o}{2}\cdot \frac{1}{2\pi }{\int }_{-\infty }^{\infty }|H(\omega ){|}^{2}d\omega }\\ & =\frac{E_s}{\frac{N_o}{2}}\end{array}$$
等号成立的条件是：
$$H(\omega )=K{S}^{\ast }(\omega )e^{-j\omega t_0}$$
显然可以看出实际输出的信号功率和噪声功率与下面的式子有关：
$$\frac{1}{2\pi }{\int }_{-\infty }^{\infty }|H(\omega ){|}^{2}d\omega$$
不过一方面，我们完全可以通过调整$K$，来使上述积分为1；另一方面，它同时乘在输出信号功率和噪声功率上，所以可以不考虑它的影响，而直接认为采样时刻有：
$$\begin{array}{rl}{P}_s& =E_s\\ P_n& ={\sigma }_n^2=\frac{N_o}{2}\end{array}$$
不过注意上述推导是在一个维度下进行的，或者简单认为就是2PAM的情况。
对于两个维度，即IQ调制时，
$$P_n={\sigma }_n^2={\sigma }_{nI}^2+{\sigma }_{nQ}^2=N_o$$