简介：这是一份面向通信工程专业学生与MATLAB仿真学习者的技术文档，系统整理了QPSK、8PSK、16PSK和16QAM四种调制方式在高斯白噪声信道下的信道容量曲线推导与计算过程。文档以《Channel Codes: Classical and Modern》为参考，从二维M元信号集合的复数表示和星座图出发，结合复高斯噪声概率密度函数，逐步推导出M进制调制容量公式，并详细解释了容量公式中期望项的计算方法。随后给出了从信噪比到Eb/N0的转换思路，以及在MATLAB中通过遍历-10~35 dB信噪比范围、使用大样本采样来估计容量值的完整程序段，包含四种调制方式的星座图设置，可直接运行或二次修改。包内共1个doc文件，大小约345KB，既有理论推导也有可复现的源代码，适合用于通信原理课程设计、毕业设计或调制理论自学。目前已有1352人学习下载，是一份兼顾理论深度与工程实践的实用参考资料。
1. 不只有香农公式：这份资源把 QPSK、8PSK、16PSK、16QAM 的信道容量曲线一次算清
信道容量曲线是通信系统设计里绕不开的参照系：给定信噪比，一个调制方式到底能撑起多少 bit/s/Hz 的可靠传输，直接决定编码方案和调制阶数的选择。很多人上手就写C = 0.5*log2(1+SNR)，但那是连续高斯输入下的极限值，换成 QPSK、8PSK、16PSK、16QAM 这些离散星座，容量公式要基于互信息重新推一遍。这份课程设计资源把推导过程和 MATLAB 实现都打包了，核心是用蒙特卡洛采样估计二维复高斯信道下的互信息，然后画出容量与 Eb/N0 的关系曲线。适合正在做数字调制课程设计的学生、需要估算系统编码余量的通信工程师，以及任何想把"星座图长什么样"和"容量上限是多少"对应起来的人。
2. 高斯信道下 M 元调制的容量：从连续互信息到可计算的蒙特卡洛形式
2.1 为什么离散星座不能直接套香农公式
香农容量公式给出的 AWGN 信道容量，是假设输入信号服从连续高斯分布时取得的理论上限。实际系统里输入信号取自有限集合，比如 QPSK 只有四个复平面上的点，携带的信息量受限于星座点数，也受限于这些点之间的距离。在高信噪比下，M 元调制的容量上限趋近 log2(M) bit/符号，而不是像高斯输入那样随 SNR 无界增长。这就是为什么需要重新推导针对离散输入集合的容量表达式。
二维 M 元信号集合 {s_m} 中的每个符号都可以用两个归一化正交函数的线性组合表示，在基带等价于一个复数。高斯信道中二维 M 进制调制的容量公式，可以看成二进制情况的直接推广。核心关系是互信息：
I(X;Y) = h(Y) - h(Z)其中 h(Y) 是接收信号 y 的微分熵，h(Z) 是复高斯噪声的微分熵。因为噪声是二维复高斯，方差为 N0，其微分熵是固定的：
h(Z) = log2(pi * e * N0)所以容量计算的关键落到了 h(Y) 上。接收信号 y 的概率密度不是简单的高斯分布，而是 M 个高斯分量的混合，每个分量对应一个星座点 s_m。这个混合密度直接决定了互信息的大小。
2.2 容量公式的分步推导：p(y) 的混合高斯结构
在等概率发送的前提下，接收信号的概率密度是 M 个复高斯密度的平均：
p(y) = (1/M) * sum_{m=1}^M p(y|s_m) p(y|s_m) = (1/(pi*N0)) * exp(-|y - s_m|^2 / N0)噪声每一维的方差是 N0/2，所以复噪声总的方差是 N0，这解释了密度函数中系数是 pi*N0 而不是 sqrt 形式。把 h(Y) 写成期望形式：
h(Y) = -E[log2 p(y)]这个期望没有闭式解，但可以用蒙特卡洛估计。生成 L 个带噪样本 y_l，每个样本都独立地从一个随机选中的星座点加复高斯噪声得到，然后用样本均值近似期望：
C ≈ - (1/L) * sum_{l=1}^L log2( (1/M) * sum_{m=1}^M (1/(pi*N0)) * exp(-|y_l - s_m|^2 / N0) ) - log2(pi * e * N0)这个形式就是代码里log2py和Capacity两行表达式的来源。公式（4）里的第一项是 log2(p(y)) 的均值，用样本平均计算；第二项是噪声熵，直接带入 pieN0。
2.3 蒙特卡洛估计的三个关键参数：样本数、信噪比范围和能量归一化
蒙特卡洛估计的精度取决于样本数 L。代码里sample=66666是一个偏大的取值，理论上每个信噪比点都有六万多个随机样本，均值估计的方差足够小，曲线会很光滑。但代价是计算量：每个样本都要遍历 M 个星座点计算 exp，四层循环（信噪比 × 样本 × 星座点 × 四种调制方式）跑起来很慢，这一点在后面的避坑章节展开。
信噪比范围设为snr=-10:1:35，覆盖了从容量几乎为零的低信噪区到容量接近 log2(M) 的高信噪区。注意这个 snr 是 Es/N0 意义上的信噪比，不是 Eb/N0。代码里通过PsPn=10.^(0.1*snr)得到功率比，然后用No=Es/PsPn(i)算出每个符号对应的噪声方差。Es 是调制信号的平均能量，由星座点的模值平方均值求得。这个过程成立的前提是功率比等于能量比，因为符号周期相同，功率和能量只差一个常数因子。
能量归一化是这里最容易被忽略的点。四种调制的 Es 必须一致，画出来的曲线才有可比性。代码里 QPSK、8PSK、16PSK 的星座点都落在单位圆上，Es=1；但 16QAM 用的是坐标 1 和 3 的网格，Es=(42+810+4*18)/16=10。这导致 16QAM 曲线和三种 PSK 曲线不在同一个能量基准上比较，后面会专门讨论这个问题。
3. MATLAB 代码逐段拆解：星座初始化、噪声叠加与三层循环
3.1 信号空间与四种星座的定义方式
代码先定义信噪比范围和样本数，然后按调制阶数依次生成星座点。QPSK 的四个点直接落在复平面的 (±1, ±1) 上，8PSK 和 16PSK 都在单位圆上按等角度分布。8PSK 的坐标用手写的三角函数值给出，16PSK 也类似，都是角度从 0 开始、每次增加 pi/4 或 pi/8 的圆周点。16QAM 则是 1 和 3 组成的 4×4 网格。
% QPSK 星座：四个点落在 (±1, ±1) signal_space(1) = complex(1, 1); signal_space(2) = complex(-1, 1); signal_space(3) = complex(-1, -1); signal_space(4) = complex(1, -1); % 8PSK 星座：单位圆上 8 等分，Es=1 signal_space(1) = complex(1, 0); signal_space(2) = complex(cos(pi/4), sin(pi/4)); signal_space(3) = complex(0, 1); signal_space(4) = complex(cos(3*pi/4), sin(3*pi/4)); signal_space(5) = complex(-1, 0); signal_space(6) = complex(cos(5*pi/4), sin(5*pi/4)); signal_space(7) = complex(0, -1); signal_space(8) = complex(cos(-pi/4), sin(-pi/4));这段代码的特点是手写而非循环生成。对于 QPSK 和 8PSK，手写可读性好，但 16PSK 的 16 个点手写就很容易写错角度符号。对比 16PSK 代码里第 10 到第 16 个点，角度从 7*pi/8 开始取负值，是沿着单位圆从第一象限顺时针绕回。更稳妥的写法是用exp(j*theta)循环生成，theta 取(0:M-1)*2*pi/M，既减少笔误又方便扩展到 32PSK。16QAM 的坐标也有规律，是 {-3,-1,1,3} 两两组合，用两层 for 循环加一个索引映射表可以取代 16 行手写赋值。
三种 PSK 的星座点都在单位圆上，所以 Es 都是 1，这也让三种 PSK 曲线可以直接用同一个 No 计算公式。16QAM 就不同了，它的 Es 是 10，代码里单独算了一次Es=sum(abs(signal_space).^2)/M，这一步没有错，但埋下了和 PSK 曲线对比时能量基准不一致的问题。
3.2 蒙特卡洛主循环：随机信号生成、复高斯噪声叠加与 pysm 累加
主循环结构分为三层：外层遍历信噪比，中层生成 sample 个输入符号并叠加噪声，内层对每个接收样本计算 log2(p(y)) 的贡献。以下是核心循环代码：
for i=1:1:length(snr) No = Es / PsPn(i); % 噪声方差，由能量比等于功率比得到 for j=1:1:sample x(j) = signal_space(floor(rand*M)+1); % 等概率随机选取星座点 end y = x + sqrt(No./2) * complex(randn(1,sample), randn(1,sample)); for j=1:1:sample pysm = 0; for b=1:1:M pysm = pysm + exp(-(abs(y(j)-signal_space(b))).^2 / No); end log2py(j) = log2(1/M/(pi*No) * pysm); end Capacity(i) = -mean(log2py) - log2(pi*exp(1)*No); endx(j) = signal_space(floor(rand*M)+1)用rand生成 [0,1) 均匀随机数，乘 M 后取整得到 1 到 M 的索引。这个写法对，但每次调用rand都有函数开销，后面会提到换成randi更干净。
噪声叠加行是关键：sqrt(No./2) * complex(randn, randn)。复噪声的实部和虚部各占 N0/2 方差，合起来的复噪声方差是 N0，和公式 (3) 中的密度函数系数一致。这里容易犯的错是直接用sqrt(No) * randn只加实噪声，那相当于把二维问题退化成一维，容量结果会和正确值差一截。
内层pysm累加的是 exp(-|y-s_m|^2/No)，也就是公式里高斯分量的非归一化部分。除以1/M/(pi*No)后取 log2，得到 log2(p(y)) 的单样本估计。最后Capacity = -mean(log2py) - log2(pi*exp(1)*No)，mean直接给出期望的样本均值，减去噪声熵就是互信息。这一段逻辑和公式 (4) 完全对应，没有多余的技巧，但执行效率是最差的——每算一个 y_j 都要循环 M 次 exp，sample 取 66666 时这一层循环要跑约三十万次乘加，全程序跑完需要等待较长时间。
3.3 从 SNR 到 Eb/N0 的坐标变换：容量单位与换算公式
画图时横轴用 Eb/N0 而不是 SNR，因为不同调制阶数每符号携带的比特数不同，只有换算到每比特能量才能公平比较。换算公式是：
Eb/N0 = (Ps/Pn) / CC 的单位是 bit/符号，所以 (Ps/Pn)/C 就是每比特能量与噪声谱密度之比，取 10*log10 后得到 dB 值。代码里写作：
EbNo = 10*log10(PsPn./Capacity);这个式子成立的前提是 Capacity 严格为正。在低信噪比区域，蒙特卡洛估计的容量可能因为样本噪声变成很小的负数，直接除会出现负的 Eb/N0，导致曲线出现不合理的急剧下降。这个问题在后面避坑章有具体处理办法。
画图部分用四种线型区分调制方式，axis([-2,20,0,5])把横轴限制在 -2 到 20 dB，纵轴限制在 0 到 5 bit/符号。QPSK 理论极限是 2 bit/符号，16QAM 是 4 bit/符号，所以纵轴取 5 够用。text函数在 x=18 的位置放置图例文字，避免了 legend 在不同 MATLAB 版本下的兼容性问题。
4. 四种调制的曲线对比与公平性：能量归一化的关键细节
4.1 同一 Es 下三种 PSK 的容量排序与直觉解释
QPSK、8PSK、16PSK 在单位圆上等间隔分布，Es 都等于 1，它们的容量曲线可以直接对比。直观规律是：调制阶数越高，高信噪比下的容量上限越高，但达到上限所需的 Eb/N0 也越高。QPSK 曲线从大约 0 dB 开始快速爬升，在 10 dB 左右就接近 2 bit/符号的渐近线；8PSK 需要更高的 Es/N0 才能拉开和 QPSK 的差距，但上限是 3 bit/符号；16PSK 的 16 个点挤在单位圆上，相邻点距离最小，低信噪比下误码概率更高，容量曲线在三种 PSK 里爬得最慢，但最终趋近 4 bit/符号。
这个排序的物理意义在于：星座点数越多，每个符号携带的信息量越大，但星座点之间的欧氏距离随之减小，抗噪声能力下降。容量曲线反映的是这个折中的定量结果——在某个 Eb/N0 以下，用低阶调制反而能获得更高的可靠速率。
4.2 16QAM 的 Es=10 问题：与 PSK 对比前必须归一化
16QAM 的星座点坐标是 ±1 和 ±3 的组合，平均能量是 10，而三种 PSK 的 Es 是 1。代码里 16QAM 部分重新算了 Es，所以它自己的曲线是自洽的，但和 PSK 曲线放在同一张图里对比时，能量基准不一样。16QAM 的曲线整体会比它在 Es=1 情况下右移 10*log10(10)=10 dB，也就是说，如果不做归一化，16QAM 看起来在低信噪比下容量明显低于 QPSK，但这其实是能量没对齐造成的假象。
常见做法是把 16QAM 的星座点除以 sqrt(10) 再做仿真，让所有调制方式的 Es 都等于 1，这样四条曲线就在同一个能量基准上比较。坐标缩放不影响容量曲线的形状，只影响横轴位置。如果你用这份代码做课程设计或系统级仿真，建议把 16QAM 星座点统一除以 sqrt(Es)，既方便对比，也符合大多数通信教材里星座图归一化的习惯。
4.3 从曲线读取设计结论：编码门限与调制阶数选择
容量曲线的用途不只是展示上限，更重要的是为编码方案选择提供依据。一个实用的读法是：给定目标频谱效率（比如 2 bit/s/Hz），看哪种调制在更低的 Eb/N0 下能达到这个容量，那个调制就是能量效率更优的选择。QPSK 在 2 bit/s/Hz 附近所需的 Eb/N0 明显低于 16QAM 在同一频谱效率下的需求，这意味着用 QPSK 加码率 1/2 的编码，可能比用 16QAM 加高码率编码更容易逼近容量。
另一个读法是看曲线拐点。容量曲线在低信噪比区域缓慢爬升，接近渐近线时明显变平，拐点大致对应"星座点间距与噪声标准差相当"的临界区域。在实际系统设计中，工作点一般选在容量曲线下方 1~2 dB 处，预留编码调制余量。这份资源的曲线可以作为这种预算分析的参照。
5. 复现避坑指南：样本数、负容量、循环提速与曲线截断
5.1 sample=66666 跑一次要等几分钟
现象：程序在for j=1:1:sample和for b=1:1:M两层循环里反复算 exp，46 个信噪比点、四种调制（16QAM 还要单跑一遍），总计约 46×4×66666×16 次 exp 运算，普通笔记本跑一次要十几分钟，信噪比点再多一点就直接失去耐心。
原因：蒙特卡洛估计确实需要足够多样本才能让曲线平滑，但 66666 这个数字没有理论依据。容量估计的均方误差与样本数成反比，样本数到 1 万左右时曲线的抖动已经在可接受范围。
解决：把sample降到 20000，或用parfor替换外层信噪比循环。更彻底的办法是向量化内层循环，一次性算出 sample×M 的指数矩阵，后面章节给出具体代码。改完参数后先用sample=5000验证曲线形状，确认稳定后再放大样本数画正式图。
5.2 低信噪比下容量出现负值
现象：snr=-10 dB 附近，Capacity(i)算出负数，有时候小到 -0.3，画图时纵轴出现负值，曲线穿到横轴下方。
原因：蒙特卡洛估计本身有方差。真实容量在极低信噪比下趋近于 0，理论上不可能为负，但样本均值估计 h(Y) 时，某些样本的 log2(p(y)) 波动会把均值推到噪声熵以下，产生轻微负数。sample 越小、信噪比越低，这种现象越明显。
解决：把负容量截断为 0，Capacity(i)=max(0, Capacity(i))。这个处理不影响有效区域的曲线形状，因为负值只出现在容量接近 0 的低信噪区。另一个补充手段是增大样本数，但截断更简单直接，且符合容量非负的物理约束。
5.3 16QAM 曲线比其他调制低一大截，怎么调都对不上
现象：同样在 15 dB Eb/N0 处，QPSK 容量约 2 bit/符号，16QAM 算出来只有 3 左右，甚至更低，曲线整体右移约 10 dB。
原因：16QAM 的 Es=10，而三种 PSK 的 Es=1。代码里 16QAM 部分单独计算了 Es，所以它的 No=Es/PsPn 和 PSK 路径不统一，等于在比 PSK 高 10 dB 的 Es/N0 点上做仿真，曲线自然右移。
解决：在 16QAM 星座初始化之后、进入信噪比循环之前，加一行归一化：signal_space = signal_space / sqrt(Es);再重新计算Es=1。这样四种调制就在同一能量基准上画图，对比才有意义。这是复现这份代码时最容易踩的坑，强烈建议先固定 Es 再跑仿真。
5.4 从 SNR 转 Eb/N0 时曲线出现剧烈抖动或断点
现象：EbNo=10*log10(PsPn./Capacity)这行执行后，某些点的 Eb/N0 变成很大的正数或负数，曲线出现尖刺。
原因：当 Capacity 接近 0 时，除法结果接近无穷，取对数后数值不稳定。特别是低信噪比区，容量估计的微小波动会被除法放大，产生极端值。另外如果容量为负，log10会得到 NaN，曲线直接断裂。
解决：先做容量截断（第 5.2 条），再做 Eb/N0 转换。转换后对超出axis([-2,20,0,5])范围的数据点可以不处理，因为 MATLAB 画图时会自动裁剪，但要确认没有 NaN 点，否则整条曲线都会消失。
5.5 画图时曲线在边界被切断，图例文字被截
现象：axis([-2,20,0,5])设定后，snr 范围是 -10 到 35，对应 Eb/N0 范围可能超过 20 dB，曲线在 x=20 处被垂直切断，看上去像容量饱和后又突然掉下来。
原因：坐标轴范围写死，而不同调制方式在低信噪比下的 Eb/N0 可以远大于 20。特别在容量趋近 0 的区间，Eb/N0 会飙升到几十甚至上百 dB，超出横轴范围。
解决：把横轴范围放宽到 [-2, 25] 或者干脆不设 x 轴上限，用xlim([-2 25])配合ylim([0 5])让曲线完整显示。图例文字放在 x=18 处也可能和曲线重叠，改成text的位置按各曲线在右上角的实际走势微调，或者直接用legend在 MATLAB 2019 之后的版本里有自动防重叠功能。
6. 让曲线更快更稳：三个小改动与一条验证技巧
第一个改动是向量化内层循环。原来的三层循环可以压缩成两层，用 MATLAB 的矩阵运算一次算出所有样本对所有星座点的距离：
for i=1:1:length(snr) No = Es / PsPn(i); x_idx = randi(M, 1, sample); % 直接生成整数索引 x = signal_space(x_idx); % 抽取随机星座点 y = x + sqrt(No/2) * complex(randn(1,sample), randn(1,sample)); dist2 = abs(bsxfun(@minus, y.', signal_space(:).')).^2; % sample×M 距离矩阵 py = sum(exp(-dist2/No), 2) / M / (pi*No); % 每个样本的 p(y) 估计 log2py = log2(py); Capacity(i) = -mean(log2py) - log2(pi*exp(1)*No); endbsxfun在旧版本 MATLAB 里可用，新版本直接支持y.' - signal_space(:).'的隐式扩展。距离矩阵一次性算出来，避免逐样本逐星座点的 exp 循环，速度提升通常在十倍以上。randi(M,1,sample)比floor(rand*M)+1更简洁，也避免每次调用 rand 的冗余开销。
第二个改动是给容量加截断保护，在Capacity(i)计算后加一行Capacity(i) = max(0, Capacity(i));，防止低信噪比下的负值影响 Eb/N0 转换。这是让曲线在低信噪区看起来干净的唯一必要处理。
第三个改动是固定随机种子。在程序开头加rng(42);，这样每次跑出来的曲线完全一致，方便调参数时对比差异。调优阶段用rng('shuffle')观察不同种子下的波动范围，确认曲线形状的稳健性。
验证技巧是画一张对照图：把 QPSK 的仿真容量曲线和 AWGN 容量公式C = log2(1+snr)画在同一张图上，QPSK 在高信噪比区应该贴着这条理论上限走，因为四个点可以等效近似连续高斯输入。如果 QPSK 曲线明显低于理论线，说明样本数不够或噪声方差算错了；如果曲线在高信噪比区超过 log2(M)，说明容量估计出了问题，需要回头检查星座点坐标。另外可以验证 8PSK 在 Eb/N0=-1.59 dB 处容量不超过 0，因为 -1.59 dB 是整个 AWGN 信道的香农极限，任何离散调制在该点都不可能取得正速率。
我从那以后每次复现容量仿真，都强制先做一遍能量归一化和随机种子固定，再开始调参。这个方法帮我少走了不少弯路，数据也对得上教材曲线。希望帮到你。
本文还有配套的精品资源，点击获取
