MATLAB复现QPSK/16QAM等调制信道容量曲线

简介：面向数字通信、信息论课程学习者及需要复现仿真实验的读者，这份文档围绕QPSK、8PSK、16PSK与16QAM四种二维M元调制，推导高斯信道下的信道容量曲线，并给出可对照的MATLAB实现思路。内容从二维信号矢量、复高斯噪声模型与信道容量公式出发，结合《Channel Codes: Classical and Modern》中的推导，说明信号能量、噪声方差、信噪比与Eb/N0之间的换算，并展示不同调制星座图的构建与容量曲线绘制过程。包内共1个doc文件，约345KB，集中收录推导过程与源代码，便于按公式逐步核对、理解采样与作图环节。已有1353人学习，适合希望把香农容量理论与MATLAB仿真结合起来，完成课程作业、实验报告或通信系统入门分析的读者参考。

1. 曲线对不上：先搞清楚这份 MATLAB 资源在算什么

复现调制方式的信道容量曲线，第一版图十有八九是错的：QPSK 那条线在 Eb/N0 高段没停在 2 bit，16PSK 和 16QAM 的交叉位置也对不上教材。这份资源给的正是四种调制（QPSK、8PSK、16PSK、16QAM）在高斯信道下的信道容量曲线 MATLAB 实现，附带的推导说明把从互信息定义到可运行脚本这一段补全了，参考文献指向经典信道编码教材。星座点怎么摆、$E_s$ 和 $N_0$ 怎么定、$p(y)$ 这个 M 项高斯混合怎么用样本平均替代积分、横轴怎么从 SNR 换成 Eb/N0，四件事它都碰了。适合做通信方向课程设计、仿真作业，或者想给链路预算找一个理论上界的人。

2. 从香农公式到复高斯混合：容量表达式里每个量怎么定

2.1 直接把 log2(1+SNR) 套上去，四条曲线会重合

第一次复现最容易踩的坑，是把 $C=\log_2(1+\text{SNR})$ 抄进循环，对四种调制各画一条，结果四条线严丝合缝地叠在一起。这不是代码写错，是模型用错——公式里压根没有 $M$。$\log_2(1+\text{SNR})$ 描述的是输入服从复高斯分布、只受平均功率约束时的容量，它在所有输入分布里取最大，是上界，不是某个具体星座的容量。

一旦把输入限制在 $M$ 个离散复点上，问题就变成"离散输入、连续输出"的互信息。回到定义：

$$I(X;Y)=h(Y)-h(Y\mid X)$$

因为 $Y=X+N$，给定 $X$ 之后 $Y$ 的随机性完全来自噪声，所以 $h(Y\mid X)=h(N)$。噪声是复高斯，总方差 $N_0$，微分熵有闭式 $h(N)=\log_2(\pi e N_0)$。麻烦的是 $h(Y)$：它是一个 $M$ 项高斯混合的熵，没有解析解，只能数值算。

这张表的第三列不是摆设。后面所有"曲线整体平移""高 SNR 段翘到天上""低 SNR 出现 NaN"的问题，八成能在某一行的口径上找到原因。

2.2 p(y) 写成 M 项混合，噪声方差的口径必须自洽

$Y$ 的边缘分布是 $M$ 个等概复高斯的平均：

$$p(y)=\frac{1}{M}\sum_{m=1}^{M}\frac{1}{\pi N_0}\exp!\left(-\frac{|y-s_m|^2}{N_0}\right)$$

这里 $|y-s_m|^2$ 是复平面上的欧氏距离平方，$N_0$ 是复噪声总方差。对应到 MATLAB，噪声必须写成

n = sqrt(No/2) * (randn(1,L) + 1j*randn(1,L)); % 每维方差 No/2，复方差 No

实部和虚部各拿一半方差，合起来才等于 $N_0$。很多人的第一版写成sqrt(No)*(randn+1j*randn)，每维方差变成 $N_0$，复方差是 $2N_0$，但概率密度函数里还写着 $1/(\pi N_0)$，两处口径打架，整条曲线会整体右移大约 3 dB，看起来"形状对、位置不对"。

这个错误很难靠眼睛发现，因为曲线形状确实没问题。判断办法是单独估一次噪声功率：

