★ FEATURED ARTICLE 简介16QAM星座图与误码率仿真资料面向通信原理、数字通信课程学习者及MATLAB仿真入门者用于理解星座映射、调制解调流程与BER评估。压缩包共8个文件均为.m脚本体积约3KB内容覆盖主流程调用、星座点绘制、四进制与二进制互转、QAM调制解调、随机二进制序列生成及比特成形等环节结构紧凑便于逐函数阅读和二次修改。已有1927人学习下载在相关实验与课程设计中具有一定参考价值。读者可据此搭建端到端仿真链路观察不同信噪比下AWGN信道中的星座分布计算并绘制误码率曲线掌握16QAM系统性能分析的基本方法也可将其作为课程实验或项目原型的代码骨架。1. 从一张糊成一团的星座图说起16QAM 星座图仿真到底在算什么调完qammod顺手scatterplot屏幕上 16 个点糊成一片四个角的点还往外飘换到误码率那一步曲线整体右移三五个 dB怎么都对不上理论线。“16QAM 星座图仿真和误码率仿真”这件事难的不是写几行 matlab而是让每一步都能自圆其说出了问题能定位到具体环节。这套工程把链路拆成比特生成、格雷映射、星座归一化、AWGN 加噪、硬判决解调、BER 统计、与berawgn理论曲线对照出图每一段都能单独拎出来验证。它适合正在做通信原理课设、要交一份能复现的仿真报告的人也适合已经上手工作、想快速搭一条基线链路去核对链路预算的开发者。2. 链路骨架一份 16QAM 仿真工程从入口到出图的调用关系仿真脚本最容易踩的问题不是算法错而是参数散落。M 定义在文件头归一化开关写在循环里理论曲线又用了另一套 Eb/N0 单位改一个忘一个最后曲线不对还得回头逐行找。我一般按“入口只负责参数和出图算法函数只负责单点计算绘图函数只负责画”这三层来切好处是换调制阶数不用重写循环扫新的信噪比范围也不用碰算法。2.1 目录结构与每个文件的输入输出按这个思路工程目录大概是这样qam16_sim/ ├── main_qam16.m % 入口参数定义、调用扫描、统一出图 ├── ber_qam16_awgn.m % 扫描循环逐点调用单点函数汇总 BER ├── ber_qam16_point.m % 单点映射 → 加噪 → 解调 → 计错 ├── plot_constellation.m % 星座图与判决域可视化 └── readme.txt % 参数说明与复现步骤main_qam16.m是唯一持有全局参数的地方它把M、EbNo_dB、Nbits打包传给ber_qam16_awgn.m后者只做循环与错误数累计真正的比特到符号、符号到比特的往返交给ber_qam16_point.m。这种切法的直接收益是可测试性怀疑映射有问题单独跑一段qammod/qamdemod回环不传噪声误码率必须是 0怀疑加噪有问题单独把awgn换成手工噪声看曲线是否跟着变。另一层考虑是复现。仿真类代码最怕“昨天跑出来是对的今天跑出来不一样”所以入口脚本第一行我会固定rng(2024)之类的种子。默认的随机数发生器在不同 MATLAB 版本之间有过变化跨版本对比曲线时把种子和发生器类型一起钉死能省掉不少扯皮。2.2 主入口脚本从参数到两张图入口脚本承担三件事定参数、调扫描、出图。写法上我尽量让它短短到一眼能看出整条链路的顺序。%% main_qam16.m —— 16QAM 星座图与误码率仿真主入口 clc; clear; close all; rng(2024); % 固定随机种子保证曲线可复现 M 16; % 调制阶数 k log2(M); % 每符号比特数16QAM 为 4 Nbits 2e5; % 单次调用的比特数 EbNo_dB 0:2:14; % 每比特信噪比扫描范围 minErr 100; % 每个信噪比点至少攒够的错误比特数 maxBit 4e6; % 单点比特上限避免高信噪比跑死 sps 1; % 符号级仿真不做过采样 % 1) 归一化参考星座先确认 16 个点摆对了 const qammod(0:M-1, M, UnitAveragePower, true); fprintf(平均符号功率 %.4f\n, mean(abs(const).^2)); % 应接近 1 % 2) 误码率扫描 [berSim, berTh] ber_qam16_awgn(EbNo_dB, M, Nbits, minErr, maxBit); % 3) 星座图与理论/仿真对照图 plot_constellation(const, M); figure; semilogy(EbNo_dB, berTh, k--, EbNo_dB, berSim, ro-, LineWidth, 1.2); grid on; xlabel(E_b/N_0 (dB)); ylabel(BER); legend(理论值, 仿真值, Location, southwest); ylim([1e-6 1]);这段里几个参数值得单独说。Nbits控制单次调用的粒度太小会让循环次数暴涨太大会让高信噪比点一次跑很久才凑够错误数2e5是个折中。minErr决定统计置信度100 个错误对应的相对标准差大约 10%曲线在瀑布区不会抖得太难看。maxBit是保险丝14 dB 以上 16QAM 的 BER 已经在 1e-5 量级不放上限的话一个点能跑十几分钟。sps这里固定为 1因为只做符号级误码统计不做波形和滤波过采样除了拖慢速度没有任何收益如果后面要加根升余弦成形这个参数才需要放开。fprintf那行是刻意留的体检项。16QAM 矩形星座的 I、Q 各取 ±1、±3平均符号能量是 10归一化后必须接近 1。这一行打印出 10 或者 0.1后面的曲线就不用看了直接回头查归一化。2.3 参数面板改哪个数会带来什么后果把可变参数集中之后最好再配一张表标明每个数动过之后会影响什么。这张表我一般直接写进readme.txt省得下一次打开工程还得重新推。参数含义常用取值改动后的直接影响M调制阶数4 / 16 / 64改变每符号比特数与星座规模理论曲线需同步换k每符号比特数log2(M)影响 Es/N0 与 Eb/N0 的换算写死会平移曲线Nbits单次比特数1e5 ~ 1e6影响循环次数与单点耗时EbNo_dB信噪比范围0:2:14范围越宽高信噪比端越难收敛minErr最小错误数100 ~ 200越大曲线越平滑耗时线性增长maxBit单点比特上限1e6 ~ 1e7防止高信噪比死循环归一化开关UnitAveragePower全程true调制解调不一致会直接判错真正容易出事的是最后一行。qammod和qamdemod的归一化开关必须成对出现只在一头打开接收端拿到的判决门限就整体偏了表现是低信噪比误码率还行、高信噪比怎么也降不下去很像“地板效应”但根本不是信道问题。这类错误在星座图上不明显因为整体缩放不影响点的相对位置只有把判决边界画出来才看得出来。3. 星座图模块16 个点怎么摆、怎么归一化、怎么画才不糊星座图是这条链路里唯一“看得见”的东西也最容易让人误判。点画得好看不代表映射对点画得糊也不一定是映射错很多时候只是scatterplot的参数没调。这一章把坐标生成、格雷映射、判决域三件事拆开讲每一件都给一段能单独跑的代码。3.1 矩形 16QAM 的坐标与平均功率归一化矩形 16QAM 的构造很朴素I 路和 Q 路各自从 {-3, -1, 1, 3} 里取一个电平两两组合成 16 个点。因为每个电平出现的概率相同单路的平均功率是 (1199)/4 5两路加起来符号平均能量 Es 10。归一化就是整体除以 √10让 Es 变成 1这样后面谈信噪比时不会混进一个固定的功率因子。M 16; levels [-3 -1 1 3]; % 单路四个电平 [I, Q] meshgrid(levels, levels); % 两两组合 constRaw I(:) 1j*Q(:); % 16 个复数星座点未归一化 Es_raw mean(abs(constRaw).^2); % 平均符号能量应为 10 constNorm constRaw / sqrt(Es_raw); % 归一化到 Es 1 % 与工具箱的输出做对比两者必须一致顺序可能不同 constRef qammod(0:M-1, M, UnitAveragePower, true); fprintf(自建星座平均功率 %.4f\n, mean(abs(constNorm).^2)); fprintf(工具箱星座平均功率 %.4f\n, mean(abs(constRef).^2));meshgrid出来的顺序是 I 变化最快I(:)和Q(:)展开后得到 16 个组合前四个是 I 从 -3 到 3、Q 固定 -3正好是格雷映射里一条边的顺序。Es_raw一定是 10如果算出来不是多半是levels写成了 {-1, -1/3, 1/3, 1} 这种已经归一化过的值又除了一次 √10曲线会整体平移 10 dB这种错误比想象中常见。和工具箱结果对比时不要用isequal因为点的排列顺序不一定相同比较平均功率和最小间距更稳。最小间距归一化后是 2/√10 ≈ 0.632这个数在画判决域时要用到。3.2 格雷映射相邻星座点只差一个比特格雷映射的意义在于噪声让接收点跑到相邻星座点时只有一个比特判错从而把符号错误率近似换算成比特错误率。qammod在InputType设为bit时默认就走格雷映射但“默认”这两个字不能当作免检标签最好自己验一遍。M 16; const qammod(0:M-1, M, UnitAveragePower, true); % 每个星座点对应的 4 比特标签MSB 在前 idx qamdemod(const, M); bits dec2bin(idx, log2(M)) - 0; % 16×4 的 0/1 矩阵 % 计算每个点到最近邻的距离并统计标签的汉明距离 D abs(const - const.); % 16×16 距离矩阵 D(1:M1:end) inf; % 对角线置无穷排除自身 [~, nn] min(D, [], 2); % 每个点的最近邻索引 hamDist sum(bits ~ bits(nn, :), 2); % 与最近邻的汉明距离 fprintf(最近邻汉明距离最大值 %d\n, max(hamDist)); % 格雷映射下应为 1这段代码的逻辑是先用qamdemod把星座点还原成 0~15 的整数标签dec2bin转成 4 位二进制再对每个点找几何上最近的邻居比较两者的比特串。如果映射不是格雷码最大值会到 2 甚至 3。注意bits(nn,:)用的是行索引nn是 16×1 的列向量直接索引会得到一个 16×4 的矩阵sum(...,2)按行求和就是每个点的汉明距离。有一个细节容易被忽略qammod的比特分组是按输入向量的顺序每 k 个比特组成一个符号第一个比特是最高位。如果上游用reshape把比特流排成矩阵再转置分组顺序就变了解调回来的比特串和原始比特串对不上误码率会凭空抬高一大截而且高信噪比端不会收敛到 0。这种错误在星座图上看不出任何异常只能靠无噪回环来抓。3.3 判决域可视化把“最近邻”画出来星座图只画点看不出接收端的判决规则。加一张判决域的热力图能直观看到每个点管辖的区域也能验证归一化之后的坐标轴范围。M 16; const qammod(0:M-1, M, UnitAveragePower, true); % 在复平面上铺网格逐点做最近邻判决 lim 1.6; [Re, Im] meshgrid(-lim:0.02:lim, -lim:0.02:lim); gridPts Re(:) 1j*Im(:); D abs(gridPts - const.); % N×16 距离矩阵 [~, label] min(D, [], 2); % 最近邻的星座点索引 region reshape(label, size(Re)); figure; imagesc(-lim:0.02:lim, -lim:0.02:lim, region); axis xy; axis equal; hold on; plot(real(const), imag(const), k, MarkerSize, 10, LineWidth, 1.5); xlabel(I); ylabel(Q); title(16QAM 判决域最近邻);abs(gridPts - const.)利用隐式扩展得到一个 N×16 的矩阵第 i 行是第 i 个网格点到 16 个星座点的距离min(...,[],2)取每行最小值就是最近邻。reshape回网格尺寸后用imagesc上色axis xy保证纵轴向上不然图会上下翻转看起来像映射错了。lim取 1.6 是因为归一化后最外圈点在 ±3/√10 ≈ ±0.949留一点余量能看到边界。如果星座图上的点超出这个范围说明归一化没生效平均功率还是 10这时候判决域图的边界会把外侧点切掉一眼就能发现。4. 误码率仿真引擎蒙特卡洛循环与 Eb/N0 换算误码率部分的核心是一句话在给定的 Eb/N0 下发足够多的比特数错的比例。听起来简单但“足够多”和“给定”这两个词各自埋着坑。足够的样本量决定了曲线尾部是否可信Eb/N0 与 Es/N0 的换算决定了曲线位置是否正确。这一章把单点函数、扫描循环、理论对照分开写。4.1 单点误码统计函数单点函数只做一件事给定一个 Eb/N0返回这个点的误码率和错误比特数。错误比特数要一并返回因为上层的循环需要用它来判断是否继续加样本。function [ber, nErr, nBit] ber_qam16_point(EbNo_dB, M, Nbits) % 输入EbNo_dB 每比特信噪比(dB)M 调制阶数Nbits 本次比特数 % 输出ber 误码率nErr 错误比特数nBit 实际统计比特数 k log2(M); EsNo_dB EbNo_dB 10*log10(k); % 关键换算Es/N0 Eb/N0 × k txBits randi([0 1], Nbits, 1); % 独立等概比特流 txSym qammod(txBits, M, InputType, bit, ... UnitAveragePower, true); rxSym awgn(txSym, EsNo_dB, measured); % 复基带 AWGN rxBits qamdemod(rxSym, M, OutputType, bit, ... UnitAveragePower, true); nErr sum(rxBits ~ txBits); ber nErr / Nbits; nBit Nbits; endEsNo_dB EbNo_dB 10*log10(k)这一行是整段代码里最不能写错的地方。16QAM 的 k 410*log10(4) ≈ 6.02也就是说同一个 Es/N0 对应的 Eb/N0 要低 6 dB。如果这里忘了加曲线会整体左移 6 dB看上去“仿真比理论还好”很容易被误判成代码写得不错。awgn的第三个参数用measured表示按输入信号的实际功率测量后再加噪。因为前面已经归一化到 Es 1理论上可以直接省略这个参数但实际工程里如果中途插入了成形滤波或者增益调整功率就不再是 1measured更稳。代价是每次调用都要重新估计功率样本很少时估计偏差会让曲线轻微抖动所以Nbits不宜取得太小。qamdemod的UnitAveragePower必须和qammod保持一致这一条在前面提过但在函数里特别容易漏因为两个函数的参数列表长得不一样复制粘贴时很容易只改了一处。4.2 扫描循环与错误数下限外层循环负责遍历 Eb/N0并在每个点上累积到足够的错误数。固定比特数的写法在高信噪比端会浪费大量时间在低信噪比端又可能样本不足所以用“错误数下限 比特数上限”的双重条件更合理。function [berSim, berTh] ber_qam16_awgn(EbNo_dB, M, Nbits, minErr, maxBit) berSim zeros(size(EbNo_dB)); for ii 1:numel(EbNo_dB) nErr 0; nBit 0; while nErr minErr nBit maxBit [~, e, n] ber_qam16_point(EbNo_dB(ii), M, Nbits); nErr nErr e; nBit nBit n; end if nErr 0 berSim(ii) NaN; % 一个错都没有这个点不可信 else berSim(ii) nErr / nBit; end fprintf(Eb/N0 %4.1f dB | 错误 %6d / 比特 %9d | BER %.3e\n, ... EbNo_dB(ii), nErr, nBit, berSim(ii)); end berTh berawgn(EbNo_dB, qam, M); % 第一参数单位为 dB endwhile的条件顺序有讲究先判断错误数再判断比特上限。如果反过来低信噪比点可能因为先撞上上限而提前退出虽然影响不大但语义上不清晰。berSim(ii) NaN这个处理是为了在图上明确标出“没测到”而不是画一个 0 出来semilogy遇到 0 会直接断线看起来像曲线掉了。berawgn的第一个参数是 dB 值不是线性值。这一点和很多教科书公式不一样公式里用的是线性 Eb/N0写代码时如果传了10^(EbNo_dB/10)理论曲线会跑到天上。我见过不止一次这种翻车表现是理论线在 0 dB 处就是 1e-8仿真线还在 1e-1两条线隔着半个屏幕。打印语句不是装饰。每次运行把每个点的错误数和总比特数打出来跑完扫一眼就知道哪个点样本不够错误数刚好卡在 100 的点曲线会明显偏离相邻点这种点在报告里要么补样本要么直接标注置信区间。4.3 理论曲线对齐与出图理论曲线除了调berawgn最好再手算一遍公式两条理论线重合才能确认工具箱的调用方式没理解错。16QAM 在格雷映射下的比特误码率近似为M 16; k log2(M); EbNo_dB 0:2:14; EbNo 10.^(EbNo_dB/10); % 转成线性值代入公式 Pb_approx (4/k) * (1 - 1/sqrt(M)) .* ... qfunc(sqrt(3*k/(M-1) * EbNo)); Pb_toolbox berawgn(EbNo_dB, qam, M); % 注意这里是 dB figure; semilogy(EbNo_dB, Pb_approx, b-, EbNo_dB, Pb_toolbox, k--); grid on; xlabel(E_b/N_0 (dB)); ylabel(BER); legend(解析近似, berawgn);qfunc是 MATLAB 自带的 Q 函数3*k/(M-1)对 16QAM 是 0.8代进去就是常见的 0.75·Q(√(0.8·Eb/N0))。两条线在中高信噪比几乎重合低信噪比处会有细微差别因为解析式是高信噪比近似berawgn用的是更精确的符号错误率转换。知道这个差别看到两条线在 0~4 dB 分开一点就不会慌。出图时ylim([1e-6 1])是必须的不然某个点没测到错误会画出 0坐标轴被拉到 1e-300整张图什么都看不见。图例位置放在southwest是因为曲线从左下往右下走右上角空着图例放那儿不挡线。5. 避坑与排查曲线对不上理论值的五个现场仿真跑完对不上理论线是这条链路最常见的状态。下面五条是我自己踩过、也帮别人排查过的现场按“现象 → 原因 → 解决”记下来下次遇到直接对号入座比从头读代码快得多。5.1 曲线整体平移先怀疑信噪比定义现象仿真曲线形状和理论线几乎平行只是整体右移了 6 dB 左右或者右移了 10 dB。原因右移 6 dB 基本可以锁定 Es/N0 与 Eb/N0 混用。16QAM 每符号 4 比特10*log10(4) 6.02忘了这一步换算就会差这么多。右移 10 dB 则是归一化没做矩形 16QAM 的原始平均功率是 10awgn按 Es 10 的信号加噪等效信噪比就低了 10 dB。解决在单点函数入口打印EsNo_dB - EbNo_dB16QAM 应该稳定输出 6.02在参考星座上打印mean(abs(const).^2)应该稳定输出 1。两个数任意一个不对先修这里别去动信道模型。5.2 曲线形状异常比特分组和判决口径现象低信噪比端误码率比理论高很多高信噪比端却怎么也不往下掉曲线出现明显的“地板”。原因多半是比特分组顺序错位。用reshape把比特流排成 k 行的矩阵再转置分组顺序就变了qammod按新的顺序打包qamdemod按原顺序拆包两边对不上。另一种可能是统计口径混用拿符号错误率当比特错误率或者反过来。解决先做无噪回环把awgn去掉误码率必须是 0不是 0 就说明映射与解映射不匹配。统计比特错误率时用sum(rxBits ~ txBits)不要用sum(rxSym ~ txSym)去数符号浮点噪声下两个复数几乎不可能精确相等那样数出来的符号错误率恒为 1。5.3 跑得慢、尾部抖样本量与上限现象14 dB 那个点跑了十几分钟还没结束或者曲线在最后几个点上下跳动相邻点之间不单调。原因高信噪比下错误事件稀疏固定比特数要么跑不完要么凑不够错误数。抖动则是样本量不足的直接表现100 个错误对应的相对标准差约 10%50 个错误就接近 14%在semilogy上会被放大成肉眼可见的折线。解决用minErr maxBit双条件minErr给 100 到 200maxBit给 1e6 到 1e7。跑不到错误数的点直接标NaN在报告里注明“该点未观测到错误”比硬画一个点更诚实。另外把Nbits调大能减少循环次数但单次调用的内存占用也会上升2e5到1e6之间比较舒服。5.4 星座图糊成一团绘图参数与坐标轴现象scatterplot出来的图 16 个点挤在中心四周全是散点看不出星座结构。原因直接把几千上万个含噪符号丢进去噪声点的密度远高于星座点本身视觉上就是一团。坐标轴自适应范围被离群点拉大星座点显得更小。解决只画前 1000 到 5000 个符号并用axis手动限定在 ±1.5 左右或者先画无噪参考星座再把含噪点叠上去用不同颜色区分。想看判决效果用第 3 章那段判决域代码比散点图直观得多。5.5 换调制阶数后直接报错现象把M从 16 改成 64代码在qammod处报维度不匹配。原因InputType为bit时输入比特长度必须是log2(M)的整数倍。16QAM 时Nbits是 4 的倍数没问题64QAM 需要 6 的倍数原来的2e5除不尽最后一个符号的比特不够。解决在入口处加一句Nbits floor(Nbits/log2(M))*log2(M)把比特数向下取整到 k 的倍数。同时别忘了minErr和理论曲线的调用也要跟着换berawgn的第二个参数从qam到M都要同步64QAM 的平均符号能量是 42归一化因子跟着变。6. 二次开发把固定 16QAM 扩成可切换调制与衰落信道这套代码写顺之后最自然的改造方向是把M参数化让它能在 QPSK、16QAM、64QAM 之间切换再挂一个平坦衰落信道看看相干检测下的曲线怎么走。改造成本不高因为前面已经把所有跟阶数相关的东西集中在入口和三处换算里。调制方式每符号比特 k归一化前 Es理论曲线调用备注QPSK22berawgn(EbNo,psk,4,nondiff)与 4QAM 等价可直接对比16QAM410berawgn(EbNo,qam,16)基线配置64QAM642berawgn(EbNo,qam,64)高信噪比端收敛更慢改造时唯一需要动脑的是衰落信道那段。加性白噪声的写法可以直接复用衰落要把信道系数乘在符号上并且注意噪声功率要按发射端的 Es/N0 来定不能乘完衰落再用measured测功率否则每个符号的等效信噪比都在变测出来的曲线会偏低。function ber ber_point_fading(EbNo_dB, M, Nbits) k log2(M); Nsym Nbits / k; EsNo 10^(EbNo_dB/10) * k; % 线性 Es/N0 txBits randi([0 1], Nbits, 1); txSym qammod(txBits, M, InputType,bit,UnitAveragePower,true); % 平坦瑞利衰落逐符号独立E[|h|^2] 1 h (randn(Nsym,1) 1j*randn(Nsym,1)) / sqrt(2); % 复高斯噪声实部虚部各占一半功率 noise sqrt(1/(2*EsNo)) * (randn(Nsym,1) 1j*randn(Nsym,1)); rxSym h .* txSym noise; % 先过信道再加噪 rxSym rxSym ./ h; % 理想信道估计下的相干检测 rxBits qamdemod(rxSym, M, OutputType,bit,UnitAveragePower,true); ber mean(rxBits ~ txBits); endnoise那一行的系数是sqrt(1/(2*EsNo))因为复噪声的总方差 N0 1/EsNo实部和虚部各分一半标准差就是sqrt(1/(2*EsNo))。如果直接写成sqrt(1/EsNo)噪声功率翻倍曲线整体右移 3 dB又是一个“看着像但就是差一点”的现场。rxSym ./ h用的是理想信道估计实际系统里这一步要换成导频估计误差会让曲线出现不可消除的地板仿真阶段先用理想值把上界跑出来再逐步加估计误差比一上来就上完整链路好定位。验证方法我习惯分三步走第一步h全置 1确认衰落函数退化成 AWGN 并且和原曲线完全重合第二步关掉噪声只留衰落确认误码率随 Es/N0 单调下降且不出现异常第三步才是打开全部和理论瑞利曲线对照。三步里任何一步对不上问题范围就缩到很小一块。从那以后我每次改完调制阶数或者信道模型都强制先跑一遍 0 dB 单点、把星座图和判决域打印出来看一眼确认 16 个点摆在正确位置、平均功率是 1再放开扫描。这个习惯救过我好几次因为绝大多数“玄学”曲线问题根源都在最开始那几行参数上。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
