简介面向无线通信系统研究的一份 MATLAB 源码包重点实现脏纸编码和汤姆林森-哈拉希玛预编码两种高级干扰抑制技术。脏纸编码利用信道状态信息把已知干扰视为“脏纸”进行预补偿而汤姆林森-哈拉希玛预编码则在发射端通过线性处理降低接收端解码难度两者都能提升非理想信道下的传输可靠性但处理思路不同。资源包含两个 m 脚本文件分别用于基础取模运算和整体预编码流程演示覆盖信道建模、信号生成、预编码、传输、接收端解码以及误码率统计等关键步骤可直接运行并观察不同信噪比下的性能差异。整个压缩包只有两个文件均为 MATLAB 脚本总大小约 2KB代码精炼、易于阅读和修改适合通信专业学生、算法研究人员开展课程设计、算法验证或性能对比实验。目前已有 366 人浏览学习是快速理解脏纸编码与预编码实现细节的实用参考。1. 从已知干扰说起脏纸编码与THP的价值在下行链路里基站常常同时掌握发给多个用户的数据。如果其中一路数据对另一路用户来说就是干扰那就构成了典型的发射端已知干扰场景。脏纸编码Dirty Paper Coding, DPC和 Tomlinson-Harashima 预编码THP回答的是同一个问题既然干扰在发射端已知能不能通过预处理让接收端几乎感觉不到它的存在答案直接决定 MIMO 系统的吞吐量和 BER 表现。这里要拆解的是一套可直接运行的 MATLAB 示例包含modulo.m与Dirty_or_TH_precoding.m两个文件。看完你不仅能复现 BER 曲线还能把参数改到自己的链路上搞清楚 τ 怎么取、噪声折叠从哪来、以及为什么 DPC 在仿真里性能最好却最难落地。2. 原理辨析DPC 如何预见干扰THP 如何约束发射功率2.1 脏纸编码已知干扰不损失容量的信息论结论假设一张纸上已经写满杂乱笔迹而你恰好知道每一处笔迹的内容那么在同样的纸上重新书写新信息并不会比在白纸上写字更难。这就是著名的脏纸问题在加性高斯信道里只要干扰是发射端已知的信道容量就等于无干扰时的容量。这个结论对多用户 MIMO 下行链路意义重大——基站同时给多个用户发送数据每一路信号对其他用户而言都是干扰但理论上总容量并不会因这些互扰而损失。理解这个结论的关键在于已知二字。发射端如果不知道干扰就只能把它当作噪声去对抗如果知道就可以把干扰当作一种可逆的变换来处理。DPC 正是利用这个自由度在码本构造阶段就把干扰的影响补偿掉。代价是它隐含了一个假设发射端有足够的动态范围来处理干扰补偿后的信号。实际系统中这个动态范围受限于是就有了 THP 这种工程化的近似实现。2.2 THP 的模运算用折叠替代削波Tomlinson-Harashima 预编码并不追求完美消除干扰而是引入一个周期性的模运算。信号和干扰经过线性组合后如果幅度超出星座允许的范围就被折叠回[-τ/2, τ/2)区间内。这样发射端既做掉了干扰补偿又不至于因为产生超大功率信号而失真。接收端执行同样的模运算。由于折叠周期一致只要噪声不太大折叠操作在接收端可以完美还原原符号。这个过程类似差分编码中的相位折叠但作用在幅度域上。模运算的代价是引入了一个人为的非线性项——噪声落在边界附近时会被折叠到另一端这就是后文要讲的噪声折叠问题。2.3 三种预编码方案的取舍方案干扰补偿方式发射功率约束实现复杂度典型定位线性 ZF矩阵求逆直接抵消弱病态信道损失大低高 SNR 下的基线方案理想 DPC已知干扰预补偿无约束理论容量界高性能上界THP模运算 非线性反馈强约束在模边界内中MU-MISO/MU-MIMO 下行这三者经常出现在同一篇仿真文章的对比图里。ZF 是线性方法DPC 是理论上限THP 则是在两者之间取平衡的非线性方法。从仿真代码的角度看THP 和 DPC 的差别往往只在编码分支处的一行modulo(...)这也正是Dirty_or_TH_precoding.m把两者放进同一个脚本的原因。3. 代码拆解modulo.m 与 Dirty_or_TH_precoding.m 的关键实现3.1 用 floor 实现复数模运算modulo.mmodulo.m是整个 THP 方案的地基。它把输入信号折叠到以原点为中心的周期区间内实部和虚部分别处理。完整代码如下function u modulo(x, tau) % u modulo(x, tau) % 将输入 x 的实部与虚部分别折叠到 [-tau/2, tau/2) 区间。 % % x : 标量或向量/矩阵可为复数 % tau : 折叠周期必须为正实数 if isreal(x) u x - tau .* floor((x tau/2) ./ tau); else ur real(x) - tau .* floor((real(x) tau/2) ./ tau); ui imag(x) - tau .* floor((imag(x) tau/2) ./ tau); u complex(ur, ui); end end这里floor((x tau/2) / tau)判断当前值落在第几个折叠周期乘以tau得到需要减去的整数倍周期。实数分支直接用./算子保证输入是行向量时输出维度不变复数分支必须对实部虚部分开运算因为复数的取整没有定义。tau如果传入 0 或负数会直接触发除零错误所以调用前要做合法性检查常见做法是在主脚本参数区定义好之后再传入。提示这段代码在向量化场景下效率不错但如果你的仿真里符号数上百万可以考虑把floor换成mod并预先计算tau/2减少重复计算。3.2 主脚本结构参数区、信道区、编码区Dirty_or_TH_precoding.m按仿真链路顺序组织。参数区定义调制阶数、信噪比范围、每信道符号数和信道实现次数信道区生成瑞利衰落系数和已知干扰编码区通过method字段切换 DPC 与 THP最后是判决和 BER 统计。骨架如下clear; close all; clc; M 4; % QPSK SNRdB 0:2:18; N 2e4; % 每信道实例的符号数 nChan 30; % 独立信道实现次数 method THP; % THP / DPC sym [11i, 1-1i, -11i, -1-1i] / sqrt(2); tau 2 * sqrt(2); % 模边界见 4.2 节公式 for snrIdx 1:length(SNRdB) for ch 1:nChan h (randn 1i*randn) / sqrt(2); % 瑞利衰落 v (randn(1,N) 1i*randn(1,N)); % 已知干扰 s sym(randi(4, 1, N)); % 数据符号 % ---------- 编码分支 ---------- switch method case THP u modulo(s - v, tau); x u / h; case DPC u s - v; x u / h; end % ---------- 信道与噪声 ---------- nvar 1 / 10^(SNRdB(snrIdx)/10); n sqrt(nvar/2) * (randn(1,N) 1i*randn(1,N)); y h .* x v n; % ---------- 解码分支 ---------- switch method case THP z modulo(y, tau); case DPC z y; end % 最小距离判决与误码统计略 end endh是标量瑞利信道因子v是发射端已知的复干扰序列功率与信号的量级一致。注意这里没有对x做平均功率归一化原因在 5.2 节展开。nvar 1/SNR的写法把符号能量固定为 1所以噪声方差直接由信噪比倒数给出省去 Es/N0 换算但改成 Eb/N0 时要除以每符号比特数。3.3 THP 分支的数学验证与接收端模运算THP 分支的核心只有三行发射端u modulo(s - v, tau)信道后y h*x v n接收端z modulo(y, tau)。验证过程如下设u modulo(s - v, tau)根据模运算定义有u s - v k*tau其中k是某个整数实部虚部各自成立。发射信号x u / h接收端去掉信道后y h * (u / h) v n u v n (s - v k*tau) v n s k*tau n对y再做一次模运算modulo(s k*tau n, tau) modulo(s n, tau)。只要s n落在(-tau/2, tau/2)区间内结果就是s n干扰被完全抵消。这里必须强调顺序发射端在s - v上做模接收端在y上做模两次折叠周期相同。如果在发射端做了模、接收端却直接判决那k*tau项就会残留BER 会高到完全不可用这是初学者最容易犯的错误。3.4 DPC 分支去掉模运算之后发生了什么DPC 分支把编码改成了u s - v省掉了modulo。这意味着发射端要发射x (s - v) / h来完全消除干扰没有幅度限制。从数学上接收端y s n干扰完美去除BER 曲线看起来会非常漂亮。但实际系统里这个分支几乎不可用当信道h处于深衰落幅度接近 0或者干扰v的瞬时幅度很大时x的发射功率会飙升。你可以在代码里加一行观察avg_power mean(abs(x).^2); fprintf(method%s, 平均发射功率 %.2f dB\n, method, 10*log10(avg_power));同样信道条件下THP 分支的平均发射功率通常在 0 dB 附近波动DPC 分支则可能高出 10 dB 以上。这解释了为什么 THP 被称为可实现的脏纸编码——它用模运算换掉了无限功率需求代价是噪声折叠带来的少量性能损失。4. MATLAB 无线通信实战从 QPSK 到 16-QAM 的参数调整4.1 信道与干扰的生成方式瑞利信道用两个独立高斯变量合成h (randn 1i*randn)/sqrt(2)。除以sqrt(2)是为了让|h|^2的期望为 1这样 SNR 定义不会受到信道平均增益偏移。如果你要模拟莱斯信道或带多普勒的时变信道可以替换为comm.RicianChannel但块衰落下的研究结论基本一致。已知干扰v在这个模型里被设计为与信号同量级的复高斯随机序列。更贴近实际的做法是让v来自另一个用户的数据符号乘上某个系数比如v 0.8 * s_other这样可以观察不同干扰强度对 BER 曲线的影响。修改这个系数后你会发现 THP 的 BER 几乎不变而普通无预编码链路的 BER 会明显上升这正是模运算抗已知干扰的直接体现。4.2 τ 到底怎么取从星座边界反推模边界tau不是拍脑袋定的。它必须满足两个条件一是能包住全部星座点二是给噪声留出余量。常见的推导公式是tau 2 * (max(|Re(sym)|) d_min/2)其中d_min是最小星座点距。代码如下dist abs(sym(:) - sym(:).); dist dist eye(length(sym)) * inf; % 排除对角线上的 0 d_min min(dist(:)); tau 2 * (max(abs(real(sym))) d_min/2);对于 QPSKmax|Re(sym)| 1/sqrt(2)d_min sqrt(2)算出tau 2*sqrt(2)。对于 16-QAM星座点归一化到平均功率 1 后max|Re(sym)| 3/sqrt(10)d_min 2/sqrt(10)于是tau 2 * (3/sqrt(10) 1/sqrt(10)) 8/sqrt(10) ≈ 2.53调制方式平均符号功率max|Re(sym)|d_mintauQPSK10.7071.4142.82816-QAM10.9490.6322.530注意这里的d_min是针对功率归一化的星座算出来的。如果你用qammod(0:M-1, M)生成星座它的平均功率不是 1那么tau 2*(31) 8与表格里的数值差了一个sqrt(10)因子原理相同但 SNR 定义要跟着改。4.3 扩展到 16-QAM 的改法把M从 4 改成 16 不是只换星座表那么简单。第一符号能量变了通常在星座映射后乘1/sqrt(10)归一化第二比特映射建议用 Gray 码否则相邻符号的误码会引起多个比特错误第三判决部分如果继续用最小距离判决要注意复数距离矩阵的维度。一个干净的写法是if M 4 sym [11i, 1-1i, -11i, -1-1i] / sqrt(2); elseif M 16 pts [-3, -1, 1, 3] / sqrt(10); [I, Q] meshgrid(pts); sym I(:) 1i*Q(:); end放大N到 2e4、增加nChan到 50BER 曲线的毛刺就会明显减少。如果想把曲线画得更平滑还可以按每 SNR 点累积固定错误数的方式停止仿真比如累计到 200 个错误才跳出循环这样可以避免高 SNR 下置信区间过宽。4.4 读懂第一张 BER 曲线运行脚本后会看到一条随 SNR 下降的 BER 曲线。THP 分支的曲线在高 SNR 时斜率与理论 QPSK 在 AWGN 下的曲线基本平行但整体右移大约 1~2 dBDPC 分支的曲线则非常接近 AWGN 理论值。如果你把method改成ZF发射端用x s/h不做干扰抵消BER 会直接平在 0.5 附近——因为没有消除干扰这正是对照实验的意义。5. 模边界、噪声折叠与 CSI 误差THP 仿真中的四个坑5.1 τ 太大浪费功率τ 太小出现混叠把tau调大模运算的窗口变宽发射信号允许的幅度范围变大THP 退化成线性预编码功率优势消失把tau调小星座点可能落在折叠边界之外接收端模运算会把信号折到错误的象限误码率直接飙升。一个实用的验证方法是画出u的实部直方图如果直方图在±tau/2处出现明显的截断堆积说明 τ 选小了如果分布集中在中部且几乎没有边界附近的样本说明 τ 偏大可以适当缩小。5.2 功率归一化的两个流派与一个陷阱常见做法是在发射端对x做平均功率归一化x x / sqrt(mean(abs(x).^2))。这在一般通信仿真里没问题但放到 THP 中要格外小心归一化因子会等比例缩放整个信号接收端如果再直接做modulo(y, tau)模边界与信号幅度不再匹配折叠关系被破坏。我一般建议在纯预编码对比仿真里不做平均功率归一化而是让模运算自然约束功率。如果一定要归一化必须把缩放因子一直传送到接收端先除以缩放因子再做模运算否则不会有预期效果。5.3 接收端模运算带来的噪声折叠接收端的modulo会把落在边界附近的噪声折叠回(-tau/2, tau/2)内。直观地说如果噪声把某个符号推过tau/2边界接收端取模后它会被折回到负半轴即使离正确星座点只差一小段距离也可能判到完全不同的符号上。这就是 THP 相对 DPC 损失 1~2 dB 的主要原因。噪声折叠无法消除但可以缓解增大 τ 让折叠发生得更不频繁同时代价是发射功率上升。实践中最优 τ 往往略大于 4.2 节的理论最小值具体做法是扫描tau * [0.8, 0.9, 1.0, 1.1, 1.2]找 BER 最低点。5.4 信道估计误差对两种策略的差异化影响THP 的反馈与信道分解绑定CSI 误差会在反馈项里逐级传播所以第二路及之后的数据流对 CSI 误差更敏感。而 DPC 分支的干扰补偿是直接相减误差表现为残余干扰的线性叠加而不是非线性折叠。仿真中可以给信道加一个乘性误差h_hat h 0.05 * (randn 1i*randn) / sqrt(2); x u / h_hat; % 用带误差的估值做预均衡你会观察到 THP 的 BER 曲线出现平台而 DPC 的曲线只是整体抬高。这个差异在设计导频和信道估计算法时很关键——如果系统 CSI 质量一般THP 的优势会比理想 CSI 下变小。6. 把预编码代码推进到 MU-MIMO 与 Simulink 联调6.1 从 SISO 到 2×2 MU-MISO 的扩展框架单链路演示能讲清原理但 THP 真正的舞台是多用户 MISO/MIMO 下行。扩展思路是把标量反馈换成向量反馈信道矩阵H用户数 × 天线数用 LQ 分解得到反馈矩阵和酉矩阵。核心代码只有几行[Qm, Rm] qr(H, 0); L Rm; % 下三角 F L * diag(1 ./ diag(L)); % 单位对角下三角 P Qm * diag(1 ./ diag(L)); % 发送预编码矩阵 u1 s1; u2 modulo(s2 - F(2,1) * u1, tau); x P * [u1; u2]; % 2 天线发送向量接收端各自做缩放和模运算即可。这个框架下你可以直接比较线性 ZF、THP 与理想 DPC 三者的和速率和 BER也是目前 MU-MIMO 下行链路论文里最常见的一组对照曲线。6.2 与 MATLAB 优化工具箱、Simulink 的对接如果你习惯用 matlab 优化工具箱搜索预编码矩阵或功率分配参数可以把编码分支整理成一个函数句柄交给fmincon或者ga做联合优化。而在系统级联调场景中modulo.m可以封装成 MATLAB Function Block 直接放进 Simulink注意把tau设为模块参数而不是硬编码在函数内部这样不同调制方式可以直接复用同一个模块。大规模蒙特卡洛循环建议把内层for ch 1:nChan改成parfor同时用RandStream管理随机数流保证并行运算的结果可以复现。否则每次跑出来的曲线毛刺位置都不一样排错会很痛苦。6.3 用 bertool 快速核对仿真曲线调试 BER 仿真时先用理论曲线做锚点能省很多时间。在 MATLAB 命令窗口运行bertool选择Monte Carlo页签把 4.2 节生成的 SNR 点与 BER 结果填入再叠加 QPSK 在 AWGN 下的理论曲线。如果 THP 仿真曲线和理论曲线随 SNR 的下降斜率一致说明链路基本正确如果斜率偏差较大优先检查是否是 5.1 节的 τ 选取问题。这个方法同样适用于把method切成DPC时验证理想干扰消除后 BER 应接近 AWGN 理论值这一结论。本文还有配套的精品资源点击获取