简介面向通信原理学习与数字调制仿真实操的 MATLAB 代码包围绕 16QAM 星座图绘制与误码率分析展开适合通信工程、电子信息类专业学生及需要验证调制解调流程的开发者。压缩包共 8 个文件均为 .m 脚本或函数整体约 3KB涵盖主仿真流程、星座图可视化、四进制与二进制码字转换、调制解调、随机二进制序列生成及比特成形等环节可串起从信号生成、AWGN 信道模拟到最佳接收与 BER 统计的端到端实验。已有 1927 人学习说明其在课程实验与自学场景中具有一定参考度。读者可据此复现星座点映射观察噪声下接收点分布理解最大似然检测思路并通过不同信噪比下的误码率曲线评估系统性能适合作为通信系统仿真入门与课程设计的对照材料。1. 星座图和 BER 曲线放在同一张图上问题才藏不住16QAM 星座图仿真和误码率仿真是基带链路入门绕不开的一对组合前者看符号摆得对不对后者看比特错得有多少。很多人第一次做直接在脚本里写scatterplot(qammod(...))加berawgn(...)两条线一叠就算交差结果答辩时被问一句“你这条仿真曲线为什么比理论高 3 dB”就答不上来。真正能用的仿真是把映射表、归一化因子、噪声方差、判决门限这四件事各自闭环星座图只是它们的可视化出口。这套 MATLAB 脚本面向三类人做通信课程设计要交曲线和散点图的学生、刚转基带算法想搞清楚 Eb/N0 与 Es/N0 换算的工程师、以及需要给接收机指标做预算、想用仿真反查 EVM 门槛的人。2. 16QAM 的数学底子格雷映射决定了 BER 曲线能压多低动手写代码之前先把星座点的位置、能量和比特顺序定死。这三样任何一样定错后面画出来的图都只是好看不能用来下结论。这一章不写代码先把“为什么这么摆”说透。2.1 星座点为什么排成方阵I/Q 两路 4-PAM 的笛卡尔积16QAM 最主流的排布是矩形方阵同相支路 I 取 4 个电平 {-3, -1, 1, 3}正交支路 Q 也取同样的 4 个电平两两组合得到 16 个星座点。这样做的直接好处是解调可以拆成两个独立的一维 4-PAM 判决——I 路只看实部Q 路只看虚部每个维度三个门限-2、0、2判决复杂度是线性的硬件实现里就是三组比较器。选矩形而不是圆形是因为同样的平均符号能量下矩形星座的最近邻距离更大。未归一化时矩形 16QAM 的最小距离 d_min 2平均符号能量 Es 10归一化最小距离平方 d_min²/Es 0.4。而 16PSK 相邻点夹角 22.5°d_min²/Es 2(1 − cos22.5°) ≈ 0.152。两者差了约 4.2 dB这就是同样误码率下 16PSK 要多花 4 dB 功率的原因。星形 16QAM 能做到比矩形略好一点但映射和判决都变复杂工程上除了特定场景很少用。所以做仿真、做课程设计、做链路预算默认拿矩形方阵起步是合理的。还有一点容易被忽略矩形星座的符号错误几乎全部来自“错到最近邻”。水平或垂直相邻的四个点距离是 d_min对角相邻的点距离是 √2·d_min在中等以上信噪比时贡献的错误概率低一个量级。这个性质是后面 BER 近似公式成立的前提也是格雷映射能起作用的前提。2.2 归一化因子 0.3162 从哪来为什么必须写死一维 4-PAM 的四个电平是 ±1、±3等概率出现时平均功率是 (1 1 9 9)/4 5。I、Q 两路独立所以 16QAM 的平均符号能量 Es 5 5 10。为了让后续所有公式里的 Es 都等于 1需要乘一个归一化因子 1/√10 ≈ 0.3162。这个数看着不起眼但它是整个仿真里最容易翻车的地方。一旦符号功率不是 1你手算的噪声标准差就对不上画出来的 BER 曲线会整体平移甚至形状都变——低信噪比区看着还行高信噪比区直接“触底”不降。我一般会在映射函数里就做完归一化并在脚本开头加一行自检sym my_qammod16(bitsMat); assert(abs(mean(abs(sym).^2) - 1) 1e-3, 符号功率未归一化检查 1/sqrt(10));这行断言比任何注释都管用。符号数少的时候均值会有抖动所以阈值放到 1e-3样本量取到 1e5 以上就稳了。2.3 格雷映射让一个符号错误只赔一个比特格雷映射的定义是任意两个相邻星座点对应的比特组只差一位。在矩形 16QAM 里可以分解实现——I 路的 2 个比特按格雷顺序映射到 4 个电平Q 路同理二维最近邻就自然只差 1 bit。一维 4-PAM 的格雷表只有四行记住它比记公式有用b1 b2电平0 0-30 1-11 131 01注意表里的顺序不是二进制自然序11 排在 10 前面这样才能保证相邻电平只差一位。二维展开后16 个星座点的比特分配如下行是 b1b2列是 b3b4b1b2 \ b3b40001111000-3-3j-3-1j-33j-31j01-1-3j-1-1j-13j-11j113-3j3-1j33j31j101-3j1-1j13j11j格雷映射的收益是实打实的符号错误率 SER 在格雷映射下近似等于 4 倍的比特错误率 BER也就是 BER ≈ SER/4非格雷映射时一个最近邻错误可能带 2 个甚至 3 个比特错误BER 曲线整体右移看起来像“仿真做错了”其实只是映射表排错。判断方法很简单把映射表画成网格逐点检查上下左右的邻居汉明距离必须都是 1。2.4 一个不依赖工具箱的映射与解映射实现MATLAB 的 Communications Toolbox 里有qammod/qamdemod用起来省事但不少教学环境和轻量部署没有这个工具箱而且qammod的比特输入顺序、归一化开关在不同版本里行为有差异。我一般自己写一对函数逻辑透明出问题能逐行查。function sym my_qammod16(bitsMat) % 输入 bitsMat: N x 4每行一个符号第 1 列为 MSB % 输出 sym: N x 1 复符号平均功率归一化为 1 pam [-3, -1, 3, 1]; % 格雷电平表索引 2*b1b21 idxI 2*bitsMat(:,1) bitsMat(:,2) 1; % I 路 2 bit - 电平索引 idxQ 2*bitsMat(:,3) bitsMat(:,4) 1; % Q 路 2 bit - 电平索引 sym (pam(idxI). 1j*pam(idxQ).) / sqrt(10); end function bitsMat my_qamdemod16(rx) % 输入 rx: N x 1 复符号含噪 % 输出 bitsMat: N x 4与调制端比特顺序一致 lv [-3, -1, 3, 1]; % 判决电平顺序必须与 pam 一致 lut [0 0; 0 1; 1 1; 1 0]; % 电平索引 - 比特对 rI real(rx) * sqrt(10); % 先还原到 ±1/±3 域再判决 rQ imag(rx) * sqrt(10); [~, iI] min(abs(rI - lv), [], 2); % 最小欧氏距离判决 [~, iQ] min(abs(rQ - lv), [], 2); bitsMat [lut(iI,:), lut(iQ,:)]; endmy_qammod16的关键在pam数组的顺序它是格雷序不是升序写成升序会让 11 和 10 的电平互换BER 曲线立刻恶化。idxI的算法2*b1b21把两位二进制映射到 1~4正好对上数组下标。my_qamdemod16里先乘回 √10 再判决这一步很多人省掉直接对归一化后的 ±0.316 和 ±0.949 判决结果因为门限没跟着缩放而大面积误判。提示如果环境里有 Communications Toolbox可以用qammod(bits, 16, InputType, bit, UnitAveragePower, true)做交叉验证两条路径出来的 BER 曲线应该逐点重合误差在统计涨落范围内。对不上就说明自己的映射表有问题。3. 把星座图从理想画到带损伤发射链路与 EVM 诊断星座图的价值不在“好看”在于它能在 BER 曲线出问题之前先给出线索点云是圆的还是椭圆的、有没有旋转、四个角的点是不是比中心的点散得更开。这一章把发射端、信道、接收端串起来每一步都能单独验证。3.1 比特源与符号功率自检第一步只做映射不加噪确认星座点位置和你脑子里的图一致。比特生成用randi注意矩阵形状——reshape是列优先的直接reshape(bits, [], 4)会把比特按错误的顺序分组这是一个高频翻车点。rng(20240501); % 固定种子保证每次跑出来的图一样 N 20000; % 符号数 k 4; % 16QAM 每符号 4 bit bits randi([0 1], N*k, 1); bitsMat reshape(bits, k, []).; % 关键先 reshape 成 4 x N 再转置 sym my_qammod16(bitsMat); fprintf(平均符号功率 %.4f\n, mean(abs(sym).^2)); scatter(real(sym), imag(sym), 40, filled); grid on; axis equal; xlabel(I); ylabel(Q); title(理想 16QAM 星座图);reshape(bits, k, [])得到 4×N 的矩阵每一列是一个符号的 4 个比特转置后每一行才是符号。如果写成reshape(bits, [], k)等于把前 4 个比特拆成 4 个符号的第一个比特映射出来的星座点看着还是 16 个但比特和符号的对应关系全乱BER 会稳定在 0.5 附近——这时候别怀疑噪声先怀疑 reshape。跑完看输出的平均符号功率应该是 1.0000 附近。如果输出 10说明归一化因子漏了如果输出 0.1说明除了两次 √10。3.2 噪声方差自己算别把信噪比定义交给函数awgn函数有三个默认行为容易踩默认按measured估计信号功率、默认输入是实数、默认 SNR 指的是符号功率与噪声功率之比而不是 Eb/N0。一旦你的符号没归一化或者混用了 Es/N0 和 Eb/N0曲线就会平移。我习惯自己算标准差代码里写清楚每个量的定义。复基带下噪声总功率 N0 Es / (Es/N0)实部和虚部各分一半所以每一维的标准差是 √(N0/2)。符号已归一化时 Es 1EbN0_dB 10; % 每比特信噪比 EsN0_dB EbN0_dB 10*log10(k); % 16QAM: Es/N0 Eb/N0 6.02 dB EsN0 10^(EsN0_dB/10); sigma sqrt(1/(2*EsN0)); % 每一维噪声标准差 noise sigma * (randn(N,1) 1j*randn(N,1)); rx sym noise;参数说明10*log10(k)这一项是 16QAM 的核心换算k4 时等于 6.02 dB漏掉它曲线整体左移 6 dB看上去“仿真比理论还好”实际是偷了功率。sigma前面的 1/2 来自复噪声两维平分写成sqrt(1/EsN0)会让噪声功率翻倍曲线右移 3 dB。这两个常数我见过太多人写反建议在代码里加注释固定下来。3.3 EVM 与散点图先看形态再看数值加噪之后画散点同时算 EVM。EVM 的 RMS 定义是误差矢量功率与参考符号功率之比的平方根在只有 AWGN 的情况下它和 Es/N0 是一一对应的EVM √(1/(Es/N0))。evm_rms sqrt(mean(abs(rx - sym).^2) / mean(abs(sym).^2)); evm_theory sqrt(1/EsN0); fprintf(EVM 仿真 %.2f%%理论 %.2f%%\n, evm_rms*100, evm_theory*100); idx randperm(N, 3000); % 抽样画图避免百万点糊成一团 scatter(real(rx(idx)), imag(rx(idx)), 8, .); hold on; plot(real(sym), imag(sym), r, MarkerSize, 10, LineWidth, 1.5); grid on; axis equal; xlabel(I); ylabel(Q);EVM 仿真值和理论值应该吻合到小数点后一位样本量两万时误差在 0.3% 以内。如果 EVM 对得上但 BER 对不上问题在判决或比特映射如果 EVM 就对不上问题在噪声方差或归一化。这个二分法能省掉大量排查时间。散点图的形态信息量很大点云呈标准圆形扩散说明只有 AWGN沿某个方向拉长说明有 IQ 不平衡整体旋转一个角度说明有载波相位偏差四个角上的点比中心散得更开说明有非线性失真比如功放压缩。这几种形态在只看 BER 曲线时是完全看不出来的。3.4 接收端判决与误比特统计判决用最小欧氏距离和 2.4 节的解映射函数配合。统计误比特时要注意比较的两个数组展开顺序必须一致否则sum(a ~ b)算出来的是无意义的数。bitsHat my_qamdemod16(rx); errBits sum(bitsMat(:) ~ bitsHat(:)); ber errBits / numel(bitsMat); fprintf(BER %.3e%d / %d bit\n, ber, errBits, numel(bitsMat));bitsMat(:)和bitsHat(:)都是按列优先展开只要两者形状相同都是 N×4逐元素比较就是对的。如果中间某一步把bitsHat转置过形状变成 4×NMATLAB 不会报错会静默地做广播比较得到一个完全错误的 BER——这种错误最难查因为代码不崩。单点验证做完就可以进循环扫信噪比了。这一步的意义是在投入长时间仿真之前先确认一个点的 BER 落在理论值附近避免跑了两小时才发现映射表排错。4. 误码率仿真蒙特卡洛循环、理论曲线与置信区间BER 仿真的本质是蒙特卡洛计数发 N 个比特数错了多少个比值就是估计值。难点不在这个比值在于“N 取多大才可信”和“曲线末端为什么总是断掉”。4.1 理论曲线先有一个能对照的靶子矩形 16QAM 在 AWGN、格雷映射、理想同步条件下的 BER 有一个好用的近似闭式解BER ≈ (3/8) · erfc( √( (2/5) · Eb/N0 ) )这个式子的来源是一维 4-PAM 的符号错误概率近似为 1.5·Q(√(0.4·Es_pam/N0))两维合计得到符号错误率再除以每符号 4 比特并利用最近邻只差 1 bit 的性质。它在 BER 1e-1 的区间内与精确值吻合得很好做工程对照完全够用。function ber ber16qam_theory(EbN0_dB) % 矩形 16QAM 格雷映射在 AWGN 下的 BER 近似 EbN0 10.^(EbN0_dB/10); ber (3/8) * erfc(sqrt((2/5) * EbN0)); end参数说明2/5这个系数等价于3/(2(M-1))·log2(M)/...化简后的结果M16 时就是 0.4别自己改成3/(2*(M-1))之外的形式。erfc是 MATLAB 基础函数不需要任何工具箱。这条曲线在 BER1e-3 处对应 Eb/N0 ≈ 10.5 dB可以作为你判断仿真对错的锚点——如果你的仿真曲线在 1e-3 处落在 7 dB基本可以断定 Es/N0 和 Eb/N0 混用了。4.2 循环骨架按错误数停不按符号数停固定符号数的循环在低信噪比区浪费算力在高信噪比区又攒不够错误。正确的做法是双条件要么攒够目标错误数要么达到符号数上限。EbN0_dB 0:2:14; targetErr 200; % 每个点至少 200 个错误比特 maxSym 2e6; % 单点符号上限防止高信噪比死循环 berSim nan(size(EbN0_dB)); for idx 1:numel(EbN0_dB) EsN0 10^((EbN0_dB(idx) 10*log10(k))/10); sigma sqrt(1/(2*EsN0)); errBits 0; nSym 0; while errBits targetErr nSym maxSym Nb 2e5; % 每批符号数 b randi([0 1], Nb*k, 1); bm reshape(b, k, []).; s my_qammod16(bm); r s sigma*(randn(Nb,1) 1j*randn(Nb,1)); bh my_qamdemod16(r); errBits errBits sum(bm(:) ~ bh(:)); nSym nSym Nb; end berSim(idx) errBits / (nSym*k); fprintf(Eb/N0%2d dB, 符号%d, 错误%d, BER%.3e\n, ... EbN0_dB(idx), nSym, errBits, berSim(idx)); end semilogy(EbN0_dB, berSim, o-, LineWidth, 1.2); hold on; semilogy(EbN0_dB, ber16qam_theory(EbN0_dB), k--, LineWidth, 1.2); grid on; xlabel(E_b/N_0 (dB)); ylabel(BER); legend(仿真, 理论近似, Location, southwest);targetErr 200对应 BER 估计的相对标准差约 1/√200 ≈ 7%画在对数坐标上几乎看不出抖动。maxSym 2e6是为了在 BER 低于 1e-5 时主动收手——那个区间要攒够 200 个错误得跑几十亿比特单机不划算用理论曲线外推更务实。每批 2e5 个符号是内存和循环开销的折中太小则循环次数多太大则单次分配慢。批次内的四个步骤顺序不能变生成比特、映射、加噪、判决。中间任何一步用了上一批的残留变量比如 sigma 没更新曲线就会出现台阶状异常。4.3 置信区间让曲线末端有话说曲线在高信噪比处断掉不是仿真错了是错误数不够。错误数 Nerr 对应的 BER 估计其 95% 置信区间半宽约为 1.96/√Nerr。这个数字决定了你该不该相信曲线末端那个点错误比特数相对标准差95% 置信区间半宽1032%±62%10010%±20%4005%±10%10003.2%±6.2%工程上我一般要求每点至少 100 个错误重要节点比如用来定指标的 1e-4攒到 400 个以上。如果某个点达到符号上限时错误数还是 0正确做法是把这个点标成“无错误”而不是画成 0——对数坐标下的 0 会被 MATLAB 丢掉看起来像曲线提前结束实际上是你没有数据。4.4 加速把画图移出循环仿真慢通常不是慢在数学上是慢在每轮都调scatter、semilogy或者动态扩容数组。几条经验绘图命令全部移到循环外用nan预分配结果数组批量生成比特而不是逐比特循环有条件时开parfor按信噪比点并行每个点用独立随机流注意rng的用法别让各 worker 拿到同一个种子。还有一个容易忽略的开销randn的调用次数。把randn(Nb,1) 1j*randn(Nb,1)写成一次调用生成 2×Nb 再拆分比两次调用略快但代码可读性下降除非真的卡在性能上否则不建议为了这点收益牺牲清晰度。5. 避坑与排查五条来自实际仿真的记录这一章记录的是我调到凌晨还在盯屏幕的那几次每条都按现象、原因、解决写能对号入座就省几小时。5.1 归一化与信噪比定义上的三个坑第一条BER 曲线整体左移 6 dB看着比理论还好。现象仿真曲线在 1e-3 处对应的 Eb/N0 比理论低 6 dB 左右形状倒是对的。 原因把 Es/N0 直接当成 Eb/N0 用了漏掉10*log10(4) 6.02 dB的换算。噪声按 Eb/N0 算但符号功率按 Es 走等于凭空多给了 6 dB 能量。 解决在代码里固定一个换算函数EsN0_dB EbN0_dB 10*log10(k)所有噪声计算只接受 Es/N0不允许两套变量混用。第二条低信噪比吻合高信噪比曲线“触底”不再下降。现象BER 到 1e-4 附近就压不下去加大仿真量也没用散点图上看点云也不算太散。 原因符号没有归一化平均功率是 10 而不是 1但噪声方差是按 Es1 算的等效信噪比被抬高了一个固定量判决门限又按归一化后的电平设的两头不匹配。 解决在映射函数里做归一化脚本开头加功率断言判决函数里先乘回 √10 再比较。第三条EVM 对得上BER 差一个数量级。现象EVM 仿真值 12.5%理论值 12.4%看着没问题但 BER 比理论高十倍。 原因判决电平数组的顺序与调制电平表的顺序不一致。调制用格雷序[-3,-1,3,1]判决用了升序[-3,-1,1,3]导致 1 和 3 两个电平被互换。 解决把电平表和比特映射表定义成同一个函数里的常量或者用调制端同一份数组做判决禁止两处各写一遍。5.2 比特映射与代码实现上的两个坑第四条BER 稳定在 0.5 附近几乎没有正确比特。现象不管信噪比多高BER 都在 0.4 到 0.5 之间星座图看着还是 16 个点。 原因reshape的维度顺序错了写成reshape(bits, [], 4)比特和符号的对应关系被彻底打乱。映射表本身没错所以星座点位置正常但接收端恢复出来的比特顺序和发送端对不上。 解决统一用reshape(bits, k, []).并在第一次运行时用isequal验证“映射后再解映射、无噪声”能 100% 还原。第五条脚本跑了一整夜没结束。现象任务管理器里 MATLAB 占满一个核进度条不动低信噪比点早就跑完了。 原因while 循环只设了错误数条件高信噪比区永远攒不够目标错误又没设符号上限同时循环内每批都调用scatter画图句柄越积越多。 解决加maxSym上限绘图移出循环把每批符号数从 2e4 提到 2e5 减少循环次数必要时按信噪比点parfor并行。注意排查顺序建议固定为“单点无噪验证 → 单点加噪验证 → 扫曲线”。无噪时要求 BER 严格为 0加噪单点要求与理论差在统计涨落内这两步过了再扫曲线否则就是在错误的基线上做长时间计算。6. 把 BER 反推成 EVM 预算仿真结果怎么用到指标上6.1 EVM 与 SNR 的换算关系仿真跑通之后最有用的一步是把 BER 目标反推成 EVM 指标——因为测试仪器上能直接读的是 EVM不是 Eb/N0。AWGN 占主导时两者的关系很干净EVM_rms √(1/(Es/N0))换成 dB 就是 SNR_dB −20·log10(EVM_linear)。EVM [0.20 0.15 0.125 0.10]; SNR -20*log10(EVM); % 对应 Es/N0单位 dB EbN0 SNR - 10*log10(4); % 换算到 Eb/N0 ber ber16qam_theory(EbN0); disp(table(EVM., SNR., EbN0., ber., ... VariableNames, {EVM,EsN0_dB,EbN0_dB,BER}));6.2 一张反查表目标 BER 对应多少 EVM按 4.1 节的近似式和上面的换算矩形 16QAM 在 AWGN、理想同步、格雷映射下的对应关系如下。这些数字可以当作链路预算的起点实测时再叠加实现余量。目标 BEREb/N0 (dB)Es/N0 (dB)EVM_rms1e-2约 7.9约 13.9约 20%1e-3约 10.5约 16.6约 15%1e-4约 12.2约 18.2约 12%1e-5约 13.4约 19.4约 11%用法是反的先定系统目标 BER比如 1e-4查表得到需要的 Es/N0 约 18.2 dB再折算成 EVM 约 12%。实测 EVM 如果只有 15%说明还差 2 dB 左右的余量要么改调制阶数要么在前端找损耗。注意这张表只算 AWGN真实系统里的相位噪声、IQ 不平衡、功放非线性都会让同样 EVM 下的 BER 更差所以工程上一般再留 1.5 到 3 dB 的实现余量。6.3 验证套路三条曲线必须同时成立我现在验收一个 16QAM 仿真工程会同时看三样东西理想星座图16 个点位置精确、功率为 1、某一信噪比下的散点图与 EVM与理论值差在 0.5% 以内、BER 曲线与理论近似在统计涨落内重合。三样里任何一样对不上就回到对应的单点验证去查而不是继续扫曲线。还有一个小习惯每次改完参数先用rng固定种子跑一遍短样本把关键中间量平均功率、sigma、EVM、错误数打印出来存进日志再放开跑长样本。早期我嫌这一步麻烦直接上 2e6 符号结果一个reshape的顺序错误让我对着一条“怎么调都和别人不一样”的曲线折腾了整晚。从那以后我每次动映射表或者噪声方差都强制走一遍“无噪还原 → 单点 EVM → 扫曲线”这三步再急也不跳过。希望帮到你。本文还有配套的精品资源点击获取