多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

16QAM软解调与LDPC编译码及FFT频偏估计的MATLAB仿真实现

16QAM软解调与LDPC编译码及FFT频偏估计的MATLAB仿真实现 简介基于Matlab2024b环境面向通信工程、电子信息等领域提供一套结合16QAM调制软解调、LDPC编译码与FFT频偏估计的通信系统误码率仿真资源适合课程设计与科研预研帮助需要研究误码率性能的学生和工程师快速理解完整链路。整套资源共14个文件包括9个M脚本、4个MAT数据文件和1个文本说明压缩包仅112KB这些M脚本覆盖主函数、LDPC编译码、调制解调与频偏估计等模块并附中文注释与程序操作视频。仿真链路完整覆盖随机二进制序列生成、LDPC编码、16QAM调制、AWGN信道、FFT频偏估计与补偿、软解调、LDPC译码及误码率统计等环节可直接运行主程序脚本复现结果。另有对比程序便于分析不同模块对误码率的影响操作视频还演示了MATLAB路径设置等关键步骤。目前已有70人学习适合需要系统掌握通信同步与编译码仿真细节的进阶用户。1. 基于16QAM软解调LDPC编译码FFT频偏估计的通信系统仿真这份MATLAB资源解决什么问题做16QAM软解调LDPC编译码FFT频偏估计这条完整链路时最容易卡住的地方不是单个模块而是几个模块拼起来之后BER曲线永远比理论门限差几个dB。这份MATLAB误码率仿真资源把发射端、加噪信道、FFT频偏估计与补偿、16QAM软解调、LDPC译码和误码率统计完整串成一条可运行的链路自带中文注释和程序操作视频。它解决的是三类具体问题LDPC软判决译码需要的LLR怎么从16QAM符号里算出来FFT频偏估计的峰值检索和补偿相位怎么跟帧结构配套以及Eb/N0与Es/N0换算错误导致的整条曲线偏移。适合通信方向做课程设计、毕业设计或者打算把基带仿真链路快速跑通的工程师。2. 链路架构与模块选型16QAM软解调、LDPC译码器和FFT频偏估计的配合关系2.1 发射端数据通路从信息比特到16QAM符号这条链路的源头是随机信息比特经过LDPC编码变成带冗余的码字比特再按4个比特一组映射成16QAM符号。为什么中间要插一个bit2int步骤因为qammod函数只接受0到15的整数索引不接受比特向量必须先做比特到符号索引的转换。常见做法是让LDPC编码输出列向量然后用bit2int(codeword, 4)把每4比特并成一个整数这样qammod才能直接吃进去。16QAM符号映射选用Gray映射而不是自然二进制映射。Gray映射保证相邻星座点之间只有1个比特不同硬判决时单符号错误只造成1个比特错误这对后面LDPC译码的迭代收敛帮助明显。MATLAB里qammod(0:15, 16, gray, UnitAveragePower, true)一行就能生成归一化星座点注意这里必须加UnitAveragePower否则星座点平均功率是10而不是1后面噪声方差全部要重新换算这是一个非常隐蔽的坑。发射端还要考虑前导训练序列。FFT频偏估计不能拿随机数据硬算帧结构里要放一段收发两端都知道的训练符号。我一般会在每一帧开头拼接128个固定训练符号接收端拿本地训练序列和接收前导做共轭相关才能把调制信息消掉留下纯粹的频偏相位旋转项。整条数据通路可以归纳成下面这张表模块输入输出关键参数随机数发生器无信息比特列向量长度32400均匀分布LDPC编码器信息比特码字比特DVB-S.2码率1/2码长64800bit2int码字比特符号索引0~15每4比特一组MSB优先16QAM调制符号索引复基带符号Gray映射平均功率1帧拼接训练序列数据完整发送帧训练符号128个参数选择上码长64800的DVB-S.2 LDPC码是卫星通信标准的经典配置纠错性能好但仿真耗时高。如果只是验证链路正确性第一次跑建议把帧循环次数设小一点先跑到能画出BER曲线趋势再逐步加大。训练序列长度128个符号对FFT频偏估计来说已经够用频率分辨率取决于训练序列长度和FFT补零点数这个在第5章会展开讲。2.2 为什么接收端要用软解调而不是硬判决接收端从16QAM符号恢复比特有两种做法。硬判决是算出每个符号离哪个星座点最近直接判成对应4比特然后丢给LDPC译码器。软判决则是计算每个比特的对数似然比LLRLLR的正负号代表这个比特倾向0还是1幅度大小代表置信度。LDPC的置信传播译码算法本质上是概率域的迭代喂给它带置信度的软信息收敛速度和误码性能都明显优于硬判决。在工程仿真里软解调带来的增益通常在1到2dB左右这在LDPC纠错链路上是质变。如果接收端用硬判决LDPC译码器拿到的是0/1硬值相当于把概率信息全部丢掉了迭代译码的时候外信息更新只能靠翻转性能大打折扣。16QAM每个符号携带4个比特这4个比特的置信度并不相同实部和虚部的符号位决定象限幅度位决定象限内的精细位置幅度位的LLR在星座中心附近天然偏小这个差异只有软判决才能体现出来。具体实现上LLR计算可以走两条路。一条是精确的对数似然计算要遍历16个星座点做log-sum-exp运算Matlab写起来直观但仿真慢另一条是max-log近似只找离接收符号最近的0比特星座点和最近的1比特星座点用两者距离差除以噪声方差作为LLR。通信系统仿真里95%的工程实现都用max-log性能损失在0.1dB以内代价是代码简单好几倍。后面的3.1节会给出完整函数。2.3 FFT频偏估计在整个同步链中的位置频偏的来源很现实收发晶振不可能完全同频终端移动还有多普勒反映到基带就是接收符号附加了一个随符号索引线性增长的相位旋转。如果不做补偿16QAM星座图会呈圆形旋转软解调算出来的LLR全是错的LDPC译码器再强也白搭。同步链路的常见处理顺序是先做粗定时找到帧起点再做频偏估计最后开始解调译码。FFT频偏估计的基本思路是利用训练序列的周期重复结构。接收前导和本地参考前导共轭相乘之后调制信息被消掉剩下exp(j2πΔf n)的复指数序列对它做FFT频谱峰值所在频率就是频偏。这个方法的优势是不需要像锁相环那样迭代收敛一次FFT就能得到整个搜索范围内的粗估结果特别适合突发通信机制。频偏估计放在软解调之前还有个工程原因LLR计算依赖星座点位置任何残余相位旋转都会让距离差计算失真。FFT估计出来的频偏先补偿掉大部分剩余的小残差对LLR的影响就有限了。但需要注意FFT频偏估计有一个天然的分辨率限制fEst的精度大致等于fs除以FFT点数补零只能改善显示效果不能创造真实分辨率这个误区和改进办法在5.3节专门讲。3. 核心模块MATLAB实现LLR计算、LDPC译码器配置与FFT频偏估计3.1 16QAM软解调基于Gray映射的LLR函数软解调函数的核心输入是接收符号列向量和噪声方差输出是4×N的LLR矩阵每行对应一个比特位置。先通过qammod生成Gray映射的16星座点再用de2bi生成0到15索引对应的4比特表最后遍历4个比特位对每个比特找最近的0集合星座点和最近的1集合星座点距离差除以噪声方差得到LLR。function llr softDemod16QAM(rxSym, noiseVar) % 16QAM软解调max-log LLR近似 % rxSym : 接收符号列向量 % noiseVar: 复高斯噪声总方差等于10^(-EsN0dB/10) const qammod(0:15, 16, gray, UnitAveragePower, true); bits de2bi(0:15, 4, left-msb); llr zeros(4, numel(rxSym)); for k 1:4 idx0 find(bits(:, k) 0); idx1 find(bits(:, k) 1); d0 min(abs(rxSym. - const(idx0).).^2, [], 1); d1 min(abs(rxSym. - const(idx1).).^2, [], 1); llr(k, :) (d0 - d1) / noiseVar; end end这段代码里最关键的是量纲匹配。qammod加UnitAveragePower之后星座点平均功率为1所以添加噪声时复高斯总方差noiseVar就是每符号的噪声能量。距离差在星座空间里是能量单位除以噪声方差后才是无量纲的LLR直接喂给LDPC译码器才能被正确解释。很多翻车案例就是把noiseVar写成了N0而不是双边带方差或者把星座点功率当成1实际是10导致LLR整体缩放错误。de2bi用left-msb表示每行4比特里左边是最高位这个顺序必须和发射端bit2int的MSB优先保持一致。如果两端一个用MSB一个用LSB整个比特序会错乱LDPC译码后误码率会卡在0.5附近而且不随SNR下降这是典型的黑匣子问题。3.2 LDPC编码器/译码器配置DVB-S.2 1/2码率为例LDPC模块直接依赖Communications Toolbox里的comm.LDPCEncoder和comm.LDPCDecoder校验矩阵用dvbs2ldpc生成。这段代码只负责配置对象真正的编译码在每帧循环里调用。H dvbs2ldpc(1/2); % DVB-S.2标准码率1/2码长64800 ldpcEncoder comm.LDPCEncoder(H); ldpcDecoder comm.LDPCDecoder(H, ... DecisionMethod, Soft decision, ... MaximumIterationCount, 50, ... OutputValue, Whole codeword);为什么选DVB-S.2的校验矩阵因为dvbs2ldpc是MATLAB内置函数不需要额外下载校验矩阵数据而且DVB-S.2的LDPC是系统码结构编码输出前32400位就是原始信息位译码后直接取前32400位对比即可不需要额外提取信息位的下标。MaximumIterationCount设为50是性能和耗时的折中迭代太多次数低SNR下收敛慢但每帧译码时间成倍上升迭代太少10次以内往往在错误平台上跑不出来BER曲线会异常陡峭甚至不收敛。译码器配置里的Soft decision必须和3.1节的LLR输出配套。如果写成Hard decision输入需要是0/1硬值与LLR完全不兼容误码率会直接飙升。OutputValue选Whole codeword是为了方便对比信息位如果选Information part输出长度就是32400也能用但不同MATLAB版本的默认行为略有差异显式指定更稳妥。这里有一个耗时敏感的提醒码长64800的LDPC每译一帧要做50次迭代MATLAB单帧耗时可能零点几秒到几秒不等。整条BER曲线如果跑17个SNR点、每点20帧总共340帧低SNR段会明显感觉慢。首次跑通链路时建议把numFrames先改为2确认无误再拉大量帧。3.3 FFT频偏估计与补偿峰值检索与相位旋转频偏估计函数利用前导的共轭相关去掉调制信息然后对相关序列做FFT并找峰值。实现上有三个细节要注意FFT点数要适当补零fftshift把零频移到中间频率轴用归一化符号速率fs来标定。function fEst fftFreqEst(preambleRx, preambleRef, fs) % 基于FFT的频偏估计 % preambleRx : 接收端截取的前导符号 % preambleRef: 本地已知训练符号 % fs : 符号速率基带仿真里通常取1 corrSeq preambleRx .* conj(preambleRef); % 去调制剩exp(j2pifn) Nfft 2^nextpow2(length(corrSeq) * 4); % 补零到4倍做插值 spec fftshift(fft(corrSeq, Nfft)); [~, idx] max(abs(spec)); fBins ((0:Nfft-1) - Nfft/2) / Nfft * fs; % 频率轴单位Hz fEst fBins(idx); end共轭相关的物理含义是接收前导符号取共轭后乘以本地参考星座调制相位被抵消输出序列只剩下频偏引入的旋转e的j2πΔf n次方以及残留噪声。对这样一条纯复指数序列做FFT频谱峰值的位置就是Δf。Nfft取4倍长度补零属于频率插值它让峰值看起来更精细但真实分辨率仍然由corrSeq的原始长度决定这一点必须心里有数。补偿时要用累计相位不能只补偿一个固定相位。每个符号经历的相位旋转是2π乘以fEst再乘以符号序号所以补偿指数要逐符号变化n (0:length(rxData)-1).; rxDataComp rxData .* exp(-1j * 2 * pi * fEst / fs * n);参数fs在这里是符号速率不是过采样率。如果仿真链路里做了成型滤波和过采样那么FFT频偏估计的输入应该是过采样后的符号fs也要相应换成采样率否则频率轴标定全部偏差。这个细节在基带仿真里最容易忽略我见过不少人在过采样链路里把fs直接取1结果频偏估计值比真实值大好几倍。4. 完整误码率仿真主程序参数设置、蒙特卡洛循环与BER曲线绘制4.1 仿真参数表Eb/N0范围、帧数、码长怎么定仿真参数要把调制阶数、码率、帧长、噪声换算一次性定清楚否则后面每改一个参数都容易出连锁错误。这张表是我常用的起点配置可以直接照抄参数取值说明调制方式16QAM每符号4比特Gray映射LDPC码率1/2DVB-S.2码长64800信息比特长度32400码字比特长度64800每符号信息比特log2(16) × 1/2 2Eb换算Es时的核心系数Eb/N0范围0~8 dB步进0.5 dB覆盖LDPC门限附近每SNR点帧数20首跑建议2帧起步训练符号数128固定序列用于FFT频偏估计归一化频偏0.01每个符号旋转0.01个周期每符号信息比特这个值很关键。16QAM每个符号原本携带4个码字比特码率1/2之后每个符号实际只携带2个信息比特所以Es/N0和Eb/N0相差10×log10(4×0.5)3.0103dB。如果把码率因子漏掉整个BER曲线会右移3dB看起来像LDPC没工作似的。4.2 主循环代码结构与关键行说明主程序的核心是双层循环外层扫Eb/N0内层做蒙特卡洛帧累加。每个SNR点先把Eb/N0换算成Es/N0再换算成总噪声方差然后逐帧走完编码、调制、加频偏、加噪声、频偏估计补偿、软解调、LDPC译码、误比特统计的完整流程。EbN0dB 0:0.5:8; numFrames 20; infoLen 32400; bitsPerSym 4; codeRate 1/2; ber zeros(size(EbN0dB)); trainIdx mod((0:127).^2, 16); % 固定训练序列收发已知 trainSym qammod(trainIdx, 16, gray, UnitAveragePower, true); for snrIdx 1:numel(EbN0dB) EsN0dB EbN0dB(snrIdx) 10*log10(bitsPerSym * codeRate); noiseVar 10^(-EsN0dB/10); errBits 0; totalBits 0; for fr 1:numFrames data randi([0 1], infoLen, 1); codeword ldpcEncoder(data); symIdx bit2int(codeword, bitsPerSym); txSym qammod(symIdx, 16, gray, UnitAveragePower, true); frame [trainSym; txSym]; fOffset 0.01; % 归一化频偏 nFrame (0:length(frame)-1).; rxFrame frame .* exp(1j*2*pi*fOffset*nFrame); noise sqrt(noiseVar/2) * (randn(size(rxFrame)) ... 1j*randn(size(rxFrame))); rxFrame rxFrame noise; fEst fftFreqEst(rxFrame(1:128), trainSym, 1); rxFrameComp rxFrame .* exp(-1j*2*pi*fEst*nFrame); rxData rxFrameComp(129:end); llr softDemod16QAM(rxData, noiseVar); decoded ldpcDecoder(llr(:)); errBits errBits sum(data ~ decoded(1:infoLen)); totalBits totalBits infoLen; end ber(snrIdx) errBits / totalBits; end semilogy(EbN0dB, ber, o-); grid on; xlabel(Eb/N0 (dB)); ylabel(BER);噪声生成那行需要解释清楚。noiseVar是复噪声总方差包含实部和虚部两个维度每一维的方差是noiseVar/2所以生成噪声时要乘以sqrt(noiseVar/2)。这里如果用sqrt(noiseVar)等效噪声功率会翻倍等效SNR下降3dB跟漏掉码率因子的效果叠加起来就是6dB偏差非常致命。频偏估计用的是rxFrame前128个训练符号估计出fEst之后对整个帧做相位补偿。nFrame从0开始累计与发送端加频偏的相位完全对齐。这里有个细节接收端没做定时同步直接取前128个符号当作训练序列前提是仿真里帧起点已知。真实系统中必须在帧同步之后才能用这个方案5.4节会展开讲相关坑。4.3 误码率统计口径信息比特还是码字比特误码率统计必须用信息比特不能统计码字比特。原因是Eb的定义就是信息比特能量如果用码字比特数做分母曲线和理论参照不在一个坐标系下。比如同样1000个信息比特码率1/2对应2000个码字比特分母翻倍之后BER数值会偏小一半曲线看起来变好了实际上是统计口径错误。我建议每帧统计完误比特数后用errBits/totalBits的方式累计而不是每帧单独算一个BER再平均。单独平均会给低误码帧过高权重导致曲线抖动。每SNR点至少积累几十个误比特再停止这个也可以用while循环实现比如while errBits 100这样低SNR段自动多跑几帧高SNR段误码较少时也不会跑太久。首跑阶段可以不管统计精度先把链路跑通看趋势。5. 仿真避坑实录LDPC16QAM链路最常见的五个坑5.1 误码曲线整体右移Eb/N0与Es/N0换算写错现象是整条BER曲线往右偏大约3dBLDPC在某个理论门限附近本该急剧下降的地方还在高位徘徊。原因就是直接把Eb/N0当成了Es/N0来算噪声方差漏掉了调制阶数4和码率1/2这两个因子。解决方法是先统一换算EsN0dB EbN0dB 10×log10(bitsPerSym × codeRate)。我自己的习惯是把每符号信息比特数写成独立变量这样换QPSK或者码率2/3时不用回去改公式。5.2 LLR量纲不匹配噪声方差或星座功率写错现象是BER曲线下降趋势对但比理论性能差1到2dB而且降到某个平台后不再往下走。原因多半是qammod没加UnitAveragePower星座点平均功率是10而不是1距离差算出来偏大除以一个偏小的noiseVarLLR整体偏大或者噪声生成时用了sqrt(noiseVar)而不是sqrt(noiseVar/2)。解决方法是发射端qammod固定加UnitAveragePower噪声生成用sqrt(noiseVar/2)然后放一个固定无噪声单帧检查LLR的符号和幅度分布是否合理。这个坑最隐蔽因为星座图看起来完全正常只有译码性能悄悄变差。5.3 FFT频偏估计峰值找不正补零不等于提高分辨率现象是频偏估计值总是和真实值差一两个频率bin而且误差不随SNR增大而消失。原因是补零只做插值不带来新信息训练序列长度只有128时FFT真实分辨率是fs/128即使补到512点主瓣宽度还是由128决定峰值位置量化误差就在半个bin左右。解决方法是先保证频偏大于fs除以4倍训练序列长度再用抛物线插值细化峰值。插值方式是在峰值附近取三个点拟合二次曲线然后找曲线顶点代码很简单mag abs(spec); p idx; alpha mag(p-1) - 2*mag(p) mag(p1); delta 0.5 * (mag(p-1) - mag(p1)) / alpha; fEst fBins(p) delta * fs / Nfft;加了插值之后估计精度通常能提升一个数量级代价是几行代码。频谱泄露对这个估计影响也很大如果训练序列不是整周期截断主瓣会展宽峰值偏移严重时甚至漏检所以训练符号序列设计要避免周期性与FFT点数高度相关。5.4 频偏补偿方向搞反相位旋转越补越歪现象是补偿后星座图不但没聚拢反而旋转得更严重BER比不补偿还差。原因是补偿相位符号写反了发送端加的是exp(j2πfOffset·n)补偿端必须用exp(-j2πfEst·n)一正一负才是逆运算。这个错误在代码里特别容易粗心因为MATLAB不报错只有星座图能看出来。解决方法是写完补偿代码后立刻画scatterplot对比补偿前后的星座图如果旋转方向反了会非常直观。从那以后我每次都在补偿之后检查一次前导段的残余相位前导星座点应该贴在标准16QAM位置上。5.5 LDPC迭代次数或帧数太少BER曲线毛刺和断崖现象是BER曲线在低SNR段上下跳动高SNR段突然变成零画不出semilogy。原因是MaximumIterationCount设得太小比如10次译码在部分帧上没收敛同时每SNR点只跑5帧统计样本不够低SNR段误码率方差极大。解决方法是迭代次数保底50次每SNR点至少积累50到100个误比特再停止。高SNR段误码为零本身是好事但显示上要人为给一个下限比如bermax(ber,1e-6)否则semilogy会报警告。首跑调链路时故意把SNR范围放宽到0到10dB先看整体趋势再逐步收缩。6. 进阶验证FFT频偏估计精度、BER门限自检与扩展玩法6.1 用单帧无频偏链路验证代码正确性拿到程序后第一件事不是跑整条BER曲线而是关掉频偏和噪声先验证发射机和接收机闭环正确。把fOffset设为0噪声方差设一个很小的值比如1e-6跑一帧检查解码比特和原始信息比特是否完全一致。如果这一帧都错问题基本出在bit2int的比特序、LDPC配置、或者LLR量纲上跟频偏估计无关。这个单帧测试是整个仿真链路的后悔药能在5分钟内定位大多数低级错误。6.2 FFT频偏估计的MSE验证把频偏估计单独抽出来验证精度方法是固定SNR在10dB让真实频偏从0.001到0.02等间隔变化每个频偏值跑几十次统计估计值的均方误差。你会发现误差呈现出阶梯状这是FFT栅栏效应导致的加上抛物线插值之后误差会平滑很多。注意频偏不能超过±0.5因为FFT频谱峰值检索在±fs/2范围内是单值的超过这个范围会出现折叠模糊估计值会跳到镜像位置。6.3 扩展换调制阶数、换码率、换LDPC标准这套架构最值钱的地方在于模块化。想换QPSK只需要写一个2比特的软解调函数或者直接复用16QAM的代码改成max-log二维判决想换码率把dvbs2ldpc(1/2)换成dvbs2ldpc(3/4)同时把infoLen改成对应的N×(1-rate)想换5G NR的LDPCMATLAB里ldpcQuasiCyclicMatrix可以根据基矩阵生成校验矩阵编码器接口完全一样。每次改版之后都重新走一遍6.1的单帧闭环验证再看BER趋势是否往理论方向移动。有一次我为了省事直接在一个已有链路上改码率忘了同步更新infoLen结果编码器报错排查了半天才发现是变量没跟着码率变。从那以后我每次换参数都强制走一遍单帧闭环验证再跑BER曲线这个习惯帮我省掉了大量翻车时间。做通信仿真慢一点没关系每一步都验证到位最终曲线才值得信。希望帮到你。本文还有配套的精品资源点击获取
返回列表