L = 1e6; n = sqrt(No/2)*(randn(1,L) + 1j*randn(1,L)); fprintf('E|n|^2 = %.6f, 期望 N0 = %.6f\n', mean(abs(n).^2), No);

估出来的均值应该贴着 $N_0$，误差在 $1/\sqrt{L}$ 量级。这一步几十毫秒，但能省掉后面半天。

2.3 h(Y) 用样本平均替代积分，采样分布必须是 p(y)

$h(Y)=-\mathbb{E}[\log_2 p(Y)]$ 是个积分，没有闭式。蒙特卡洛的做法是：从真实分布 $p(y)$ 里抽 $L$ 个样本 $y_1,\dots,y_L$，然后用 $-\frac{1}{L}\sum_l \log_2 p(y_l)$ 估计。关键在于"从真实分布抽"——必须先等概随机选星座点，再叠噪声，不能直接从一个高斯里抽。

% 正确：先随机选符号，再加噪声，样本服从 p(y) idx = randi(M, 1, L); x = sig(idx); y = x + sqrt(No/2)*(randn(1,L) + 1j*randn(1,L)); % 错误示范：直接把星座点当噪声均值，忽略了 p(y) 是混合分布 % y = sig(1) + sqrt(No/2)*(randn(1,L) + 1j*randn(1,L));

第二段注释那种写法算出来的是单点高斯的熵，不是混合熵，容量会明显偏大，而且随着 $M$ 增大偏差越来越离谱。

最后合成容量：

Capacity = -mean(log2py) - log2(pi*exp(1)*No);

有个很好用的自检：把 $M$ 设成 1，也就是只发一个符号、不携带任何信息。此时 $p(y)=p(y\mid s_1)$ 是单高斯，代入上式：

$$-\mathbb{E}[\log_2 p(Y)]=\log_2(\pi N_0)+\log_2 e$$

$$C=\log_2(\pi N_0)+\log_2 e-\log_2(\pi e N_0)=0$$

结果精确为零。任何 $N_0$ 口径写错、$1/(\pi N_0)$ 归一化因子漏乘、pi*exp(1)位置放错的问题，都会让 $M=1$ 时算出非零值。这个自检比盯着曲线看十分钟管用，我后来每次改完公式都先跑这一版。

3. 星座点与能量归一：四种调制的 signal_space 怎么写

3.1 手写星座 vs 向量化生成

原资源里四种星座都是逐点手写的，QPSK 写 4 行、8PSK 写 8 行、16PSK 写 16 行，16QAM 又是 16 行。手写的好处是直观、方便对照星座图，坏处是改一个点容易漏、错一眼看不出来，而且 PSK 那部分用的是 $\cos/\sin$ 的正负对称写法，complex(cos(7*pi/8), sin(-7*pi/8))这种表达式读起来费劲——虽然结果是对的（第二象限角度 202.5°，实部为负、虚部为负），但审核成本高。

用等角公式生成更省事，也更难写错：

% M-PSK：单位圆上等角分布，Es = 1 M = 16; m = 0:M-1; sig_psk = exp(1j*2*pi*m/M); % 16QAM：方形星座，实部虚部各取 {-3,-1,1,3} aI = [-3 -1 1 3]; [A,B] = meshgrid(aI, aI); sig_qam = A(:) + 1j*B(:);

meshgrid出来的A(:)是列优先展开，每个实部会和四个虚部配对，正好 16 个点。注意这个顺序和原资源手写的顺序不同，但容量只和星座点集合有关、和编号无关，所以不影响结果；如果后面要接格雷映射或者误码率统计，编号就必须重新排。

四种调制的未归一化平均能量差别不小，这是很多人疑惑"为什么曲线之间能直接比"的原因：

$E_s$ 不同并不影响最终曲线，因为每个 $M$ 都是用自己的 $E_s$ 反推自己的 $N_0$，横轴又换成了 Eb/N0，能量尺度已经被约掉。但如果要叠加一条高斯输入上界 $\log_2(1+E_s/N_0)$ 做对比，就必须先统一归一化，否则上界会整体偏几个 dB，看着像"仿真结果超过理论上界"。

统一归一化很简单，一行：

