ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

QPSK误码率仿真:信噪比对齐与Eb/N0换算的MATLAB实现

QPSK误码率仿真:信噪比对齐与Eb/N0换算的MATLAB实现 简介这是一份由达摩老生整理的Matlab程序源码专门用于计算QPSK调制下的误码率随信噪比变化趋势并绘制对应曲线。资源适合通信方向的新手以及有一定开发经验的工程师使用可作为课程设计、毕业设计或科研入门时的参考模板。包内共2个文件1个m脚本完整实现了从QPSK调制、信道加噪、相干解调到误码率统计的流程关键步骤均配有注释便于读懂和修改另1个docx文档则补充说明原理、参数设置与运行方式辅助快速上手。整个压缩包仅13KB结构精简、下载便捷不需要额外安装大型依赖。截至目前已有854人学习使用运行后可直接得到误码率曲线也能对照代码理解数字调制、信噪比与误码率之间的关系适合快速验证、对比实验以及二次开发。1. 从零写 QPSK 误码率仿真信噪比没对齐曲线永远对不上QPSK 误码率随信噪比变化的仿真十个新手跑第一版七八个会对着屏幕怀疑人生明明代码没报错BER 曲线和课本理论值就是差一截。这通常不是 MATLAB 用法问题而是大家嘴里都叫「信噪比」的那几个量实际上各说各话。以 matlab 计算 QPSK 误码率随信噪比变化为主题的这类程序源码核心其实只有三件事把随机比特映射成 QPSK 符号、在指定 Eb/N0 下叠上复数高斯噪声、判决后统计误码个数。适合通信原理课设、毕业设计里需要出 BER 曲线图的同学也适合做链路预算前先用脚本验证调制方案的人。下面按我平时最顺手的写法把这些拆开讲代码可以直接拷去跑顺带把最常见的几类翻车原因一次说清。2. 搭一条最小 QPSK 基带链路发射、加噪、判决的完整代码2.1 格雷映射把随机比特变成 QPSK 符号通信系统工具箱里现成的pskmod一行就能完成调制但自己手写映射有两个好处一是能清楚看到每个符号对应哪两个比特二是后面查误码统计时不会对着工具箱的黑匣子发懵。先固定一组基础参数% 基础参数 M 4; % QPSK 调制阶数 k log2(M); % 每符号携带 k2 比特 numBits 1e6; % 参与统计的比特总数 numSym numBits / k; % 对应 5e5 个 QPSK 符号 rng(0); % 固定随机流保证仿真可复现 data randi([0 1], numBits, 1); % 原始 0/1 比特流 dataSym reshape(data, k, []).; % 每行 2 比特对应一个符号 symIdx dataSym(:,1) * 2 dataSym(:,2); % 高位在前得到 0~3 的十进制索引 % 格雷映射星座相邻星座点只差 1 个比特 constellation (1/sqrt(2)) * [11i; -11i; -1-1i; 1-1i]; tx constellation(symIdx 1); % MATLAB 索引从 1 开始这里有几个关键点。symIdx的位序必须和映射表一致dataSym(:,1) * 2 dataSym(:,2)把「高位乘 2 加低位」变成了 0~3 的索引顺序是 00、01、11、10正好对应格雷码的相邻关系。星座点乘1/sqrt(2)是为了把符号能量归一化到 1这样每个星座点的模平方等于 1Es 的计算会非常干净。格雷映射不是玄学它保证相邻判决区域的两个星座点只差 1 比特。相邻符号判错是最常见的事件若用自然二进制映射一个符号判错可能带错 2 个比特BER 会比理论值高一截。这个细节在后面比对理论和仿真曲线时会明显体现出来。2.2 复数 AWGN 信道噪声功率怎么加才和理论一致基带等效仿真的 AWGN 信道加的是复数高斯白噪声不是实数噪声。很多对不上的曲线都是栽在这一步噪声实部、虚部的功率分配错了。对复数噪声n x jy要求总功率方差等于 N0那么实部和虚部分别占一半% 在某个 Eb/N0 下加噪声 EbN0 10; % 目标 Eb/N0单位 dB Es mean(abs(tx).^2); % 理想星座下就是 1 Eb Es / k; % 每比特能量QPSK 下是 0.5 N0 Eb / (10^(EbN0/10)); % 由 dB 值换算噪声功率谱密度 n sqrt(N0/2) * (randn(numSym,1) 1i*randn(numSym,1)); rx tx n;randn产生的是单位方差高斯噪声乘以sqrt(N0/2)后实部方差是 N0/2虚部方差也是 N0/2合成总功率才是 N0。如果图省事写成sqrt(N0) * randn实部虚部各多了一倍功率等效信噪比直接低了 3 dB这条曲线左移的 3 dB 会让你找半天原因。注意这里把 Es 换成mean(abs(tx).^2)而不是直接写 1是个好习惯。如果以后改成 16QAM、8PSK星座能量变了这段代码不必改。Eb 的换算也一并放在这里后面扫信噪比时直接复用。2.3 判决与误码统计一段能直接跑的完整脚本判决规则取决于星座怎么排。上面格雷映射里I 路实部正负号决定第一个比特Q 路虚部正负号决定第二个比特。real(rx) 0对应第一比特为 0imag(rx) 0对应第二比特为 0。把整个信噪比扫描、加噪、判决、统计串起来就是主循环EbNoVec (0:0.5:8); % 扫描 0~8 dB步长 0.5 berSim zeros(size(EbNoVec)); txFixed tx; % 发射序列固定只换噪声 for idx 1:length(EbNoVec) N0 Eb / (10^(EbNoVec(idx)/10)); n sqrt(N0/2) * (randn(numSym,1) 1i*randn(numSym,1)); rx txFixed n; b1 real(rx) 0; % I 路判决第一个比特 b2 imag(rx) 0; % Q 路判决第二个比特 est [b1 b2]; errTotal sum(sum(est ~ dataSym)); berSim(idx) errTotal / numBits; end发射序列在整个循环里保持不变每次只重新生成噪声这是蒙特卡洛仿真的常见做法能省掉重复生成随机比特的开销也不会引入发射端的额外随机性。sum(sum(est ~ dataSym))先按位比较得到逻辑矩阵两次求和后就是总误比特数。提示numBits一定要能被k整除。QPSK 的 k2所以 1e6 没问题如果你改成 1e51 这种奇数reshape会直接报错或悄悄丢掉最后一个比特。主循环跑完后可以先肉眼瞥一眼berSim0 dB 附近应该在 0.1 上下4 dB 附近应该到 1e-3 量级8 dB 附近应该在 1e-5 左右。如果数量级对不上先回去查噪声功率再查判决极性。3. 蒙特卡洛统计与 Eb/N0 换算仿真次数、步长和置信区间3.1 仿真比特数怎么定看误码个数而不是看比例误码率统计的可靠性不取决于「跑了多少比特」而取决于「统计到了多少个误码」。误码率 1e-4 时读到的比例看起来精确实际上可能只错了几十个比特噪声的随机性让这个点抖得厉害。经验法则是每个信噪比点至少积累 30~100 个误码相对误差才压得住。误码数 e 和相对标准差近似满足1/sqrt(e)e100 时相对偏差约 10%e30 时约 18%。目标误码率至少需要的比特数期望误码数1e-21e41001e-31e51001e-41e61001e-51e7100这张表的意思是想画到 8 dB 的 1e-5单点就要跑 1e7 比特。脚本里numBits 1e6跑到 1e-5 就明显吃力实际仿真到 1e-4 就该提示自己这个点以下只是趋势参考。低误码率点跑得慢是蒙特卡洛方法的物理限制不是 MATLAB 慢。与其硬跑不如在第 6 章讲的半分析法里找思路。另一点是 Eb/N0 的扫描步长。0.5 dB 步长既能看清曲线拐弯又不至于让高信噪比区间占用过多计算时间。0.1 dB 步长除了让曲线变密不会带来更多信息还让 1e-5 以下的点成倍增加耗时。3.2 Eb/N0、Es/N0 与 SNR三种叫法的换算关系这组换算关系是这类源码里最容易出错的点单独拎出来说。Eb/N0 是每比特能量与噪声功率谱密度之比Es/N0 是每符号能量与 N0 之比两者的差就是每符号的比特数EsN0_db EbN0_db 10*log10(k); % QPSK 下 k2差 3.01 dB对这个脚本里的单位能量星座Es 1Eb 0.5所以线性数值上 Es/N0 2 * Eb/N0dB 表示正好差 3 dB。很多人直接把 Eb/N0 的 dB 值当成符号信噪比传给awgn函数awgn内部按 SNR 处理结果整条曲线右移 3 dB理论和仿真对不上时很难看出来。还要分清「符号信噪比」和「比特信噪比」在误码率公式里的位置。理论 BER 公式用的是 Eb/N0理论 SER 公式用的是 Es/N0画图时横轴必须对应用哪个量再算理论值。混用横轴定义是曲线错位的头号原因。3.3 循环结构、随机数种子与 parfor 可复现性如果信噪比点比较多想用 parfor 并行加速这里有个常见的坑randn在 parfor worker 上的随机流和串行执行时不同每次运行结果会有差异虽然统计意义上曲线形状一致但具体数值不可复现。调试阶段建议先串行跑固定rng(0)确认代码没问题后再考虑并行。rng(0); % 放在所有随机数生成之前同一版本 MATLAB 下曲线可复现固定随机流的价值在于你发现 6 dB 这个点和理论偏差偏大时重新跑一次还是同一个值可以安心排查是不是统计波动而不是怀疑代码随机出错。换 MATLAB 版本或换机器后rand 的默认算法可能不同同样的rng(0)得到不同序列但曲线统计特征不变这不算 bug。4. 把仿真曲线和理论曲线叠在一起验证方法4.1 理论公式为什么是 0.5*erfc(sqrt(Eb/N0))QPSK 可以拆成两路正交的 BPSKI 路和 Q 路各解调各的比特互不干扰。每一路的噪声方差是 N0/2每路比特能量是 Eb所以单路误比特率就是标准 BPSK 的Q(sqrt(2*Eb/N0))。两路合起来总误比特率不变这就是 QPSK 和 BPSK 在相干解调下 BER 曲线重合的原因。MATLAB 里qfunc在通信工具箱里erfc是基础函数用erfc写更通用。两者的关系是qfunc(x) 0.5 * erfc(x/sqrt(2))代入 BPSK/QPSK 的 BER 公式就变成berTheory 0.5 * erfc(sqrt(10.^(EbNoVec/10))); % 线性的 Eb/N0这个公式的前提是格雷映射、理想相干解调、无脉冲成型、无符号间干扰。如果你的仿真加了根升余弦成型滤波和匹配滤波理论公式本身不变但噪声功率的定义要对齐到采样点否则理论曲线和仿真曲线会出现固定偏移。4.2 两条曲线怎么叠semilogy 绘图与图例BER 曲线横跨多个数量级线性坐标会压扁低误码率部分必须用 semilogy。绘图代码放在主循环之后figure; semilogy(EbNoVec, berSim, o-, MarkerFaceColor, auto); hold on; semilogy(EbNoVec, berTheory, r-, LineWidth, 1.2); grid on; legend(仿真 BER, 理论 BER, Location, southwest); xlabel(E_b/N_0 (dB)); ylabel(误比特率 BER);仿真点用空心圆加实线理论值用纯实线图例里明确写「仿真 BER」和「理论 BER」避免回头忘掉哪条是哪条。grid on对对数坐标尤其必要网格线横平竖直能快速看出仿真点相对理论线的偏移方向和大小。注意log 坐标的 Y 轴下限不要设成 0误码率不可能等于 0只是统计不到而已。Y 轴下限设 1e-7 或 1e-8 更合理。4.3 仿真点离理论线多远算正常统计波动区间曲线叠好后下一个问题是「差多少算有问题」。每个点的误码数 e errTotal相对统计偏差约±1/sqrt(e)到±2/sqrt(e)对应 68% 和 95% 置信区间。e100 时约 20% 的上下浮动都是正常的e10 时浮动可能超过 60%这个点基本没有参考价值。看曲线的正确姿势是低信噪比区0~2 dB误码率高e 大点应该贴理论线贴得很紧高信噪比区7~8 dB误码率低e 小点应该围着理论线抖。如果反过来了——低信噪比区就开始乱飘高信噪比区反而整整齐齐——那说明代码有系统性问题不是统计波动。这时候检查顺序是噪声功率定义、Eb/N0 换算、判决极性、误码统计口径按这个顺序排查比对着曲线猜快得多。5. 避坑与排查QPSK 误码率仿真最常见的 5 个翻车现场5.1 现象整条曲线左移约 3 dB —— 原因Eb/N0 被当成 Es/N0 用仿真的曲线和理论线形状完全一致但整体往左挪了大概 3 dB这是这类脚本里最典型的错位。原因是加噪声时用的是N0 1 / (10^(EbN0/10))这里把 Es1 当成了 Eb而 QPSK 下 Eb Es/2 0.5真实的 N0 应该更小等效信噪比被抬高了。解决一定要走「Es → Eb → N0」的换算链Eb mean(abs(tx).^2) / k; N0 Eb / (10^(EbN0/10)); % EbN0 是 dB 值顺带提醒如果后面改成 16QAMk4同样的问题就会变成 6 dB 的偏移。遵循「先算 Es再除 k 得 Eb」的写法换调制方式时不容易再踩一次。5.2 现象曲线在理论上方面积差约一倍 —— 原因把符号误码当比特误码仿真点比理论 BER 整体偏高而且高出不到一倍看起来像「有点问题但又不严重」。典型的错误是把每个符号判错计一次误码再除比特总数。一个 QPSK 符号判错时可能错 1 个比特也可能错 2 个比特格雷映射下绝大多数相邻错误只错 1 个比特但按符号计数会把「错 1 比特」和「错 2 比特」都算成 1 次错误再除以总比特数BER 被明显夸大。解决判决后按比特逐位对比得到的是真正的误比特数。或者单独算误符号率再用格雷映射的近似关系ber ≈ ser / k去核对。如果用的是symerr函数它的返回值是符号差错率不是比特差错率直接拿来当 BER 画图就是上面这个现象。5.3 现象高信噪比点掉到 0semilogy 画不出来跑到 7 dB、8 dB 时误码率算出来是 0log 坐标轴画不出 0 这个值曲线上出现一个断崖式的空洞。原因是统计到的误码数本来就是 0比如 1e6 比特在 8 dB 下期望误码只有几十个这轮运气好一个没碰上不是真的一点错都没有。解决给每个信噪比点设定最小误码数门槛误码数不够就继续加比特直到凑够再算或者设置固定的最低比特数并且明确标注这个点的统计下限是1/numBits。画图时把误码为 0 的点跳过或者人为改成 1/numBits 并画成空心标记不让曲线断掉但心里要清楚这个点只是下界。5.4 现象两次运行曲线完全不同 —— 原因随机数流未固定同一套代码连着跑两次某些点的 BER 差好几倍尤其是高信噪比区。原因是每次启动 MATLAB 时 rand 流默认状态不同噪声序列完全不一样。你以为代码有 bug其实只是统计波动被放大了。解决所有随机数生成之前加一行rng(0)或rng(default)。固定随机流后至少在同一环境里曲线完全可复现方便排查问题。如果要用 parfor 加速要意识到并行 worker 的随机流和串行时不同曲线形状一致但细节值对不上属于预期行为。5.5 现象程序报错「reshape 元素数目不匹配」—— 原因比特数不是 2 的整数倍numBits 999999这类数除以 2 有余数reshape(data, k, []).直接报错。还有一种更隐蔽的换 Eb/N0 数组长度时把numBits顺手改了但忘了保证它能被 k 整除。解决定义比特数时按符号数倒推numSym 5e5; numBits k * numSym;这样永远不会出现余数。如果代码里多处用到 numBits在reshape之前加一行断言assert(mod(numBits, k) 0, numBits 必须能被 k 整除);脚本开头加这个断言只花一行省的是改完参数后莫名奇妙报错的时间。6. 从误码率曲线往前走一步眼图与星座图验证根因BER 曲线告诉你「错了多少」但没告诉你「为什么错」。我习惯在跑完主循环后把中间变量的星座图直接画出来花十秒能看出很多问题figure; plot(rx, b.); % 加噪后的接收符号 hold on; plot(constellation, ro, MarkerSize, 8); axis equal; grid on; xlabel(I 路); ylabel(Q 路);正常的接收符号应该以四个星座点为中心形成四团圆形噪声云。如果云的中心不在星座点上检查判决前的归一化如果云是拉长的椭圆说明 I/Q 不平衡或噪声功率分配错了如果四团云边界模糊成一片说明 Eb/N0 太低曲线大概率吻合但点很密。眼图则需要过采样和脉冲成型才有意义用eyediagram(rxSignal, sps)能看到码间串扰和最佳采样时刻这是教科书仿真里很少真去画、工程上却很实用的一步。低误码率点跑不动的问题我的处理是蒙特卡洛只跑到 1e-5 量级更高信噪比点改用半分析法或直接标注「统计下限」。别为了凑一条完整的 1e-7 曲线让 MATLAB 跑一晚上这不划算。每次跑完先看 errTotal 到底统计到了几个误码少于 30 就说明这个点的可信度有限截图前心里先有数。这是我在被曲线坑过几次之后养成的习惯希望帮到你。本文还有配套的精品资源点击获取
返回列表