dsfds' blog

学海无涯，回头是岸

0%

MPSK 和 MQAM 的误码率仿真

发表于 更新于 分类于 Valine：

拖延几日，不觉课程已结课几天了，在此趁着脑中还有些关于课程内容的印象，简单写一篇记录，关于如何使用MATLAB构建一个调制技术的仿真程序。

BPSK

首先来看最简单的二进制仿真，以 BPSK 为例，我们来看看如何搭建一个仿真仿真程序，从中得到些一般的仿真设置方法。

在进行仿真之前，先看看我们期待需要得到的是什么——BPSK的误码率 Pe 随信噪比 SNR 的变化曲线，那么这个图像的横坐标就是信噪比，纵坐标就是误码率。这里面还有一个问题需要提前考虑，老师也在上课的时候强调过让我们去思考，即 SNR 和 $\frac{E_b}{N_0}$ 的关系，我们的结论是 $SNR=2\frac{E_b}{N_0}$，这里面有一个从带通等效到低通的过程，同时把基带波形的因素也进行一般化，具体的我只是有一个模糊的概念，等我真的理解了，有机会再来补上吧。最终我们要画的就是 Pe 随 $\frac{E_b}{N_0}$ 的变化曲线。

这里直接给出理论值计算公式，具体推导并不复杂。
$$
P_e=\frac{1}{2}erfc(\frac{E_b}{N_0})
$$

接下来的问题在于我们如何通过 $\frac{E_b}{N_0}$ 定量描述噪声和信号的大小，在仿真中我们需要实际比较噪声和信号幅度的大小确定最终的判定结果。

先来看一个比较笨的办法，这个办法并不具有一般性，它没有使用到一个特别重要的东西——星座图，后续我会利用星座图给出一个更加简单的方法。考虑如下这个公式：
$$
SNR = \frac{P_s}{P_n}=\frac{E_b/T_s}{N_0\cdot B}\times log_2(M)=\frac{E_b}{N_0}=\frac{a^2}{2\sigma^2}
$$

信号幅度a=1，则立即可以得到噪声功率

$$
\sigma^2=\frac{1}{2E_b/N_0}
$$

BPSK 的基带信号是双极性的，用 -1 和 +1 来表示两种不同极性的信号，前面已经得到噪声的功率，那么信号加上噪声幅度就得到接收信号。这就是关键之处，确定信号和噪声的定量关系，已经噪声该如何加到信号上，明白了这个，其余都很简单。最终我们可以得到如下的代码：

BPSK.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
EbN0dB = 0:1:10;
EbN0 = 10.^(EbN0dB/10);
N = 2*10^6;                     % 发送比特数即码元数
Pb = 0.5*erfc(sqrt(EbN0));      % 理论值
ber = zeros(1, length(EbN0dB)); % 仿真误码率 预分配
a = sign(randn(1, N));          % 基带信号
noise = randn(1, N);            % 标准噪声
sigma = sqrt(1./EbN0/2);        % 噪声功率
for n = 1:length(EbN0dB)
    rk = a + sigma(n)*noise;    % 加噪声得到接收信号
    dec_a = sign(rk);           % 判决
    ber(n) = sum(abs(a-dec_a)/2)/length(a);
end

semilogy(EbN0dB, Pb, 'LineWidth', 1);
hold on;
semilogy(EbN0dB, ber, 'rd');
legend('理论结果', '仿真结果')
xlabel('\itE_b/N_0\rm(dB)'); ylabel('Pb');
title('BPSK误码率仿真')
grid on

得到的结果如下所示：

第二种方法更具有一般性，理解了它，其他调制方式也是一样的可以实现。这需要星座图的帮助，通过星座图，我们可以得到最小欧式距离和平均比特能量的关系。在 BPSK 中，有
$$
d_{min}=\sqrt{4E_b}
$$

同时，$\frac{d_{min}}{2}$ 也正是信号的幅度，令 $d_{min}=2$，则信号的幅度为 $a=1$，由 $\frac{E_b}{N_0}$ 和 $d_{min}$ 与 $E_b$ 的关系，我们可以得到
$$
N_0=\frac{E_b}{E_b/N_0}=\frac{d_{min}^2}{4}\cdot\frac{1}{E_b/N_0}
$$