sig = sig / sqrt(mean(abs(sig).^2)); % 归一化到 Es = 1 Es = 1;

我一般习惯在生成星座之后立刻归一化，后面 $N_0 = E_s/(P_s/P_n)$ 就恒等于 $1/(P_s/P_n)$，少一个变量、少一处口径。

3.2 Es、N0、噪声方差这条链子要闭合

整条链路是这样的：给定 SNR（单位 dB，含义是 $E_s/N_0$）→ 换算成功率比 → 由星座能量反推噪声能量 → 采样噪声 → 估计容量。

snr_dB = -10:1:35; PsPn = 10.^(0.1*snr_dB); % 功率比 = 能量比 = Es/N0 Es = mean(abs(sig).^2); % 平均符号能量 for i = 1:length(snr_dB) N0 = Es / PsPn(i); % 噪声复方差 ... end

N0 = Es/PsPn(i)这一步的前提是"功率比等于能量比"。符号等概发送、噪声零均值，所以平均功率就是平均符号能量除以符号周期，两者只差一个常数，比值完全相同。这也是为什么整个流程里不需要出现符号率、带宽、采样率这些量——容量曲线的横轴是归一化的。

提示：如果你的 SNR 定义是"每比特信噪比"而不是"每符号信噪比"，这里要换成N0 = Es./(log2(M)*EbN0_linear)，否则横轴会差一个 $\log_2 M$ 的因子，8PSK 和 16QAM 的曲线会整体平移 1.5 dB 和 2 dB，而 QPSK 恰好不动，很容易被误判成"只有 QPSK 是对的"。

3.3 把双重 for 循环换成矩阵运算

原资源的写法是对每个 SNR 点、每个样本点、每个星座点三层循环。sample=66666、46 个 SNR 点、$M=16$，内层exp要算将近 5000 万次，MATLAB 里跑几分钟很正常，改一次参数就要重来一遍，调试体验很差。

把最内两层换成矩阵运算，速度快一个数量级，代码也更短：

