简介OTFS调制是面向高速移动场景的先进无线通信技术这份MATLAB代码包聚焦其仿真实现与峰均功率比PAPR分析适用于通信方向研究者、工程师及高年级学生进行算法验证与性能评估。压缩包共7个文件均为.m脚本大小约9KB核心代码覆盖OTFS调制、时频网格构建、多径信道生成与输出、MP检测及解调等完整环节另有gaussianFilter.m用于信号预处理帮助降低噪声、提升时频映射准确度。已有2567人学习下载。借助这些代码可以厘清OTFS从发射到接收的仿真链路直观观察PAPR计算过程与控制效果理解不同抑制方法对放大器效率和误码性能的影响所有脚本均可直接运行且支持按需调整调制阶数、信道参数等变量适合作为二次开发的起点支撑课程设计、科研实验或论文复现。在高速移动场景下PAPR的合理控制有助于提升系统能效本包可直接用于相关课题的快速验证与分析。1. OTFS 与 PAPR为什么高速移动场景下这个 MATLAB 代码包值得拆开看正交时频空间调制OTFS在近年来的物理层研究中频繁出现核心卖点是把调制符号直接铺在时延-多普勒Delay-Doppler, DD域上而不是传统的时频域。这样做的好处很直接高速移动场景下信道随时间快速变化OFDM 的子载波正交性容易被多普勒扩展破坏而 OTFS 在 DD 域里把每个符号经历的信道变成近似时不变的稀疏响应接收端用消息传递MP类算法就能把符号完整捞回来。代价是波形本身的峰均功率比PAPR偏高放大器回退深度加大功耗和线性度都受影响。这个压缩包里给的正是 OTFS 发射链路到接收检测的完整 MATLAB 实现并明确带了一组 PAPR 分析和处理的代码适合两类人一类是做物理层算法仿真的研究生和工程师想搭一条完整的 OTFS 基带链路作为对照基线另一类是关注波形峰均比问题、想验证 PTS 等抑制手段在 DD 域波形上是否有效的人。后面所有代码分析都基于包内实际文件组织展开重点落在怎么跑通、怎么读懂中间变量、以及 PAPR 曲线怎么画才严谨。2. 从 DD 域映射到波形OTFS 调制与解调的 MATLAB 实现拆解2.1 时频网格与 DD 域符号放置的逻辑OTFS 的发射端不是直接把 QAM 符号扔到子载波上而是先构造一个 M×N 的 DD 域网格。M 对应时延维度的分辨率通常等于子载波数N 对应多普勒维度的分辨率通常等于一个 OTFS 帧内的 OFDM 符号数。物理参数上子载波间隔 Δf 决定多普勒覆盖范围±Δf/2符号周期 T 决定时延覆盖范围最大可分辨时延为 1/Δf。选择 M64、N16、Δf15kHz 时一个帧的时间跨度是 16×66.7μs≈1.07ms多普勒分辨率约为 1/(N·T)≈937.5Hz最大可测多普勒约为 ±7.5kHz。这个配置在 3.5GHz 载频下对应约 ±770km/h 的移动速度正好覆盖高铁场景。代码里OTFS_modulation.m做的事情就是把随机生成的 QAM 符号矩阵 X_ddM×N通过逆辛有限傅里叶变换ISFFT映射到时频域再用海森堡变换生成时域发送波形。具体可以拆成三步理解。2.2 ISFFT 与海森堡变换的代码实现第一步是 ISFFT对 DD 域符号矩阵做 N 点 DFT 和 M 点 IDFT注意方向不能搞反。MATLAB 里最直接的写法是function X_tf isfft(X_dd, M, N) % 输入 X_dd: M x N 的 DD 域 QAM 符号矩阵 % 输出 X_tf: M x N 的时频域符号矩阵 % 先对每一列做 M 点 IDFT时延-频率 X_tf ifft(X_dd, M, 1); % 再对每一行做 N 点 DFT多普勒-时间 X_tf fft(X_tf, N, 2); % 归一化因子按能量保持处理 X_tf X_tf * sqrt(N / M); end逻辑说明第一维列方向是时延轴经 IDFT 从时延域变换到频率域第二维行方向是多普勒轴经 DFT 从多普勒域变换到时间域。归一化因子保证变换前后能量一致否则符号幅度会随 M、N 变化直接影响后续 PAPR 计算的绝对值。这里容易犯的错是把两个变换方向写反结果是星座图完全散开、误码率平台下不来而且从频谱上看不到明显的带外辐射变化。第二步是海森堡变换把时频域符号矩阵转换成时域发送波形。标准实现是用一个 M 点的 OFDM 调制器逐列处理每列对应一个 OTFS 符号周期内的所有子载波function s heisenberg_transform(X_tf, M, N, cp_len) % 输入 X_tf: M x N 时频域符号矩阵 % cp_len: 循环前缀长度采样点 % 输出 s: 完整的 OTFS 时域发送信号 s []; for n 1:N % 第 n 个 OFDM 符号M 点 IFFT ofdm_symbol ifft(X_tf(:, n), M); % 加循环前缀 ofdm_symbol_cp [ofdm_symbol(end-cp_len1:end); ofdm_symbol]; s [s; ofdm_symbol_cp]; end % 这里是列向量输出s 的长度 N * (M cp_len) end逻辑说明gaussianFilter.m在真实发射机里通常承担脉冲成形的作用即对每个子载波上的符号做时域加窗抑制频谱旁瓣。代码包里的gaussianFilter.m如果按典型实现是生成一个高斯窗系数序列然后对 IFFT 输出做卷积或逐点相乘。实际仿真里是否启用这个滤波器对 PAPR 的影响不小不加窗时矩形脉冲的频谱旁瓣高PAPR 的统计分布相对集中加窗后时域信号包络被平滑峰值被削掉一部分CCDF 曲线的尾部会明显下降。在跑 PAPR 对比实验时建议把加窗与否作为独立变量控制。2.3 接收端维纳滤波与匹配滤波解调OTFS 接收端的解调不是简单取 FFT 就完事。OTFS_demodulation.m里第一个动作是做匹配滤波即在时域上对接收信号做与发射脉冲匹配的相关运算。如果发射端用了高斯滤波器接收端必须用同一个滤波器做匹配否则等效信噪比会下降。实现上可以用一发一收两个滤波器系数做卷积后抽取function Y_tf matched_filter_demod(rx_signal, M, N, cp_len, tx_filter) % rx_signal: 接收时域信号 % tx_filter: 发射端使用的滤波器系数列向量 % 输出 Y_tf: M x N 时频域接收符号 Y_tf zeros(M, N); offset 0; for n 1:N % 去掉循环前缀 symbol rx_signal(offsetcp_len1 : offsetcp_lenM); % 匹配滤波频域相乘等效于时域卷积 symbol_f fft(symbol, M); filter_f conj(fft(tx_filter, M)); symbol_filtered ifft(symbol_f .* filter_f, M); Y_tf(:, n) symbol_filtered; offset offset M cp_len; end endfunction X_dd_hat symplectic_fft(Y_tf, M, N) % 输出: M x N 的 DD 域估计符号矩阵 % 先对每一行做 N 点 IDFT时间-多普勒 Y_delay ifft(Y_tf, N, 2); % 再对每一列做 M 点 DFT频率-时延 X_dd_hat fft(Y_delay, M, 1); X_dd_hat X_dd_hat * sqrt(M / N); end逻辑说明这里先逐列做 N 点 IDFT 把时间维换成多普勒维再逐行做 M 点 DFT 把频率维换成时延维与发射端 ISFFT 的方向正好对偶。匹配滤波的意义在于最大化每个采样点的信噪比减少后续消息传递检测器输入的噪声方差估计误差。代码中容易踩的坑是如果发射端没有启用高斯滤波器但接收端仍然调用了filter_f相当于额外引入了码间干扰BER 曲线会出现错误平层而不是缓慢下降。如果你只想验证 OTFS 基础链路可以把滤波器系数设为全 1 的矩形窗如果要做实际的 PAPR 优化对比再把高斯窗系数提上来。3. 时变多径信道建模与 MP 检测器的关键实现3.1 离散时延-多普勒信道生成原理OTFS_channel_gen.m负责生成 DD 域信道矩阵。物理上高速移动环境下的多径信道可以建模为 P 条路径的叠加每条路径有自己的复增益、时延和多普勒频移。OTFS 的巧妙之处在于这些路径在 DD 域里表现为集中在若干格点上的冲激而非时频域中随时间连续变化的复杂响应。仿真的标准做法是先在连续域定义路径参数再换算到离散网格function [h_dd, channel_params] otfs_channel_gen(M, N, P, delay_spread, doppler_spread, fs) % M: 子载波数时延维 % N: 符号数多普勒维 % P: 多径数 % delay_spread: 最大时延扩展秒 % doppler_spread: 最大多普勒扩展Hz % fs: 采样率 % 输出 h_dd: M x N 的 DD 域信道响应矩阵 delta_f fs / M; T M / fs; h_dd zeros(M, N); for p 1:P % 每条路径的时延在 [0, delay_spread] 内均匀分布 tau_p rand * delay_spread; % 多普勒在 [-doppler_spread, doppler_spread] 内均匀分布 nu_p (2 * rand - 1) * doppler_spread; % 复增益瑞利衰落模型 g_p (randn 1i*randn) / sqrt(2); % 换算到离散索引时延-采样点索引多普勒-多普勒分辨率索引 l_p round(tau_p * fs); k_p round(nu_p * N * T); % 越界保护 if l_p M, l_p M - 1; end if abs(k_p) N/2, k_p sign(k_p) * (N/2 - 1); end h_dd(l_p1, mod(k_p, N)1) h_dd(l_p1, mod(k_p, N)1) g_p; end end参数说明表参数典型值含义与调整方向M64~256时延维分辨率子载波数增大可分辨更小时延但 PAPR 统计特性变化N8~32多普勒维分辨率符号数增大可测更小多普勒步进P4~12多径数越大信道越接近连续散射MP 检测收敛越慢delay_spread1~5μs城市宏站典型值超过 M/fs 会折叠doppler_spread500~2000Hz对应移动速度超过 Δf/2 会混叠需要注意mod(k_p, N)1这种循环移位写法是为了把多普勒频移映射到 N 点 DFT 的周期性输出上。这在实现上是对的但意味着仿真结果只在多普勒扩展不超过 Δf/2 时物理有效。代码包里的OTFS_channel_output.m负责把发射信号和信道响应做二维循环卷积。二维卷积直接用conv2是不对的因为 OTFS 帧结构要求在时延维和多普勒维上都是循环的标准做法是用 FFT 加速function r_dd otfs_channel_output(X_dd, h_dd) % X_dd: M x N DD 域发送符号矩阵 % h_dd: M x N DD 域信道响应矩阵 % 二维循环卷积通过两次 FFT 实现 X_f fft2(X_dd); H_f fft2(h_dd); Y_f X_f .* H_f; r_dd ifft2(Y_f); end逻辑说明时延维度的循环卷积对应 OFDM 循环前缀要覆盖最大时延扩展多普勒维度的循环卷积则依赖于 OTFS 符号块的整体设计。如果 N 太小导致多普勒分辨率不够多普勒域的循环卷积会引入帧间干扰表现为检测器输出的残余误差无法随迭代下降。3.2 消息传递检测器的迭代逻辑与初始化MP 检测器的输入是接收到的 DD 域符号矩阵和信道矩阵核心思路是把 M×N 个符号的联合检测分解为逐符号的置信度传播。OTFS_mp_detector.m的实现框架是一个迭代的干扰消除过程第一步先算每个观测节点对每个符号的似然第二步更新符号节点的概率质量函数第三步做硬判决。关键代码结构如下function [X_hat, iter_err] otfs_mp_detector(Y_dd, h_dd, constellation, noise_var, max_iter) [M, N] size(Y_dd); % 将星座点展开成列向量 const_points constellation(:).; Q length(const_points); % 初始化符号概率均匀分布 prob ones(M*N, Q) / Q; % 用信道能量做归一化 H_energy abs(h_dd).^2; noise_var_eff noise_var * mean(H_energy(:) 1); X_hat zeros(M, N); for iter 1:max_iter % 计算期望符号和方差 sym_exp prob * const_points.; sym_var prob * (abs(const_points).^2). - abs(sym_exp).^2; % 干扰消除用上一次的符号估计重构接收信号并减去 interference h_dd .* reshape(sym_exp, M, N); residual Y_dd - interference h_dd .* reshape(sym_exp, M, N); % 计算每个星座点的似然 likelihood zeros(M*N, Q); for q 1:Q diff abs(residual - h_dd * const_points(q)).^2; likelihood(:, q) -diff(:) / noise_var_eff; end % 更新概率并归一化 prob exp(likelihood - max(likelihood, [], 2)); prob prob ./ sum(prob, 2); % 硬判决 [~, idx] max(prob, [], 2); X_hat reshape(const_points(idx), M, N); iter_err(iter) mean(abs(X_hat(:) - sym_exp(:)).^2); if iter 1 abs(iter_err(iter) - iter_err(iter-1)) 1e-6 break; end end end逻辑说明这个实现里每个符号节点只连接一个观测节点因为在 DD 域信道稀疏矩阵中每个格点最多被少数几条路径影响所以消息传递退化成了迭代的逐符号最大后验估计。干扰消除那行residual Y_dd - interference h_dd .* reshape(sym_exp,...)的写法其实是在复原“仅当前符号有价值的接收分量”即从总接收信号中减去其他所有符号的贡献。noise_var_eff里加上mean(H_energy)是个工程近似把信道估计误差的方差也折算进了噪声防止高信噪比下概率更新过冲。迭代收敛后X_hat就是恢复出来的 DD 域符号矩阵直接做 QAM 解映射即可。如果 max_iter 设成 5 以下BER 性能会明显劣于线性均衡器一般 10~20 次迭代足够收敛。代码里的iter_err向量可以用来画收敛曲线观察算法在几轮后进入稳定状态这是判断信道条件好坏的一个辅助窗口。4. PAPR 计算与降低CCDF 曲线绘制和 PTS 方法落地4.1 从时域波形计算 PAPR 的完整流程PAPR 的定义是时域信号的峰值功率与平均功率之比。OTFS 的 PAPR 计算和 OFDM 有相似之处但有一个关键差异OTFS 的一个帧包含 N 个 OFDM 符号计算 PAPR 时要把整个帧的时域信号连在一起统计还是对每个符号分别统计再取平均会得出不同结论。学术界更常见的是对整个帧的时域波形求 CCDF因为实际功率放大器的峰值限制是作用于连续信号流的。MATLAB 中标准的计算流程如下function papr_db calculate_papr(s) % s: 时域发送信号列向量 % 返回 PAPR 值单位 dB peak_power max(abs(s).^2); avg_power mean(abs(s).^2); papr_db 10 * log10(peak_power / avg_power); end但跑单次仿真的 PAPR 值没有统计意义正确的做法是蒙特卡洛多次仿真取 1e4~1e5 个 OTFS 帧对每个帧算出一个 PAPR 值然后画 CCDF 曲线互补累计分布函数即 PAPR 超过某个门限的概率。代码框架如下function ccdf_curve generate_ccdf(M, N, QAM_order, num_frames) papr_all zeros(num_frames, 1); for idx 1:num_frames % 生成随机 QAM 符号 data randi([0 QAM_order-1], M*N, 1); syms qammod(data, QAM_order, UnitAveragePower, true); X_dd reshape(syms, M, N); % 调制链路省略滤波器 X_tf isfft(X_dd, M, N); s heisenberg_transform(X_tf, M, N, cp_len); papr_all(idx) calculate_papr(s); end % 按门限统计 thr_db 0:0.1:15; ccdf_curve zeros(size(thr_db)); for i 1:length(thr_db) ccdf_curve(i) sum(papr_all thr_db(i)) / num_frames; end semilogy(thr_db, ccdf_curve); xlabel(PAPR threshold (dB)); ylabel(CCDF); grid on; end参数说明QAM_order16时4bit/符号未做任何处理的 OTFS 信号 CCDF 曲线在 1e-3 概率处的 PAPR 大约在 10.8dB 附近QAM 阶数升高时 PAPR 略微增大。num_frames少于 1e4 时曲线的尾部1e-3 以下晃动很大不够严谨。画 CCDF 时注意纵轴用对数坐标横轴门限从 0dB 开始10dB 处截止比较合理。高斯滤波器开启后时域波形被平滑PAPR 会下降 0.5~1dB这对放大器效率是个可感知的改善。4.2 部分传输序列的 MATLAB 实现与复杂度权衡PTSPartial Transmit Sequence是降低 OFDM/OTFS PAPR 的经典手段思路是把时频域符号矩阵按子载波或者按符号分成 V 个子块每个子块乘以一个相位旋转因子然后搜索一组相位使叠加后的信号峰值最小。PTS 在 OTFS 上的一个直接做法是在 ISFFT 之前对 DD 域符号矩阵做分块旋转而不是在时频域分块。这样做的原因是避免改写海森堡变换的核心结构只改动发射端预处理部分。实现如下function s_pts otfs_pts(X_dd, V, phase_set, M, N) % X_dd: M x N DD 域符号矩阵 % V: 子块数V 2 表示分成上下两块V 4 表示分成四块 % phase_set: 可选相位集合例如 [1, -1, 1i, -1i] % 按行分成 V 个子块 block_size floor(M / V); sub_blocks cell(V, 1); for v 1:V row_idx ((v-1)*block_size1) : (v*block_size); sub_blocks{v} X_dd(row_idx, :); end % 搜索最优相位组合 num_phases length(phase_set); best_phases ones(V, 1) * phase_set(1); best_papr inf; % 遍历所有相位组合注意实际中可以用迭代搜索降低复杂度 for b 0:(num_phases^V - 1) phases_idx dec2base(b, num_phases, V) - 0 1; phases phase_set(phases_idx); % 重建加权的 DD 域矩阵 X_pts zeros(M, N); for v 1:V row_idx ((v-1)*block_size1) : (v*block_size); X_pts(row_idx, :) phases(v) * sub_blocks{v}; end % 走一遍调制并计算 PAPR X_tf_tmp isfft(X_pts, M, N); s_tmp heisenberg_transform(X_tf_tmp, M, N, cp_len); papr_tmp calculate_papr(s_tmp); if papr_tmp best_papr best_papr papr_tmp; best_phases phases; end end % 用最优相位生成最终输出 X_pts zeros(M, N); for v 1:V row_idx ((v-1)*block_size1) : (v*block_size); X_pts(row_idx, :) best_phases(v) * sub_blocks{v}; end s_pts heisenberg_transform(isfft(X_pts, M, N), M, N, cp_len); end逻辑说明dec2base那行是把 0 到 num_phases^V-1 的整数转换成 V 位的 num_phases 进制序列用来遍历所有相位组合。这段代码的问题是复杂度随 V 和相位数量呈指数增长V4、相位集合大小为 4 时就要遍历 256 种组合每个组合做一次完整的 OTFS 调制仿真时间会拖很长。工程上常见做法是改成随机搜索或迭代贪心即固定其他子块相位逐个优化单个子块效果接近穷举但复杂度大幅下降。文中给出的穷举版本逻辑清晰作为验证基线合适实际跑大规模仿真时可以换成贪心版本。PTS 对 OTFS 的 PAPR 抑制增益受分块方式影响。分在时延维行方向时每块覆盖全部多普勒范围符号相关性较高PTS 增益约为 1.5~2.5dB1e-3 概率点分在多普勒维则增益略低。另一个细节是接收端必须知道发射端选的相位组合才能正确解调所以 PTS 需要一个边信息信道在实际系统设计中要预留比特开销。4.3 样值钳位与峰值加窗的补充实现PTS 是概率类方法不保证每个帧的峰值都被压低到目标门限之下。想要严格限制峰值就得用确定性方法比如样值钳位Clipping或峰值加窗Peak Windowing。OTFS_PAPR相关代码里如果包含硬钳位实现其逻辑通常是在时域信号上做非线性映射function s_clipped peak_clipping(s, papr_target_db) avg_power mean(abs(s).^2); peak_limit sqrt(avg_power * 10^(papr_target_db / 10)); s_clipped s; exceed_idx abs(s) peak_limit; s_clipped(exceed_idx) peak_limit * exp(1i * angle(s(exceed_idx))); end这里只压缩幅度、不改相位避免额外的相位失真。钳位比Clipping Ratio, CR定义为限幅门限与均方根幅度之比CR 越低 PAPR 越低但带内失真增大、误码率上升。CR1.4约 3dB 回退时16QAM 的误码率损失约为 0.5~1dB这个代价对很多场景是可以接受的。更好的方案是做迭代滤波的钳位即钳位后带通滤波再钳位收敛后带外辐射和带内失真都有改善。代码包中如果只给了普通钳位版本建议自行叠加一个designfilt低通滤波器做两轮迭代。5. 把整个链路跑通主脚本执行顺序、调试技巧与参数边界OTFS_sample_code.m是主脚本它的执行顺序决定了仿真结果是否可信。建议的调用顺序是先运行OTFS_channel_gen.m生成信道参数并打印路径信息再做调制、加信道、解调、MP 检测最后计算误码率和 PAPR。这一步看起来基础但很多人会跳过信道可视化直接跑完整链路导致后续 BER 异常时不知道问题出在调制还是信道还是检测器。先把独立的模块输出分别打印和绘制出来每个模块验证无误后再组合。调试时有一个区分问题来源的技巧把信道矩阵设为恒等矩阵h_dd eye(M)此时 OTFS 链路简化为 AWGN 信道检测器输出的 BER 应该与理论 QAM 误码率曲线几乎重合。如果连这一关都过不了说明调制/解调的 FFT 方向或者归一化有误如果能过说明问题出在信道生成或 MP 检测器与信道的配合上。这个简单的消融测试可以在两分钟内定位 80% 的链路问题。关于参数边界有三点值得注意。第一M 和 N 的取值直接影响多普勒分辨率和时延分辨率M64、N16 的参数组合下最大可测多普勒约为 7.5kHz超过这个值后信道矩阵出现多普勒混叠MP 检测器的迭代误差会随信噪比升高而收敛到一个不为零的平台。此时应该增大 Δf 或者减小移动速度而不是盲目增加 MP 迭代次数。第二PTS 的相位搜索空间与调制阶数无关但受 M 影响M 越大分块后每块包含的符号越多PTS 增益越稳定但穷举时间也越长。第三gaussianFilter.m的滤波器长度不要超过 CP 长度否则滤波引入的时延扩展会超出循环前缀的保护范围导致本应被消除的符号间干扰重新出现。如果你要在这个代码包基础上做扩展研究通用路径有两条。一是替换检测器把 MP 换成基于近似消息传递的 OAMP 或者直接调cvx做稀疏重构对比不同算法在高速场景下的 BER 性能和计算复杂度。二是做多用户 OTFS 的 PAPR 分析代码里的 PAPR 计算函数稍加修改就能输出多用户叠加后的 PAPR 分布观察用户数增加时 CCDF 曲线的退化趋势。建议在跑批量的 PAPR 实验前先把蒙特卡洛帧数定在 2 万以上同时把rng(seed)在每次实验前固定下来否则不同方案之间的 PAPR 差异会被随机性淹没画出来的对比曲线很难说服审稿人。本文还有配套的精品资源点击获取