做物理层仿真的人这两年对OTFS应该不陌生了。高铁、低轨卫星、无人机这类高动态场景把通信系统的多普勒容限逼到了极限OFDM在时变信道下的子载波正交性被破坏误码率惨不忍睹。OTFSOrthogonal Time Frequency Space正交时频空间调制把符号直接调制在时延-多普勒域上信道矩阵变得稀疏且近似时不变一下子打开了局面。不过问题也来了单天线OTFS还好一旦升级到MIMO-OTFS接收端的信号检测复杂度就像滚雪球一样涨上去。MP检测算法消息传递检测、ZF均衡、低复杂度LU分解、误差纠正这几样东西组合起来才是工程上能把MIMO-OTFS接收机真正跑起来的关键。这篇文章把我最近一次MIMO-OTFS链路级仿真里做的算法选型、实现细节和踩坑过程完整写出来特别是MP算法那块带详细注释的代码以及LU分解怎么把矩阵求逆的复杂度压下来给正在做OTFS物理层仿真或者准备入坑6G候选波形的研究生和预研工程师一个可复现的参考。1. 为什么高速移动场景下OFDM会“掉链子”而OTFS不会1.1 OFDM的多普勒脆弱性与ICI问题OFDM之所以成为4G/5G的基石靠的是子载波间的正交性。接收端做FFT的时候只要子载波频率间隔Δf远大于信道相干带宽每个子载波上的衰落就是平坦的均衡只需要一个复数除法。但这一切都建立在信道在一个OFDM符号周期内近似不变的前提上。一旦移动速度提上来多普勒频移f_d v·f_c/c不再可以忽略。举个例子载频3.5GHz高铁速度360km/h多普勒频移大约是1167Hz。如果子载波间隔是15kHz这个多普勒只占子载波间隔的7.8%左右还在容忍范围内。但要是到了低轨卫星场景多普勒频移能到几十kHz甚至上百kHz相对子载波间隔的比例可能超过0.5这个时候子载波之间的正交性就彻底保不住了。数学上看OFDM接收端第k个子载波上的信号是y_k H_{k,k} x_k ∑_{l≠k} H_{k,l} x_l w_k那个∑项就是载波间干扰ICI。信道时变越严重H矩阵越偏离对角阵ICI越强。传统的办法是做频域均衡把H当成一个满阵去处理但OFDM的子载波数往往是1024、2048这种量级对满矩阵求逆的复杂度是天文数字。所以实际系统只能在高速场景下大幅度牺牲性能或者引入复杂的ICI消除算法效果还不稳定。1.2 OTFS把信道“赶”到了时延-多普勒域OTFS的思路换个角度既然时频域的信道时变特性难处理那能不能找一个新的域让信道在这个域里是稀疏的、准静态的这个域就是时延-多普勒域Delay-DopplerDD域。物理上想想很直观。电波传播经过P条路径每条路径有一个特定的时延τ_i和多普勒频移ν_i。在DD域看每条路径就是平面上的一个点τ_i, ν_i信道响应本身几乎不随时间变化——变化的只是每条路径的相位。所以DD域的信道矩阵本质上是由P个抽头决定的稀疏度极高。OTFS发射端的工作流程可以这样理解先把QAM符号放在一个N×M的DD域网格上X_dd[n, m]里的n对应多普勒轴m对应时延轴。然后通过逆辛有限傅里叶变换ISFFT把二维数据块映射到时频域X_tf F_N · X_dd · F_M^H这里的F_N就是N点DFT矩阵F_M^H是M点IDFT矩阵。映射完再走经典的海森堡变换Heisenberg transform变成时域信号发射出去。接收端做反向操作维格纳变换Wigner transform把时域信号转回时频域再做SFFT就是ISFFT的逆变换最终得到收端的DD域信号。关键点在于因为信道在DD域是稀疏的经过整个收发链路后接收端的DD域信号和发射端DD域信号之间不是复杂的时变卷积而是一个由循环移位组成的稀疏卷积。每条路径在DD域引起的效果是发送符号在时延轴方向上平移了l_τ个格子在多普勒轴方向上平移了k_ν个格子再乘上一个复数增益。这个性质直接决定了后续可以用低复杂度的稀疏检测算法而不是去啃满矩阵。1.3 MIMO-OTFS引入的检测复杂度挑战单天线OTFS的接收信号能写成y H_eff · x w这里的H_eff是NM×NM的等效信道矩阵N是多普勒维度M是时延维度。因为信道稀疏所以H_eff每一行只有P个非零元素。P是路径数通常是个位数到十几的量级。但如果发射端有N_t根天线接收端有N_r根天线情况就变了。每个接收天线上收到的是所有发射天线发出的信号叠加y^{(r)} ∑_{t1}^{N_t} H_{eff}^{(r,t)} · x^{(t)} w^{(r)}把所有天线合在一起看系统的等效数据维度变成N_t·NM接收维度变成N_r·NM。当N_t N_r 4NM 8192时数据维度是32768接收维度也是32768。这个规模下想用最大似然ML检测复杂度是Q^(N_t·NM)级别的星座点穷举直接不现实。连MMSE检测都需要对32768×32768的矩阵求逆运算量同样顶着天花板。所以要跑MIMO-OTFS仿真低复杂度检测算法不是“可选项”而是“必选项”。这里的主流技术路线基本就是两条一条是MP算法这类基于稀疏矩阵的消息传递方案另一条是结合ZF/MMSE均衡和LU分解这类线性代数技巧的低复杂度方案。标题里提到的几个关键词正好对应这两条路线的核心环节。2. MIMO-OTFS链路模型从DD域符号到等效稀疏信道矩阵2.1 发射端完整处理链路先统一符号约定后面所有代码示例都以这个模型为准。假设一个OTFS帧有N个符号周期对应多普勒维度M个子载波对应时延维度总共NM个QAM符号。这NM个符号构成DD域发射矩阵X_dd维度是N×M。发射端处理分为三步首先是ISFFT。把DD域符号矩阵X_dd通过二维变换映射到时频域X_tfX_tf F_N · X_dd · F_M^H其中F_N是N点DFT归一化矩阵F_M^H是M点IDFT归一化矩阵。这里一定要确认归一化因子不同的文献定义差一个√(NM)的系数。我习惯用DFT矩阵的单元化形式即F_N的第(n,k)个元素是(1/√N)exp(-j2πnk/N)这样F_N·F_N^H I_N能量守恒SNR定义不会出错。ISFFT这一步本质是把DD域符号“打散”到时频域让每个时频格点上的符号都是所有DD域符号的加权叠加这是OTFS获得分集增益的数学基础。第二步是海森堡变换。把时频域符号矩阵X_tf通过脉冲整形滤波器变成时域信号。离散实现时常用的做法是用OFDM调制器对X_tf的每一行对应一个OFDM符号周期做M点IFFT加上循环前缀CP后按时间顺序发送。每个符号周期的采样点数是MGG是CP长度所以一帧时域采样点总数是N(MG)。第三步是MIMO扩展。如果有N_t根发射天线每一根本来可以独立发送不同的数据流但在OTFS帧结构设计中常用做法是所有天线共享同一个时频资源各自发送不同的DD域符号块。也就是说第t根天线上的发射矩阵X_dd^{(t)}放在DD域经过各自的ISFFT和海森堡变换后发射出去。天线之间没有时分/频分隔离空间维度完全打开。2.2 接收端处理与DD域等效信道接收端第r根天线收到所有发射天线信号的叠加。每条收发链路(r,t)之间有自己独立的时延-多普勒信道响应。接收处理流程是先做维格纳变换。对时域接收信号去CP按每个符号周期做M点FFT得到时频域接收矩阵Y_tf^{(r)}维度是N×M。再做SFFT。对Y_tf^{(r)}做二维变换回到DD域Y_dd^{(r)} F_N^H · Y_tf^{(r)} · F_M做完这两步如果忽略CP引起的边缘效应可以把Y_dd^{(r)}写成Y_dd^{(r)} ∑_{t1}^{N_t} H_{dd}^{(r,t)} * X_dd^{(t)} W_dd^{(r)}这里的*是二维循环卷积。H_{dd}^{(r,t)}是第t根发射天线到第r根接收天线之间的DD域信道响应矩阵维度是N×M。由于每条路径对应一个时延偏移和一个多普勒偏移H_{dd}^{(r,t)}只有在路径对应的(ν, τ)格点上有非零值。把二维卷积按列拉直column-wise vectorization就能得到前面提到的向量形式y H_eff · x wx的维度是N_t·NM × 1所有发射天线DD域符号的串接y的维度是N_r·NM × 1H_eff的维度是N_r·NM × N_t·NM。因为每条收发链路的DD域信道都稀疏H_eff整体也是稀疏的非零元素总数大约是N_r·N_t·P·NM而不是N_r·NM × N_t·NM个。这里的稀疏比大概是如果NM8192P8N_tN_r4非零元素总共约104万而矩阵总元素数是(32768)²≈10.7亿。稀疏度只有0.1%左右。MP算法能高效收敛靠的就是这个稀疏性。2.3 稀疏结构MP算法能成立的根基H_eff的稀疏结构具体长什么样值得多花点时间理解。拿单天线OTFS来说DD域卷积可以写成如下形式Y_dd[n, m] ∑_{i1}^{P} h_i · X_dd[(n - l_i) mod N, (m - k_i) mod M]这里h_i是第i条路径的复增益l_i是时延轴偏移单位格子数k_i是多普勒轴偏移单位格子数。注意一个容易踩坑的细节时延偏移和多普勒偏移在离散DD域里是循环移位不是线性移位。也就是说符号从时延轴的一端“绕到”了另一端。这个循环特性直接来自OFDM调制里的循环前缀设计。如果信道时延扩展超过了CP长度循环特性就会被破坏产生符号间干扰这也是OTFS仿真里CP长度选择不当会导致BER地板效应的原因之一。到了MIMO-OTFSH_eff的每个块H_{eff}^{(r,t)}都是N·M × N·M的稀疏矩阵每个块内部非零元素的位置由收发链路(r,t)各自的多径时延-多普勒参数决定。这些块拼成大矩阵后因为天线间信道传播条件常常不同非零元素的分布也是千奇百怪但整体稀疏性几乎不受影响。做仿真的时候不建议直接构造完整的H_eff大矩阵然后做运算——NM8192时每行还有8个非零元全矩阵存储要占几百MB内存太浪费了。正确做法是用稀疏矩阵结构MATLAB里就是sparse类型存储或者干脆连矩阵都不存只存信道参数在算法里按索引动态计算。后面讲MP算法实现时会详细说。3. MP检测算法逐行拆解3.1 为什么最终选MP而不是ML/MMSE很多刚开始做OTFS检测的人第一个想法就是ML或者MMSE结果一仿真就卡住了。ML检测的最优性没问题但复杂度是Q^(N_t·NM)Q是星座大小。哪怕N_t1NM1024Q4QPSK搜索空间也有4^1024这个数字连超级计算机都算不过来。所以ML在MIMO-OTFS里只能作为小规模对比的上界参考不是工程方案。MMSE检测复杂度看起来“只是”矩阵运算但它的核心是求解x_hat (H^H H σ² I)^(-1) H^H y这里涉及对N_t·NM × N_t·NM维矩阵求逆。前面算过NM8192、N_t4时矩阵维度是32768。直接求逆在这个规模下无论是存储还是运算量都是不现实的。就算用稀疏矩阵求逆技巧矩阵的填入量fill-in也会把稀疏结构破坏得面目全非。MP算法看中的恰恰是H_eff的稀疏性。它的复杂度是O(N·M·P·I)其中I是迭代次数。代入典型值8192×8×10≈65万次乘加运算和MMSE求逆动辄10^11 flops相比差了五个数量级。从工程角度看这才是能实时处理的量级。当然MP不是没有代价。它假设干扰项服从高斯分布在高信噪比和强多径环境下高斯近似会变差导致BER平台效应。所以实际系统通常是MP加误差纠正级联而不是单打独斗。3.2 因子图与消息更新均值、方差和阻尼MP算法的核心思想来自置信传播Belief Propagation。系统模型y H_eff·x w里每一个观测值y_i和若干个发射符号x_j相关联。以观测节点y_i为一方、变量节点x_j为另一方可以画一张因子图。因为H_eff稀疏每个y_i只连接P·N_t个x_j因子图非常稀疏消息传递收敛很快。第一次接触MP容易在公式上绕晕我建议按下面这个流程理解第一步初始化。每个符号x_j的先验均值为0方差为星座符号功率归一化后通常是1。然后进入迭代。第二步观测节点更新。对每个观测y_i它收到来自符号x_j的信息后可以把y_i写成y_i H_{i,j} x_j ∑_{k≠j} H_{i,k} x_k w_i把干扰项∑_{k≠j} H_{i,k} x_k w_i当作一个高斯随机变量它的均值和方差由上一轮迭代中其他符号的估计值决定。有了这两个参数就能算出给定x_j时y_i的条件概率也就是观测节点发给变量节点的消息。第三步变量节点更新。每个符号x_j连到多个观测节点对应它参与贡献的所有接收天线上的观测。x_j收到所有外部观测的消息后把这些消息和先验信息相乘归一化后得到x_j的后验分布。后验分布的期望和方差就是这一轮x_j的软估计。第四步也是最容易忽略的一步阻尼。直接把新估计的期望传给下一轮迭代算法在高SNR下很容易振荡一会儿收敛一会儿发散。我试过的做法是加一个阻尼系数β0.4~0.7更新规则变成μ_new β·μ_computed (1-β)·μ_old阻尼系数的作用相当于给迭代过程加了“惯性”防止消息在几个状态之间来回跳。β的具体取值要看信道条件推荐在0.5附近做几次预仿真再定。整个算法迭代I次后对每个x_j的后验分布做硬判决取星座上距离期望最近的点就是最终的检测输出。如果要送信道译码器就把后验分布转成LLR。3.3 一段带详细注释的MATLAB核心代码下面这段是我在MIMO-OTFS仿真里实际用过的MP检测核心代码做了一定简化以突出主要逻辑。注释我写得很细因为仿真代码这东西过三个月自己回来看都容易忘何况别人要复现你的结果。function x_hat mp_detector_mimo_otfs(y_vec, H_cells, N, M, Nt, Nr, ... constell, max_iter, damping) % y_vec : Nr*N*M 维接收向量DD域按列拉直 % H_cells : Nt*Nr 个 cell每个 cell 是稀疏矩阵 H_eff{r,t}维度 N*M x N*M % N, M : DD域网格维度多普勒、时延 % Nt, Nr : 发射、接收天线数 % constell : 星座符号向量 % max_iter : MP最大迭代次数 % damping : 阻尼系数 (0,1)建议 0.5 附近 % x_hat : 检测输出拉直后的 DD 域符号向量 K N * M; % 每根天线上的符号数 Nsym Nt * K; % 总发射符号数所有天线 Nobs Nr * K; % 总观测数所有接收天线 Q length(constell); % 星座点数 % 把稀疏信道的非零元素位置存下来避免每次迭代重复查找 row_idx cell(Nt, Nr); col_idx cell(Nt, Nr); val_idx cell(Nt, Nr); for t 1:Nt for r 1:Nr [I, J, V] find(H_cells{t, r}); row_idx{t,r} I; col_idx{t,r} J; val_idx{t,r} V; end end % 初始化变量节点符号的均值和方差 mu_x zeros(Nsym, 1); var_x ones(Nsym, 1); % 归一化星座功率默认为1 % 观测节点的噪声方差 sigma2_n 1 / (10^(snr_dB/10)); % 根据SNR设置实际代码中作为参数传入 for iter 1:max_iter % 第1步观测节点计算干扰统计量 % 对每个观测 y_i计算所有相连符号的干扰均值/方差 xi_hat zeros(Nobs, 1); % 干扰噪声的均值 xi_var sigma2_n * ones(Nobs, 1); % 干扰噪声的方差 for t 1:Nt for r 1:Nr H H_cells{t, r}; % 对每根发射天线符号位置偏移是 (t-1)*K sym_offset (t-1) * K; obs_offset (r-1) * K; [I, J, V] find(H); idx (I - obs_offset - 1) * Nsym (J - sym_offset); % 这行MATLAB代码只是为了演示索引计算逻辑 % 实际实现时建议直接用稀疏矩阵乘法完成均值更新避免循环 end end % 实际高效的实现是利用稀疏矩阵乘 % 1) 构建去掉目标符号后的符号均值向量 % 2) 对每个观测节点用 H_eff 乘以该向量得到干扰均值 % 3) 对干扰方差做类似计算 for t 1:Nt for r 1:Nr sym_offset (t-1) * K; obs_offset (r-1) * K; H H_cells{t, r}; % 干扰均值计算 xi_hat_block H * mu_x(sym_offset1 : sym_offsetK); xi_hat(obs_offset1 : obs_offsetK) ... xi_hat(obs_offset1 : obs_offsetK) xi_hat_block; % 干扰方差计算 xi_var_block (H.^2) * var_x(sym_offset1 : sym_offsetK); xi_var(obs_offset1 : obs_offsetK) ... xi_var(obs_offset1 : obs_offsetK) xi_var_block; end end % 第2步变量节点更新后验均值/方差 mu_new zeros(Nsym, 1); var_new zeros(Nsym, 1); for t 1:Nt sym_offset (t-1) * K; for r 1:Nr obs_offset (r-1) * K; H H_cells{t, r}; [I, J, V] find(H); % 对每个非零元素计算其对自己对应观测节点的贡献 for nz 1:length(V) i_obs I(nz) obs_offset; % 全局观测索引 j_sym J(nz) sym_offset; % 全局符号索引 h_val V(nz); % 剔除符号j_sym自身贡献后的观测残差 resid y_vec(i_obs) - (xi_hat(i_obs) - h_val * mu_x(j_sym)); % 对应的方差 resid_var xi_var(i_obs) - (h_val^2) * var_x(j_sym); resid_var max(resid_var, 1e-12); % 防止数值出负 % 计算该观测对符号j_sym贡献的消息 msg_var 1 / ( (abs(h_val)^2) / resid_var 1/var_x(j_sym) ); msg_mean msg_var * ( ... conj(h_val) * resid / resid_var ... mu_x(j_sym) / var_x(j_sym) ); % 累加到后验统计量 var_new(j_sym) var_new(j_sym) 1/msg_var; mu_new(j_sym) mu_new(j_sym) msg_mean / msg_var; end end % 转换回均值/方差形式 for j 1:K idx_sym sym_offset j; if var_new(idx_sym) 0 var_new(idx_sym) 1 / var_new(idx_sym); mu_new(idx_sym) mu_new(idx_sym) * var_new(idx_sym); else var_new(idx_sym) 1; mu_new(idx_sym) 0; end end end % 阻尼 mu_x damping * mu_new (1 - damping) * mu_x; var_x damping * var_new (1 - damping) * var_x; end % 硬判决找到星座中距后验均值最近的点 x_hat zeros(Nsym, 1); for j 1:Nsym [~, q] min(abs(constell - mu_x(j))); x_hat(j) constell(q); end end这段代码里有个关键点要特别提醒变量节点更新那一段如果写成内层循环遍历所有非零元素NM8192、N_tN_r4时非零元素约百万量级还会循环I次MATLAB跑起来会非常慢。我的习惯是先写清楚逻辑再逐步向量化。实际仿真时把内层循环改成批次处理或者用gather/scatter操作能快两个数量级。另外阻尼的初始可以从0.5开始调。我见过有人把阻尼设成0.9结果迭代20轮还在振荡BER曲线跟心电图似的。这个参数真不是越大越好。4. ZF均衡与低复杂度LU分解4.1 ZF的数学本质与噪声放大问题MP算法虽好但不是所有场景都适合用它做第一级检测。当信道矩阵的条件数很大、多普勒扩展又特别严重时MP的收敛速度会下降这时候很多实现选择先用线性均衡做一个粗估把符号的大致位置定下来再交给非线性算法细化。ZF迫零均衡就是这个粗估环节里最常见的方案。ZF的思想很直接不考虑噪声只求让干扰完全消失。对单天线系统均衡矩阵取信道矩阵的伪逆W_zf H^† (H^H H)^(-1) H^H作用在接收向量上x_hat W_zf · y x (H^H H)^(-1) H^H w好处是干扰项被彻底消除坏处也暴露在这里噪声项被H的特征值倒数放大了。当H_eff存在很小的奇异值时噪声放大效应尤其明显。在低SNR区间ZF的BER比MMSE差不少中高SNR区间两者的差距会缩小。MIMO-OTFS里用ZF还有另一个麻烦之前说了H_eff维度可以到32768×32768直接算伪逆中间要算H^H H这个Gram矩阵的逆。这个矩阵虽然不是满的但它的逆矩阵会变得很稠密存储都是问题。所以要在这里谈低复杂度核心思路不是去“求逆”而是把“求逆”改写成“解线性方程组”。4.2 LU/Cholesky分解怎么把“求逆”变成“求解”线性代数里有个基本的观念转换你几乎永远不需要显式求一个矩阵的逆。要算出x A^(-1) b真正需要做的是解线性方程组Ax b。解方程组的效率比求逆高得多。LU分解的思路是把矩阵A拆成下三角矩阵L和上三角矩阵U的乘积A L · U然后Ax b变成L·U·x b先解L·z b前代forward substitution再解U·x z回代backward substitution。三角矩阵的求解复杂度是O(n²)非常便宜。整个过程的工程量集中在LU分解本身不管待多少个b向量分解只需做一次。如果A是Hermitian正定矩阵还可以用Cholesky分解比通用LU再省一半运算量A L · L^HMMSE均衡里需要处理的矩阵是A H^H H σ²I它天然就是Hermitian正定矩阵所以Cholesky是首选。在MATLAB里用\运算符就会自动选择Cholesky路径比如x A \ b完成的就是分解加回代的组合流程。这么说可能还是抽象直接看规模差异。假设系统等效矩阵维度是n32768MIMO-OTFS的典型值直接求逆约需2n³次浮点运算大约7×10^13 flops。LU分解需要(2/3)n³次约2.3×10^13 flops。Cholesky又砍一半约1.2×10^13 flops。虽然LU在绝对数值上还是很大但相比直接求逆已经是数量级的进步而且分解可以复用信道在一个OTFS帧内不变分解一次后续所有符号块都只需要做两次回代。在MIMO-OTFS仿真里ZF均衡的具体实现流程是构造A H_eff^H · H_eff这一步用稀疏矩阵乘法H_eff本身稀疏A的稀疏度会下降但仍有大量零元对A做Cholesky分解得到L对每个b H_eff^H · y解L·z b和L^H·x_hat z输出x_hat作为ZF的软估计由于H_eff本身是稀疏的在实际实现中可以走稀疏LU/Cholesky路径MATLAB的lu函数支持稀疏矩阵会做重排序优化比如列近似最小度算法进一步把分解复杂度降到接近O(n^1.5)甚至更低。我仿真时对比过稀疏Cholesky在NM8192、P8的场景下分解时间从直接法的几分钟降到了十秒以内效果非常明显。4.3 复杂度对比从flops数量级看差距把前面几种方案的复杂度放在一起看更直观。以下按单天线OTFSNM8192P8MP迭代10次估算算法主要运算量相对数量级ML穷举QPSKQ^(NM)4^8192不可行MMSE直接求逆O(NM)³ ≈ 5.5×10^11 flops最高ZFCholeskyO(NM)³/3 ≈ 1.8×10^11 flops比求逆省2/3ZF稀疏Cholesky实际约10^9~10^10 flops利用稀疏性再降1~2个量级MP检测迭代10次O(NM·P·I) ≈ 6.5×10^5 flops最低中间差了大概五个数量级。这意味着在实时系统里MP能做到纳秒层级延迟而矩阵求逆需要专门的加速器件配合。这也是为什么OTFS的接收机研究路线里MP以及配套的低复杂度变体会成为主流。但做个总结性判断MP并不总是优于ZFLU。在信道稀疏度极低P很大比如大于50、或者符号帧尺寸很小的场景LU分解的固定开销占比上升两者差距缩小。反过来MIMO天线规模增大时H_eff维度以N_t·N_r倍数增长MP的优势会更明显。工程上最稳妥的是实现两条路子按信道条件切换。5. 误差纠正检测之后如何再做一次“清理”5.1 误差来源梳理噪声放大、信道估计误差与硬判决误差仿真跑多了会发现MIMO-OTFS检测器的BER曲线在高SNR区域往往降不下去像一个平台一样横在那里。平台高度通常就是误差纠正环节要处理的“底线”。先理清误差从哪里来。第一类是噪声放大。ZF均衡在信道奇异值小的方向上会把噪声放大几倍甚至几十倍。MP算法在高斯近似失效时也会产生类似效应。这类误差在低SNR下尤其突出。第二类是信道估计误差。OTFS的导频设计通常是在DD域放一个脉冲导频接收端通过检测导频的延迟-多普勒扩展来估计信道。但MIMO-OTFS里导频会被多个天线共享或复用互相产生干扰估计出来的H_eff和真实信道有偏差。信道估计误差会在整个检测过程中被放大。第三类是硬判决误差的传播。MP和ZF输出的软估计在判决后可能出错如果后续有迭代干扰消除错误的判决会参与下一轮干扰重建把误差继续传递下去。5.2 迭代干扰消除与MP的联合方案针对硬判决误差传播常规做法是迭代干扰消除Iterative Interference Cancellation, IIC。以前做MIMO检测时我用过串行干扰消除SIC和并行干扰消除PIC在OTFS里同理但顺序选择上有些讲究。SIC的做法是先对当前符号按可靠度排序从最可靠的那个开始硬判决后重建它的接收贡献从接收信号里减掉。减完再处理下一个次可靠的符号。这个方案每一步都在缩减后续符号受到的干扰收敛通常很快但有一个致命问题如果排在最前面的符号判决错误错误会被一路传播下去。所以排序依据不能只看瞬时SNR最好结合MP输出的后验方差来判断可靠度。PIC的做法是所有符号同时重建干扰减去后再统一更新判决。它不依赖判决顺序天然适合并行计算但每轮迭代每个符号的判决错误都会影响其他符号的更新高SNR下容易振荡。我实际用的方案是“MP后验方差加权”的中间路线MP收敛后每个符号都会有一个后验方差方差小的说明该符号的判决可信度高。第一轮只让方差低于某个阈值的符号参与干扰重建其余符号的干扰用软估计均值代替不做硬判决。这样既避免了SIC的误差传播问题又不至于像PIC那样振荡。阈值随迭代逐步放宽最后一轮全部符号做硬判决。信道译码的Turbo式迭代也属于误差纠正的大范畴。检测输出LLR后送给LDPC译码器译码得到的软信息反馈回检测器替代先验信息重新跑MP。这种“检测-译码-再检测”的外迭代一般做两三轮即可再增加轮数收益非常有限但复杂度成倍上涨。实测数据是从第1轮到第2轮外迭代BER能改善1~2个dB第3轮以后改善不超过0.2dB。实际权衡后仿真和产品实现基本都停在两轮。5.3 收敛判断与异常BER曲线的排查MP迭代什么时候停不能光靠固定轮数。我在代码里加了两个指标辅助判断。第一个是符号翻转率。每轮迭代后统计硬判决结果相对上一轮变化的符号比例。如果翻转率低于某个阈值比如0.1%说明算法已经收敛继续迭代意义不大。第二个是残差能量。把当前估计的符号重建接收信号算接收残差的能量相对初始残差能量的下降比例。这个指标更接近“信号域”的收敛判断比符号域的翻转率更稳因为符号可能处于星座判决边界上来回变化但残差其实已经很小了。如果BER曲线出现异常比理论线高很多或者存在非预期地板我开始排查时的顺序是第一步看信道估计是否正确检查导频位置和CP长度第二步看MP是否真的收敛了打印每轮迭代的残差能量看是不是振荡第三步看SNR定义是否一致发射信号功率归一化有没有对。这三个坑覆盖了我在OTFS仿真里遇到的绝大多数异常现象。曾经有一次我BER高得离谱查了两天最后发现是SNR定义里忘了把MIMO的天线数算进去发射功率按单天线归一化了接收端多了3dB噪声。这种错误特别隐蔽因为链路每一级看起来都是对的缝合起来就错了。6. 仿真踩坑与参数设置实践6.1 波形参数N和M怎么定很多入门的仿真正是从“随便填一个N和M”开始的结果仿真结果要么边界效应严重要么多普勒分辨率不够反映在BER曲线上就是怎么调都很难看。这里分享一套我调试总结出来的参数选择思路不一定最优但至少不会让仿真结果被参数“坑”掉。N对应多普勒域的格子数。DD域的多普勒分辨率是Δν 1/(N·T)T是一个符号周期。要看清最大多普勒频移f_dmax至少需要N·T ≥ 1/f_dmax也就是N·Δf要多普勒分辨率覆盖到最大频移。典型的做法是让最大多普勒频移对应的格子数k_max不超过N的几个百分比最多30%。如果N太小多普勒分辨率粗糙两条相隔很近的多普勒路径在DD域里会混叠成一个点信道矩阵的稀疏性被破坏MP的收敛性立刻恶化。M对应时延域的格子数。时延分辨率和子载波间隔挂钩Δτ 1/(M·Δf)。需要覆盖最大时延扩展τ_max即M ≥ τ_max·Δf 1。同时M通常也决定了FFT点数一般选2的幂次方便实现。我常用的一个组合是N128M64子载波间隔15kHz符号周期66.7μs。对应的多普勒分辨率大约117Hz能覆盖约15kHz的最大多普勒范围对应极端高铁场景绰绰有余时延分辨率约1.04μs能分辨约66μs的最大时延扩展。覆盖5G城市信道模型的典型参数完全没问题。6.2 稀疏信道矩阵构造的三种易错点第一个易错点就是线性卷积和循环卷积搞混。前面说过OTFS的DD域系统模型是循环卷积这个循环性来自CP。但很多人在构造H_eff的时候下意识地用了线性卷积的索引方式导致符号块边缘出现虚假干扰BER性能变差。正确的构造方法是对每条路径的时延偏移l和归一化多普勒偏移k信道矩阵H_eff的非零位置用模运算计算行索引(m k) mod M列索引(n - l) mod N。我试过把这个模运算漏掉的版本误码率高到像没做均衡一样。第二个易错点是路径增益的复系数归一化。多径信道的总功率应该归一化到1每条径的功率按路径损耗模型分配。如果忘了归一化接收SNR和理论设定的SNR就对不上MP输出的后验方差也会失真。更隐蔽的是如果各径功率相差悬殊弱径那条路的干扰重建精度很低误差打消效果也差。第三个易错点是天线间信道的相关性建模。做MIMO仿真时如果不同收发天线对之间的信道完全独立生成那么空间分集增益是最理想值BER好得过分乐观。实际天线间距有限不同天线对之间的信道存在相关性。严谨一点的做法是引入天线相关系数矩阵按指数相关模型对独立信道做Cholesky变换后再使用。学术仿真可以按独立信道做但结论要注明这是理想空间相关下的上界。6.3 让仿真可复现的几个好习惯仿真代码这东西除了自己看懂还要经得起别人复现。我在整理MIMO-OTFS仿真代码时坚持以下几个习惯。固定随机种子。信道、噪声、数据比特每一层的随机生成器都要有明确的种子控制。不要用全局默认的rng因为MATLAB版本升级默认种子会变化别人和你的结果可能有微小差异。我的做法是信道生成用一个独立的随机流对象RandStream数据生成用另一个这样修改噪声种子不会影响信道和数据。SNR定义写清楚。我习惯定义归一化符号能量Es1噪声方差按SNR 10*log10(1/σ²)设置。多天线场景下要特别注意是在“单符号”粒度还是“单根天线”粒度上定义SNR两套之间差10log10(N_t)dB。在代码注释里把这个关系写明白可以省掉大量复现时的争议。记录版本信息和参数哈希。仿真参数表、算法版本号、代码commit号最好输出到结果文件里。我做对比实验时经常面对十几组BER曲线没有参数记录简直没法溯源。最简单的做法是在主程序开头把关键参数打包成一个struct存到结果mat文件里。性能计时用tic/toc还不够要精细到每个功能模块。我习惯在信道估计、MP检测、LU分解、迭代纠错四个模块各放计时点输出到日志。这样如果代码改动导致性能退化能第一时间定位是哪个环节变慢了。6.4 一份推荐的仿真工作流把整个MIMO-OTFS链路仿真串起来我推荐按下面这个顺序推进第一步先跑通单天线OTFS的MP检测确认BER曲线在低SNR区间符合预期。这一步主要是验证信道模型和发射接收链路没有方向性错误。第二步加入ZF和LU分解实现对比MP和其他均衡器的性能差距。这里的重点是验证复杂度分析的实际效果用tic/toc记录各算法运行时间。第三步扩展为MIMO-OTFS先N_tN_r2天线间用独立信道验证MP算法的天线可扩展性。然后再上4×4。第四步加入信道译码和误差纠错外迭代观察BER曲线是否出现平台以及平台高度是否低于目标值。如果平台过高回溯检查误差来源。第五步整理对比表格输出各个模块的运行时间、BER at target、复杂度统计。这套流程每一层都有明确的验证基准出了问题能很快定位是算法问题还是实现问题比直接堆一个全功能大链路再回头debug要省心得多。最后再分享一个小体会OTFS仿真里“复杂度分析”和“运行时间”是两回事。MP算法的flops数量级虽然低但如果实现里到处是循环和重复索引计算实际运行时间可能比LU分解还长。仿真时不要看理论复杂度就省略向量化优化那会让结果失真。我现在写MP检测都会保持“先正确、再优化”的顺序第一版保证逻辑对第二版再做向量化。这样即使优化过程出了bug也知道该查哪里。