L = 200000; % 每个 SNR 点的样本数 y = x + sqrt(N0/2)*(randn(1,L) + 1j*randn(1,L)); D = abs(y.' - sig).^2; % L x M 距离平方矩阵 py = mean(exp(-D/N0), 2) / (pi*N0); % 按行求 M 项平均 C = -mean(log2(py)) - log2(pi*exp(1)*N0);

y.'是 L×1 列向量，sig是 1×M 行向量，相减自动广播成 L×M，一次就把所有样本对所有星座点的距离算完。mean(...,2)沿第二维求平均，得到每个样本的 $p(y_l)$。三步之后容量就出来了。

这里还有一个数值细节：exp(-D/N0)在 $N_0$ 很小、距离又大的时候会下溢成 0，log2(0)得到-Inf，容量直接变成Inf。高 SNR 段尤其容易撞上。稳妥的写法是先做 log-sum-exp：

D = abs(y.' - sig).^2; mx = min(D, [], 2); % 每个样本的最近距离 LSE = mx + N0*log(sum(exp(-(D - mx)/N0), 2)); % 数值稳定的 log(sum(exp)) log2py = (LSE - N0*log(M) - log(pi*N0)) / log(2); C = -mean(log2py) - log2(pi*exp(1)*N0);

减去mx之后再取指数，最大项恒为 1，不会全下溢。代价是多两次矩阵运算，但换来的是高 SNR 段不会突然冒出一根竖线。

4. 横轴换成 Eb/N0：一次坐标变换里的三个口径

4.1 为什么横轴要除以 Capacity

容量 $C$ 的单位是 bit 每信道使用（每复维）。$E_s$ 是每次信道使用的平均能量，$E_b$ 是每比特能量，两者关系是 $E_s = C\cdot E_b$——这个等式只在"以容量速率传输"时成立，因为只有这时候每符号才恰好携带 $C$ 个比特。

$$\frac{E_b}{N_0}=\frac{E_s/N_0}{C}$$

对应代码就是EbNo = 10*log10(PsPn./Capacity)。这一步把不同 $M$ 的曲线拉到同一个可比坐标系下：QPSK 在 $C\to2$ 时 $E_b/N_0\to\infty$，因为它最多只能传 2 bit；16QAM 在 $C=2$ 处仍然是有限值，因为 4 bit 的星座传 2 bit 绰绰有余。

这也是我判断曲线对不对的第一眼依据：QPSK 曲线在 $C$ 接近 2 的位置应该向右发散，而 16QAM 曲线会继续向左延伸，两条线在 $C$ 略小于 2 的地方交叉。如果画出来 QPSK 一直贴在 16QAM 左边不交叉，那基本可以确定 $M=1$ 自检没做、或者噪声方差口径错了。

4.2 低 SNR 段 Capacity 接近 0，除法会炸

PsPn./Capacity在 $C\to0$ 时趋向无穷，而蒙特卡洛估出来的 $C$ 可能是个很小的正数、甚至因为样本波动变成负数。负数取log10得到复数，plot会警告；极小的正数会得到 $E_b/N_0=60$ dB 这种离谱值，把横轴撑到看不清。

处理办法是给容量加个下限，或者干脆按容量裁剪绘制区间：

Cplot = Capacity; Cplot(Cplot < 1e-3) = NaN; % 太小的点不画，避免除零 EbNo = 10*log10(PsPn./Cplot);

用NaN而不是直接剔除，好处是plot会自动断开线段，不会把不相邻的两点连起来画出一条假的斜线。如果一定要保留低 SNR 段，那就把样本数 $L$ 加到百万级，让估计方差压下去——低 SNR 下 $C$ 本身就小，相对误差被放大得厉害，样本不够时曲线抖得像心电图。

4.3 出图：用 legend，别用 text 标位置

原资源用text(18,2.2,'QPSK')这种绝对坐标标注曲线，一旦改变 SNR 范围、样本数或者归一化方式，标注就飞到图外面去了。换成legend一次解决：

figure; hold on; box on; grid on; plot(EbNo_qpsk, C_qpsk, 'b--', 'LineWidth', 1.2); plot(EbNo_8psk, C_8psk, 'r-.', 'LineWidth', 1.2); plot(EbNo_16psk, C_16psk, 'k:', 'LineWidth', 1.2); plot(EbNo_16qam, C_16qam, 'c-', 'LineWidth', 1.2); % 高斯输入上界，作为参照 snr_lin = 10.^(snr_dB/10); Cg = log2(1 + snr_lin); plot(10*log10(snr_lin./Cg), Cg, 'g-', 'LineWidth', 1); xlabel('E_b/N_0 (dB)'); ylabel('Capacity (bit/channel use)'); axis([-2 20 0 4.5]); legend('QPSK','8PSK','16PSK','16QAM','Gaussian input','Location','southeast');

纵轴上限给到 4.5 而不是 5，是因为 16QAM 饱和在 4，留半格刚好；横轴从 -2 起，是因为高斯输入的香农极限在 -1.59 dB，低于这个值的区域所有曲线都应该发散，画出来只会是一堆噪声。上界那条线最好一起画——它是判断"仿真结果有没有超过理论上界"的唯一参照，一旦某条曲线跑到绿线右侧还高出一截，说明估计过程有问题，不用再看别的。

5. 避坑与排查：五条从跑不出图到跑对的记录

5.1 曲线形状类问题

现象：高 SNR 段（比如 Eb/N0 > 15 dB）容量突然冲到 10 以上，或者图上冒出一根竖直的线。

原因：exp(-D/N0)在 $N_0$ 很小时整列下溢成 0，log2(0)得到-Inf，-mean(...)变成+Inf。或者某个样本的 $p(y)$ 极小，取对数后把均值拉爆。

解决：按 3.3 的 log-sum-exp 写法重写，或者至少加一句log2py(~isfinite(log2py)) = [];把异常点剔掉再求均值。后者治标，前者治本。

现象：四条曲线在低 SNR 段完全重合，高 SNR 段才分开。

原因：这不是 bug，是物理规律。$E_s/N_0$ 很小时所有星座点在噪声尺度下几乎重合，任何 $M$ 元星座的容量都逼近 $\log_2(e)\cdot E_s/N_0$，斜率相同。如果连高 SNR 段都重合，那才是真出问题了——回去检查sig有没有在循环外被覆盖。

5.2 数值口径类问题

现象：整张图相对教材右移约 3 dB。

原因：噪声采样写成sqrt(No)*(randn+1j*randn)，每维方差 $N_0$，复方差 $2N_0$，但概率密度里还按 $N_0$ 算。等效于实际噪声比公式假设的大一倍。

解决：改成sqrt(No/2)，然后跑一次 2.2 里的噪声功率估计确认 $E|n|^2=N_0$。这个检查我后来加进了固定流程，改任何跟噪声相关的行之后都跑一遍。

现象：$M=16$ 的曲线容量超过 4，或者 $M=8$ 的超过 3。

原因：混合分布的项数没对齐。检查sig的长度是不是真的等于 $M$——用meshgrid生成 16QAM 时如果写成[A,B]=meshgrid(aI)只传一个参数，出来是 4×4 没问题；但如果中途sig被别的分支覆盖成 8 个点，而M还写着 16，1/M就乘错了，容量会偏大 $\log_2(16/8)=1$ bit。

解决：在循环开头加一句assert(numel(sig)==M)，两秒钟的事，能挡住一整类错误。

5.3 变量与索引类问题

现象：16PSK 和 16QAM 两条曲线完全重合，或者 8PSK 那条画出来是 16 个点的形状。

原因：signal_space在循环之间没有清空。原资源里每个 $k$ 分支都完整重写了一遍数组，所以碰巧没问题；但只要改动其中一个分支的长度，短数组残留的尾部元素会被sum(...)/M算进去，$E_s$ 就错了。

解决：每个分支开头sig = zeros(1,M);或者干脆用函数封装，让每次调用都是新变量。我现在的习惯是把星座生成写成一个独立函数gen_constellation(M, type)，主循环只调它，不共享工作区。

现象：x(j)=signal_space(floor(rand*M)+1)偶尔报"索引超出数组范围"。

原因：rand理论上取不到 1，但如果换成round(rand*M)+1，rand接近 0 时round得 0，索引变成 1 没问题；接近 1 时round(M)+1 = M+1，越界。用floor是对的，用round是错的。

解决：统一用randi(M, 1, L)，语义清晰也不会越界，还能顺手指定随机种子做可重复实验：rng(2024)放在脚本最前面。

6. 交叉验证：三个基准点判断曲线到底对不对

曲线画出来只是第一步，能不能判断它对不对才是分水岭。我习惯用三个不依赖参考图的硬基准去卡。

第一个是饱和值。任何 $M$ 元星座的容量上限是 $\log_2 M$，QPSK 到 2、8PSK 到 3、16PSK 和 16QAM 到 4。看 $E_b/N_0=20$ dB 处四条线的纵坐标，只要有一个明显超过，就是 5.2 里说的 $M$ 与数组长度不匹配。反过来，如果 16QAM 在 20 dB 处还没摸到 3.9，多半是样本数不够、低概率区域没采到。

第二个是低 SNR 斜率。信噪比很小时，任何星座的容量都趋近于 $\log_2(e)\cdot E_s/N_0 \approx 1.4427\cdot E_s/N_0$。取 $-10$ dB 这个点，线性 SNR 是 0.1，容量应该约等于 0.144，四种调制几乎一样。这个数可以直接和图上读数对，偏差超过 10% 说明采样或 $N_0$ 口径有问题。

第三个是香农极限。高斯输入下 $\min E_b/N_0 = \ln 2 = -1.59$ dB，所有有限星座的曲线都必须在这条竖线的右边。如果哪条线跨过了 $-1.59$ dB，不用查代码，一定是除法的分母用错了。

这张表我贴在了脚本开头的注释里，改一次代码跑一遍，比截图比对快得多。有一回我把N0从 $E_s/(P_s/P_n)$ 改成 $E_s/(\log_2 M \cdot P_s/P_n)$，想当然觉得"每比特能量应该除以 $\log_2 M$"，结果四个基准全对、曲线整体平移 2 dB，肉眼完全看不出来——最后是 $-10$ dB 那个点从 0.144 变成 0.036 才暴露的。从那以后我每次动到信噪比定义，都强制把这三行基准数打成fprintf打到命令行，不看一眼不往下走。希望这套流程能帮到你。

本文还有配套的精品资源，点击获取