简介：面向高移动性通信场景的OTFS（正交时频空间调制）研究，这份代码提供低复杂度迭代LMMSE-PIC均衡器的实现。OTFS中每个符号经历相同的双选择信道衰落，配合适当均衡器可充分利用全时频分集；该实现以线性最小均方误差并行干扰消除为核心，并采用一阶诺伊曼级数近似矩阵求逆，将复杂度降至与传输符号总数呈准线性关系。代码面向通信工程与信号处理领域的研究者、研究生，尤其适合正在探索OTFS信号检测或期望在编码OTFS系统中复现均衡性能的读者。压缩包大小约15.84MB，包含可直接运行的仿真程序与相关文件。作者提供了编码误码率与EXIT图仿真所需的实现，便于对比文中提出的LMMSE-PIC均衡器与既有PIC方法、OFDM方案的性能差异。已有798人学习下载，适合作为OTFS迭代接收机设计与低复杂度算法研究的参考实现。
1. 高速移动场景下的信号检测：LMMSE-PIC均衡器为什么能拿来做OTFS解调
正交时频空间调制（OTFS）在高铁、低轨卫星这类高多普勒场景下，把调制搬到了延迟-多普勒域。信道在这个域里变化缓慢，每个发射符号经历的衰落基本一致，这是OFDM做不到的。但这个优势不是白给的——接收端要做二维均衡，复杂度极高，很多论文只能做到仿真验证。这份代码实现了LMMSE-PIC（线性最小均方误差-并行干扰消除）迭代均衡器，用「先线性滤波再逐步消去干扰」的思路，把复杂度压到可接受范围。适合正在做OTFS信道估计与信号检测的通信算法工程师，也适合想验证LMMSE-PIC迭代思想的研究者——但我得先说实话：这份代码不简单，需要你先理解信号模型，再去跑仿真，不然会看得一头雾水。
2. OTFS信号模型与均衡器选型：从延迟-多普勒域说起
2.1 OTFS收发机的基本结构
OTFS的发射端不是直接把数据符号映射到时频域，而是先通过逆辛有限傅里叶变换（ISFFT）把延迟-多普勒域的符号映射到时频域，再经过海森堡变换得到时域发送信号。接收端做对称操作，先维格纳变换回时频域，再通过辛有限傅里叶变换（SFFT）回到延迟-多普勒域。
代码里的收发链路大致长这样的核心流程：
% 发送端：延迟-多普勒域符号 -> 时频域 -> 时域 dd_grid = symbol_placement(data_bits, M, N); % M个子载波，N个OTFS符号 tf_grid = isfft(dd_grid, M, N); % 逆辛有限傅里叶变换 tx_signal = heisenberg_transform(tf_grid, M, N, cp_len); % 海森堡变换加循环前缀关键点在于调制参数：M是子载波数，N是时隙数，整个OTFS帧承载M乘以N个调制符号。延迟-多普勒域的网格分辨率由M和N的乘积决定，延迟分辨率是1/(Mdelta_f)，多普勒分辨率是1/(NT)。代码中所有矩阵维度都是M乘N，搞混维度是后面调试时最常见的翻车原因，建议先打印size验证。
2.2 为什么选LMMSE-PIC而不选其他均衡器
OTFS的均衡方案大致分三类：线性均衡器（LMMSE/ZF）、消息传递类检测器（MP/EP）、以及迭代干扰消除类（PIC/SIC）。线性均衡最简单但性能差，尤其在分数多普勒信道下误差平台很高。消息传递检测性能好，但实现复杂，每轮迭代要维护高斯消息的均值和方差，代码调试难度大。这份代码选LMMSE-PIC是权衡了性能和工程可实现性的方案——用一次LMMSE估计作为初始解，然后通过并行干扰消除逐轮逼近真实符号。相比逐符号串行消除（SIC），PIC的并行结构天然适合矩阵运算，运行速度会快不少。
复杂度方面，直接对M乘N维的有效信道矩阵做MMSE估计，矩阵求逆开销是O((MN)^3)，完全不可用。这份代码的核心技巧是把均衡拆成两步：初始段用频域LMMSE做粗略估计，后续迭代只对每个符号逐个做低维滤波，复杂度降到O(MN * L)，L是信道抽头数，通常远小于MN。实际跑仿真时，M=64、N=32的场景下，10轮迭代在普通桌面机上几百毫秒内完成，这个速度对链路级仿真完全够用。
2.3 有效信道矩阵的构造与稀疏性利用
OTFS的输入输出关系可以写成 y = H_eff * x + w，其中H_eff是时延-多普勒域的有效信道矩阵。这个矩阵长什么样直接决定均衡的复杂度。如果信道只有整数时延和整数多普勒，H_eff内部每个符号只对应少数几个非零增益元素，矩阵是稀疏的。代码利用这个稀疏性，把逐符号PIC的干扰计算限制在非零抽头范围内，避免全矩阵乘法。
% 构造有效信道矩阵（稀疏形式） H_eff = sparse(M*N, M*N); for l = 1:L for k = 1:K idx = ...; % 当前抽头对应的收/发符号索引对 H_eff(idx, :) = H_eff(idx, :) + h(l,k) * circshift(...); end end索引对的计算是这段代码里最容易蒙圈的地方。时延抽头l会沿延迟维循环移位，多普勒抽头k会沿多普勒维循环移位，移动方向取决于信道冲激响应建模时用的傅里叶变换符号约定。不同论文对SFFT正反变换的归一化因子约定不完全一致，直接套用会导致误码率发散。我习惯先把l=1、k=1的简单信道输进去，对比均衡输出和发射符号是否对齐，再做完整信道验证。
3. 工程结构拆解：一份完整仿真链路的基本构成
3.1 目录结构与模块划分
这份代码按「发端-信道-收端-评估」四段式组织，模块边界干净。拿到手先不急着跑，花几分钟把目录结构过一遍，心里有数：
|-- main.m % 主仿真脚本，控制整个流程 |-- mod/ % 调制相关 | |-- isfft.m % 逆辛傅里叶变换 | |-- sfft.m % 辛有限傅里叶变换 | |-- heisenberg_transform.m % 海森堡变换 | |-- wigner_transform.m % 维格纳变换 |-- channel/ | |-- generate_channel.m % 生成延迟-多普勒域信道 | |-- apply_channel.m % 信道卷积 |-- equalizer/ | |-- lmmse_pic_equalizer.m % 核心LMMSE-PIC均衡器 | |-- initial_mmse_est.m % 初始MMSE估计 |-- utils/ | |-- modulate_bits.m % 比特到符号映射 | |-- demodulate_bits.m % 符号到LLR | |-- ber_calculation.m % 误码率统计主脚本main.m控制参数配置、循环信噪比点、调用各模块并汇总误码率曲线。模块之间通过结构体传递参数，比如用struct('M',M,'N',N,'delta_f',delta_f,'T',T)这样的配置对象贯穿全链路。
3.2 信道模型：整数多普勒与分数多普勒的生成
OTFS信道建模和OFDM最本质的区别在于信道抽头的位置不是单个时域延时加频偏，而是一组离散的整数时延和多普勒频移。代码里支持两种信道模式：整数多普勒模式下，多普勒频移恰好是多普勒分辨率1/(N*T)的整数倍，H_eff是严格稀疏的——这种情况下OTFS的时域-多普勒域双正交性完整保留，均衡结果理想；分数多普勒模式下，多普勒频移不在网格上，会在多普勒维产生能量泄漏，H_eff不再保持完美稀疏，但实际物理信道几乎总是分数多普勒（因为真实移动设备的径向速度很难精确对应到某个整数倍频偏），所以仿真一定要跑分数模式才贴近真实。
% 信道抽头生成示例：每个抽头包含时延、多普勒和复增益 channel.tau = [0, 1, 2]; % 归一化时延抽头位置（整数） channel.fd = [0, 0.5, 1.2]; % 归一化多普勒频移（分数值） channel.h = [1.0, 0.5, 0.3]; % 复路径增益信道抽头数、多普勒值的分布、路径增益，全部在这个结构体里配置。体验时如果用默认参数跑出来误码率高得离谱，先检查信道增益有没有归一化——路径增益的平方和不为1时，接收信噪比和理论值对不上，所有误码率曲线都会整体右移。
3.3 主仿真循环与参数配置
main.m的核心不是均衡器本身，而是循环信噪比点、重复多次取平均的蒙特卡洛框架。每次循环内部大体是这样：
for snr_idx = 1:length(snr_list) for trial = 1:num_trials % 生成随机的延迟-多普勒域发射符号 tx_dd = randsrc(M, N, constellation); % 发送端处理 [tx_signal, params] = otfs_modulate(tx_dd, params); % 过信道加噪声 rx_signal = add_noise(apply_channel(tx_signal, channel, params), snr_list(snr_idx), params); % 接收端处理 rx_dd = otfs_demodulate(rx_signal, params); % LMMSE-PIC均衡 est_dd = lmmse_pic_equalizer(rx_dd, channel, params); % 统计误码 bits_tx = symbol_to_bits(tx_dd, params.QAM_order); bits_rx = symbol_to_bits(est_dd, params.QAM_order); ber_acc = ber_acc + sum(bits_tx ~= bits_rx) / numel(bits_tx); end ber(snr_idx) = ber_acc / num_trials; end参数配置基本都集中在main.m开头，需要调整的几个关键量包括子载波数M、符号数N、调制阶数、循环前缀长度、信道多普勒最大频移、迭代次数和信噪比范围。这些参数之间的关系是：M乘N决定帧大小，N又通过1/(N*T)决定多普勒分辨率，多普勒分辨率又直接影响分数多普勒的泄漏程度——改的时候最好成组调整。
4. 核心实现：LMMSE-PIC均衡器的逐行拆解
4.1 初始估计：用频域LMMSE给后续迭代提供起点
LMMSE-PIC的第一步是获取一个初步的符号估计。这个初始值质量好坏直接影响最终迭代收敛性。如果初始估计完全不靠谱，PIC后面每一轮都在消除错误的干扰项，越迭代越差，最终误码率平台比不做迭代还高。代码里用的初始方案是在时频域做LMMSE滤波，利用信道在时频域的准对角结构做低复杂度估计，然后再变换回延迟-多普勒域。
function [x_init, V_init] = initial_mmse_est(y_tf, H_tf, sigma2, params) % y_tf: 频域接收信号（M x N） % H_tf: 时频域信道频率响应 % 对每个OFDM符号独立做MMSE均衡 x_tf = zeros(M, N); for n = 1:N % 频域LMMSE滤波，对角化利用了OFDM子载波正交性 x_tf(:, n) = (conj(H_tf(:, n)) ./ (abs(H_tf(:, n)).^2 + sigma2)) .* y_tf(:, n); end % 变换回延迟-多普勒域作为初始估计 x_init = sfft(x_tf, M, N); % 估计误差方差，用于后续PIC的软信息计算 V_init = mean(abs(H_tf).^2 ./ (abs(H_tf).^2 + sigma2), 'all'); end这段代码的巧妙之处在于把二维MMSE问题拆成N个独立的频域MMSE问题。每个子载波上的信道频率响应是标量，MMSE滤波退化为逐点乘法，复杂度只有O(MN)。代价是频域MMSE没有利用延迟-多普勒域信道矩阵的稀疏结构，从这个意义上做的是次优初始估计。但作为迭代起点完全够用——LMMSE-PIC的特点就是前两轮重建精度差一些，后面迭代能追回来。改成完全LMMSE的初始估计收益不大且复杂度高很多，得不偿失。这个V_init在PIC里的地位也很重要，它表示初始估计的置信程度，值越小说明初始估计越可信，PIC的干扰消除步长走得更短。
4.2 并行干扰消除：每一轮迭代做了什么
LMMSE-PIC每轮迭代的核心动作可以用三句话概括：利用上一轮估计的符号重建所有符号的接收贡献；对目标符号减去其他符号的干扰；对剩余信号做MMSE滤波修正。所谓并行指的是所有符号的干扰消除同时进行，不依赖符号间的先后顺序，天然适合向量化实现。
function [x_new, x_soft] = pic_iteration(y_dd, H_eff, x_prev, V_prev, sigma2, params) Q = params.M * params.N; % 1. 重建信号：基于上一轮估计的软符号 x_soft = tanh(x_prev / V_prev); % 软符号映射（BPSK示例） s_recon = H_eff * x_soft; % 2. 并行干扰消除：每个符号逐个处理 x_new = zeros(Q, 1); for q = 1:Q % 目标符号的信道向量 h_q = H_eff(:, q); % 消除来自其他所有符号的干扰 r_q = y_dd - s_recon + h_q * x_soft(q); % 干扰方差：由V_prev和信道增益共同决定 var_interf = sum(abs(H_eff(:, setdiff(1:Q, q))).^2, 2) * V_prev; % 3. MMSE滤波修正 w_q = conj(h_q) ./ (abs(h_q).^2 + var_interf + sigma2); x_new(q) = w_q' * r_q; end end这里面两个细节值得反复琢磨。第一个是软符号映射。直接把硬判决结果拿来做干扰消除当然更简单，但硬判决错误会直接传播，软符号（比如x_prev除以V_prev后做tanh）让低置信度的符号保持很小的幅度，对干扰的影响也小，避免错误雪崩。当然软符号只在QAM星座的对称结构下好算，高阶QAM要按最小欧氏距离做软映射而不是简单tanh。第二个是干扰方差的计算，代码里用st维的近似——所有其他符号的干扰功率按信道增益和V_prev的乘积估计。严格做法要对每个q计算完整的协方差矩阵，复杂度直接翻倍。这种近似在信道抽头数较少时误差不大，稀疏信道下性能损耗很小。
还有一个编码细节值得学习：r_q和x_new的循环是Q次，Q等于M乘以N。Q较大时（如256）纯串行循环在MATLAB里跑得很慢，可以改成矩阵化实现——把H_eff拆成二维数组索引，一次算出所有r_q。但矩阵化版本内存占用量很大，需要权衡仿真设备的内存来选。我一般先用循环版调通逻辑，再对性能瓶颈做矩阵化优化，不建议一上来就写高维矩阵运算——调试难度大太多。
4.3 收敛条件与迭代终止策略
代码里的迭代终止条件包括固定迭代次数与残余变化阈值两种。固定次数实现简单，LMMSE-PIC在稀疏信道下迭代5到8轮后误码率基本收敛，设置更多轮数收益很小。残余变化阈值则是比较相邻两轮估计值的归一化差异，小于阈值就提前终止，在低信噪比下能省不少时间。实际仿真中固定迭代次数更常用，因为误码率曲线的对比需要每轮条件完全一致。
function [x_final, iter_used] = lmmse_pic_equalizer(y_dd, H_eff, params) % 初始估计 [x_cur, V_cur] = initial_mmse_est(...); for iter = 1:params.max_iterations [x_next, x_soft] = pic_iteration(y_dd, H_eff, x_cur, V_cur, ...); % 收敛判定：相对变化小于阈值 if norm(x_next - x_cur) / norm(x_cur) < params.tol x_final = x_next; iter_used = iter; break; end x_cur = x_next; V_cur = update_variance(x_cur, x_soft, params); % 方差随迭代更新 end x_final = x_cur; endV_cur在迭代过程中也会更新，这是容易忽略但影响收敛性能的一环。随着迭代轮数增加，软符号逐渐靠近真实符号，估计误差方差应该逐步减小——如果方差保持不变，后续迭代的干扰消除步长永远一样，收敛速度会明显变慢。代码里的update_variance具体逻辑是统计软符号x_soft与当前估计x_cur的差异均值，折算成下一轮的等效噪声方差。方差递减节奏和信道条件有关，我们曾经遇到过方差下降太快导致迭代过冲的现象——后续轮次符号被硬拉向错误星座点，误码率反而上升。看到这种情况时，可以先给方差更新加一个松弛因子，让它平滑地降。
4.4 高阶QAM的支持方式
BPSK的软符号可以用tanh直接算，但OTFS系统实际里常用16QAM和64QAM乃至更高阶调制。高阶QAM的软符号计算要对每个星座点做概率加权平均，核心是算出软符号的期望值。代码里会看到根据调制阶数分支处理的部分，16QAM的软符号计算需要遍历16个星座点。
function x_soft = soft_symbol(x_prev, V_prev, constellation) % 软符号等于后验概率加权下的星座点期望值 num_points = length(constellation); % 计算每个星座点的概率权重 dist = abs(x_prev - constellation).^2; weights = exp(-dist / (2 * V_prev)); weights = weights / sum(weights); % 加权平均得到软符号 x_soft = sum(constellation .* weights); % 同时对星座点做能量归一化 x_soft = x_soft / sqrt(mean(abs(constellation).^2)); end这个实现方法在Q等于1024甚至2048时会有性能问题——每个符号要遍历星座点并算exp，但OTFS帧内符号数大幅增加时这块会成为瓶颈。工程上更快的做法是对QAM的高位比特和低位比特分别计算软信息，利用QAM的格雷编码特性化简计算，复杂度从指数级降到线性级。代码里用的是直接加权平均，正确性没有问题，但跑大规模仿真时要做好心理准备，可能得加个进度条。后续二次开发时想提速，可以先优化这个函数。
5. 避坑指南：OTFS仿真里常见的几个坑
5.1 延迟-多普勒网格尺寸与信道抽头不匹配导致频谱泄漏
现象：误码率曲线随信噪比增高后不再下降，出现平台，而且平台高度和信道抽头数相关。
原因：OTFS帧的延迟分辨率是1/(Mdelta_f)，多普勒分辨率是1/(NT)。如果信道抽头位置不落在网格点上，能量会泄漏到相邻格点。泄漏分量在实际均衡时表现为无法消除的残余干扰。比如分数多普勒fd=1.3抽头，在N=32的帧里落在1.3倍分辨率的非整数格点上，能量扩散到第1和第2个多普勒格点，LMMSE-PIC只能消除这两个格点的干扰，丢失的部分成为底噪。
解决：仿真时要明确区分理想整数信道和实际分数信道。整数信道用于算法验证，分数信道用于性能评估。分数信道的多普勒值用确定性值并保持和帧参数匹配——让fd乘以N*T后的小数部分尽量靠近0.1以内，减小泄漏。要做完整性能测试时，在多个随机信道样本上取平均，避免某组信道参数异常。
5.2 高阶调制下软符号方差失配导致的不收敛
现象：64QAM调制时，迭代均衡的误码率在第3轮后不再下降，甚至回弹。BPSK的时候没有这个问题。
原因：LMMSE-PIC的干扰消除精度建立在软符号估计方差V_prev的准确性上。64QAM星座点间距小，相同噪声下符号判决的置信度远低于BPSK，如果V_prev更新幅度过快，把不准确的软符号当成准硬判决来消除干扰，误差会累积放大。
解决：压制V_prev的更新速度。常见做法是引入平滑因子β，方差更新改为V_new = (1-β)V_old + βV_estimated，β取0.3到0.5。另一个思路是增加一个防回弹机制——如果本轮误码率比上一轮差，跳过本轮软符号更新，保留上一轮软符号进入下一轮。这两个改法都是几行代码的事，但效果很直接。
5.3 矩阵求逆数值不稳定导致估计结果异常
现象：初始LMMSE估计偶尔出现NaN或者某个符号估计值比正常星座点大几个数量级。低信噪比时尤其频繁。
原因：频域LMMSE滤波的分母是abs(H_tf).^2 + sigma2。当信道存在深衰落，某个子载波的abs(H_tf)接近0时，分母几乎等于sigma2，这本身不会数值爆炸，但如果sigma2被错误配置——比如以dB形式传入却忘了转换、或者归一化噪声功率除以了符号能量两次——分母会出现0或负值。
解决：把噪声方差的传递链检查一遍，从SNR换算到sigma2这一步用统一函数封装，别在多个脚本里重复计算。其次在LMMSE滤波中加一个下限保护：分母取max(abs(H_tf).^2 + sigma2, 1e-12)，保证任何情况下分母为正。
5.4 蒙特卡洛仿真样本不足导致误码率曲线抖动大
现象：低误码率区间（比如1e-4以下），误码率曲线上下抖动严重，看不出真实的性能对比。
原因：误码率统计是独立事件计数，误码率降到1e-4时，每次蒙特卡洛试验平均只产生0.1个错误比特——很多轮试验零错误，偶尔一轮错一两个比特，误码率跳动几个数量级是正常现象。如果试验次数只有几十次，这个统计方差完全没法接受。
解决：按目标误码率反推试验次数。至少要积累100个错误比特才能把误码率估计的置信区间收窄到可接受范围。所以仿真框架里要加一个判断条件：累积错误比特数低于100时自动增加试验次数，而不是固定次数跑完。这个改动对仿真时间的影响很大，但数据可信度提升明显。
6. 二次开发：把LMMSE-PIC迁移到自己的接收链路
要做二次开发，第一步是看清这份代码的边界和依赖，别把它当成一个黑匣子直接塞进现成链路。它的输出是均衡后的软符号估计，不是比特级软信息。如果你的系统接的是硬判决和信道译码器，可以直接拿星座硬判决后用，但性能会差一些——因为均衡器输出的软信息包含了置信度信息，直接丢弃太可惜。常见做法是把均衡输出的软符号和等效噪声方差折算成比特级对数似然比（LLR），再送到LDPC或Turbo译码器。
% 从软符号估计折算LLR的示意代码 function llr = symbol_to_llr(x_est, V_est, constellation, bit_map) % 对每个符号位置，遍历星座点计算LLR % bit_map: 星座点到比特的映射表（格雷映射） num_bits = log2(length(constellation)); for b = 1:num_bits % 第b位为0的星座点集合 set0 = bit_map(:, b) == 0; set1 = bit_map(:, b) == 1; % 计算到集合的欧氏距离 dist0 = min(abs(x_est - constellation(set0)).^2); dist1 = min(abs(x_est - constellation(set1)).^2); % 近似LLR = (最小距离差) / (噪声方差 + 残余干扰方差) llr(b) = (dist0 - dist1) / V_est; end end这个折算方式用的是max-log近似，实际性能接近理想LLR，计算量小很多，适合作为译码器输入。把V_est换成均衡器输出的最终方差即可，但记得V_est要包含干扰消除后的残余干扰，不光是信道噪声——代码里最终方差的计算逻辑里有这一步，改造时千万别省略。
如果要把均衡器接入信道译码器的迭代系统（IDD），LMMSE-PIC的每一轮输出还需要耦合译码器回传的外信息。这个改动的复杂度高一些——PIC内部迭代和译码器迭代要协调好节奏，我建议先在单轮均衡加译码的结构上验证，再逐步增加迭代轮数。常见做法是：均衡器跑1-2轮PIC，把LLR送译码器更新一次，译码器输出的外信息再作为下一轮PIC的先验信息。这种IDD结构能让OTFS在高速移动信道下逼近信道容量上界，是这个领域一个很值得探索的方向。
改造验证的方法也比较固定：先在仿真层面对比改动前后的误码率曲线，确认性能不劣化；再在固定信道下做收敛性分析，观察每一轮迭代的均方误差下降曲线。我个人的习惯是把前面的V_prev平滑因子、软符号映射函数这些参数暴露成可配置项，在仿真里快速切换对比——不要写死在函数内部，否则每次调参都要改代码。这份代码本身不是银弹，但作为基线均衡器去做OTFS接收链路扩展是够用的。想真正跑通完整OTFS链路，很多东西要自己积累经验——就像我在最初搭建仿真链路时，完全没有意识到分数多普勒泄漏会直接把误码率平台钉死在1e-2，这个问题困扰了我整整一周。从那以后，我每次跑OTFS仿真都强制自己先把信道参数和网格尺寸的匹配关系验算一遍，再动手写均衡代码。把这步前置后，后面调试积累的问题就少多了。希望这些经验帮你在自己的调试路上少踩几个坑。
本文还有配套的精品资源，点击获取