那么噪声功率就是$\sigma^2=\frac{N_0}{2}$。然后按照第一种方法的思路去编写仿真程序即可，非常简单。再来说说它的通用性，使用星座图仿真 AWGN 信道下的调制方式，我们首先需要确定最小欧式距离和平均比特能量的关系，然后得到噪声功率即可，这样就可以得到信号和噪声的定量关系，进行仿真即可。

QPSK

利用上面提到的具有一般性的方式，我们需要使用星座图，QPSK有如下两种星座图形式：

QPSKConstellation

以第一种左边的星座图为例，说明仿真过程。对于 MPSK 调制，最小欧式距离为：
$$
d_{min}=2\sqrt{(log_2M\times sin^2\frac{\pi}{M})\varepsilon_b}
$$
在 QPSK 中，M=4，则 $d_{min}=2\sqrt{\varepsilon_b}$。注意到这和之前的 BPSK 不一样，我们可以将 QPSK 看作两路正交的相位调制叠加。仿真中我们令 $d_{min}=2$，则$\varepsilon =1$，进而可以得到两路信号的幅度都是1，同时还可以确定噪声功率为 $\sigma^2 = \frac{1}{2SNR}$。按照上面的仿真思路我们很容易可以写出代码，最终得到的图形如下所示。

16QAM 和 16PSK

同理，对16QAM和16PSK也是一样的道理，只不过星座点映射在代码实现上会稍稍复杂一点。分析 16QAM 可以将它看作两路正交的 4PAM 信号的叠加，其余分析方法和上面完全相同。16PSK 就更加简单了，不多说，直接给出代码和最终的波形图，如下所示。

QAMandPSK

QAMandPSK.mview raw
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
clear all;
close all;
clc;

EbN0dB = 0:1:15;
EbN0 = 10.^(EbN0dB/10);

M = 16; % 16 QAM
N = 64; % 64 QAM
%% 16QAM 64QAM theory
Pe_16_help = 2*(1-1/sqrt(M))*qfunc(sqrt(3*log2(M)*EbN0/(M-1)));
Pe_16 = 1- (1 - Pe_16_help).^2;
Pb_16 = Pe_16 / 4;
Pe_64_help = 2*(1-1/sqrt(N))*qfunc(sqrt(3*log2(N)*EbN0/(N-1)));
Pe_64 = 1- (1 - Pe_64_help).^2;
%% 16PSK theory
Pe_theory_PSK = 2*qfunc(sin(pi/M)*sqrt(2*log2(M)*EbN0));
%% theory plot
semilogy(EbN0dB, Pe_16, 'LineWidth', 1)
hold on
% semilogy(EbN0dB, Pe_64)
% hold on
semilogy(EbN0dB, Pe_theory_PSK, 'k', 'LineWidth', 1)
hold on
semilogy(EbN0dB, Pb_16, 'c', 'LineWidth', 1)
hold on
semilogy(EbN0dB, Pe_theory_PSK / 4, 'g', 'LineWidth', 1)
grid on

xlabel('\it{E_b/N_0} \rm(dB)')
ylabel('P_b and P_e')

% 16QAM 仿真
% 格雷码映射表
gray_code_table = [0 0 0 0; % 0
                  1 0 0 0;  % 1
                  1 1 0 0;  % 2
                  0 1 0 0;  % 3
                  0 1 1 0;  % 4
                  1 1 1 0;  % 5
                  1 0 1 0;  % 6
                  0 0 1 0;  % 7
                  0 0 1 1;  % 8
                  1 0 1 1;  % 9
                  1 1 1 1;  % 10
                  0 1 1 1;  % 11
                  0 1 0 1;  % 12
                  1 1 0 1;  % 13
                  1 0 0 1;  % 14
                  0 0 0 1;];% 15
% 星座点位置表 16QAM
star_table_1 = [-1*(-3:2:3) (-3:2:3) -1*(-3:2:3) (-3:2:3)];
star_table_2 = [ones(1, 4)*3 ones(1, 4) -1*ones(1,4) -3*ones(1,4)];
star_table = [star_table_1' star_table_2'];
star_table_T = star_table';
table_used = reshape(star_table_T, 1, M*2);
% 星座点位置表 16PSK
phase = 0:2*pi/M:(2*pi-2*pi/M);
PSK16StarTable = [cos(phase);
                  sin(phase)];
PSK_table_used = reshape(PSK16StarTable, M*2, 1);
% 基带信号
N = 4*10^4; % bit数
base = (sign(randn(1, N)) + 1) / 2;
mid_modul = reshape(base, [4 N/4]);
mid_modul = mid_modul';
% 映射到星座点上 QAM 和 QPSK
modulSignalPosi = zeros(N/4, 2);
PSK_modulSignalPosi = zeros(2, N/4);
for symbol = 1:N/4
    dude = gray2int(mid_modul(symbol, :)) + 1;
    modulSignalPosi(symbol, :) = star_table(dude, :);
    PSK_modulSignalPosi(:, symbol) = PSK16StarTable(:, dude);
end
% 噪声
sigma = sqrt(4*(M-1)/6/log2(M)./EbN0/2);
PSK_sigma = sqrt(1/log2(M)./EbN0/2);
% 过信道 解调 求信噪比
be = zeros(1, length(EbN0dB));
se = zeros(1, length(EbN0dB));
PSK_be = zeros(1, length(EbN0dB));
PSK_se = zeros(1, length(EbN0dB));
for index = 1:length(EbN0dB)
    % 加噪声
    reciPosi = modulSignalPosi + sigma(index) * randn(N/4, 2);
    
    % 最小距离解调 QAM
    distanceMid = repmat(reciPosi, [1 M]) - repmat(table_used, N/4, 1);
    distanceMid = distanceMid.^2;
    distance = distanceMid(:, 1:2:(M*2-1)) + distanceMid(:, 2:2:(M*2));
    distance = distance';
    % decodePosi 即为解调结果对应gray_code_table的位置
    [~, decodePosi] = min(distance);
    
    % 解调结果映射到基带 QAM
    toBaseGray = gray_code_table(decodePosi, :);
    toBaseSignal = reshape(toBaseGray', 1, N);
    be(index) = sum(abs(base-toBaseSignal));
    se_posi = double(mid_modul ~= toBaseGray);
    se(index) = sum(sign(sum(se_posi, 2)));
    
    % PSK
    PSK_reciPosi = PSK_modulSignalPosi + PSK_sigma(index) * randn(2, N/4);
    
    PSK_distanceMid = repmat(PSK_reciPosi, [M 1]) - repmat(PSK_table_used, [1 N/4]);
    PSK_distanceMid = PSK_distanceMid.^2;
    PSK_distance = PSK_distanceMid(1:2:(M*2-1), :) + PSK_distanceMid(2:2:(M*2), :);
    [~, PSK_decodePosi] = min(PSK_distance);
    
    PSK_toBaseGray = gray_code_table(PSK_decodePosi, :);
    PSK_toBaseSignal = reshape(PSK_toBaseGray', 1, N);
    PSK_be(index) = sum(abs(base-PSK_toBaseSignal));
    PSK_se_posi = double(mid_modul ~= PSK_toBaseGray);
    PSK_se(index) = sum(sign(sum(PSK_se_posi, 2)));
end
semilogy(EbN0dB, se/N*4, 'v', 'MarkerEdgeColor', 'r')
hold on
semilogy(EbN0dB, be/N, '*', 'MarkerEdgeColor', 'r')
hold on
semilogy(EbN0dB, PSK_se/N*4, '^', 'MarkerEdgeColor', 'm')
hold on
semilogy(EbN0dB, PSK_be/N, '*', 'MarkerEdgeColor', 'm')
legend('16QAM理论误码率', '16PSK理论误码率', '16QAM理论误比特率',  '16PSK理论误比特率', ...
    '16QAM仿真误码率', '16QAM仿真误比特率', '16PSK仿真误码率', '16PSK仿真误比特率')
title('16QAM and 16PSK')