这是 Keystone 变换系列的第四篇。前面几篇我把公式推导和物理图像都讲了一遍，重点解释了为什么运动目标的回波会“跑”出距离单元，以及 Keystone 变换为什么能通过重采样把慢时间轴“掰弯”来校正距离走动。这篇直接落地，用 MATLAB 把整个流程跑通。我会给出一组完整可复现的雷达仿真参数，从回波生成、脉冲压缩到 Keystone 变换、相参积累，全部用代码说话。内容适合做雷达信号处理、在准备相参积累相关的课程设计或者毕业设计的同学，也适合已经从公式层面理解了 Keystone，但还不知道“频域重采样”在 MATLAB 里到底怎么写的人。
1. 仿真之前，必须想明白的三个问题
1.1 距离走动在二维图上到底长什么样
仿真之前先建立直觉。雷达发射一串脉冲，每个脉冲的快时间维对应距离，脉冲之间的慢时间维对应观测时刻。目标如果静止不动，每个脉冲回波经过匹配滤波之后，峰值包络都落在完全相同的距离单元上，画成“距离-慢时间”二维图，就是一条水平的亮线。目标一旦运动，每个脉冲对应的时间延迟 (2R(t_m)/c) 在变，包络峰值就会在一个一个距离单元之间移动，二维图上看起来就是一条倾斜的线。
这条斜线就是距离走动。走动的速度直接和目标径向速度挂钩，走得越快，积累时间内跨过的距离单元就越多。仿真里最容易犯的错，是只盯着波形看而忽略距离单元的定义。距离单元宽度是 (c/(2B))，带宽 100 MHz 时也就是 1.5 米，而一个高速目标在 0.1 秒内可能飞出几十米。换句话说，如果不做任何校正，直接做慢时间维傅里叶变换去测速度，目标能量会被摊到十几个距离单元上，多普勒谱严重展宽，积累增益大幅下降。
1.2 Keystone变换的输入输出到底是谁
很多文章把 Keystone 变换讲得很玄，其实落到仿真里就是一句话：输入端是“脉冲压缩后的二维矩阵”，输出端是“重采样后的二维矩阵”。矩阵的两个维度分别是慢时间和距离频率。注意，这里说的距离频率是指对快时间维做傅里叶变换之后的频率轴，不是慢时间维的多普勒频率。
为什么要到距离频域去做？因为距离走动在时域表现为包络位置的平移，而在频域表现为相位项 ( -4\pi (f_c + f_\tau) v t_m / c )。这个相位项是距离频率 (f_\tau) 和慢时间 (t_m) 的乘积耦合，所以没法用简单的时移补偿。Keystone 的思路很直接：对每个距离频率点，把慢时间轴按比例 (f_c/(f_c+f_\tau)) 重新缩放，让耦合项中的 (f_\tau) 和 (t_m) 不再相乘，距离走动自然就被拉回来了。
1.3 为什么不能在时域直接“对齐峰值”
容易想到的一个替代方案是：既然知道包络是斜线，就估计斜率，然后把每个脉冲平移回来。这种做法叫包络对齐，在 ISAR 成像里很常用，但问题在于它需要先估计目标速度或者包络斜率，而且对多个不同速度的目标处理起来很麻烦。Keystone 变换的好处是完全不需要速度先验，它用重采样操作把“速度-频率”耦合统一消除，不管你是一个目标还是多个目标，只要信噪比不是太差，都能一次性拉直。
2. 仿真参数与信号构造
2.1 一组留有余量的参数表
仿真参数直接决定了距离走动是否明显，也决定了 Keystone 变换能不能体现出效果。我用的这组参数比较典型，覆盖了常见的 X 波段雷达场景。
| 参数 | 数值 | 说明 |
|---|---|---|
| 载频 (f_c) | 10 GHz | X 波段 |
| 带宽 (B) | 100 MHz | 距离分辨率 1.5 m |
| 脉宽 (T_p) | 10 μs | 匹配滤波增益的基础 |
| 采样率 (f_s) | 100 MHz | 取等于带宽，距离窗无过采样 |
| 脉冲重复频率 PRF | 2000 Hz | 无多普勒模糊对应的最大速度约 30 m/s，后面细讲 |
| 积累脉冲数 M | 256 | 积累时间 128 ms |
| 初始距离 R0 | 10000 m | 远场 |
| 目标速度 v | 200 m/s | 高速，明显距离走动 |
计算一下积累时间 (T_{obs} = M / PRF = 128) ms，目标在这段时间里移动约 25.6 m。距离单元宽度只有 1.5 m，所以包络会跨越约 17 个距离单元。如果这条路走通了，Keystone 变换带来的效果在图上会非常直观。
2.2 生成回波信号的MATLAB代码
下面这段代码生成基带回波。为了教学直观，我用了循环逐脉冲生成，实际工程中可以改成矩阵化运算提速。
%% 参数设置 c = 3e8; fc = 10e9; B = 100e6; Tp = 10e-6; K = B / Tp; PRF = 2000; M = 256; dt = 1 / PRF; fs = 100e6; Ts = 1 / fs; Nfast = 1024; R0 = 10000; v = 200; tm = (0:M-1) * dt; R = R0 - v * tm; % 目标距离随慢时间线性变化 %% 快时间轴 t_ref = 2 * R0 / c; % 参考时延 tau = t_ref - Tp/2 + (0:Nfast-1) * Ts; %% 发射基带信号 tx = exp(1j * pi * K * ((tau - t_ref).^2)); % 以参考时延为中心的内向线性调频 tx(abs(tau - t_ref) > Tp/2) = 0; %% 回波生成 s = zeros(M, Nfast); for m = 1:M td = 2 * R(m) / c; t = tau - td; s(m, :) = exp(1j * pi * K * t.^2) .* (abs(t) <= Tp/2); end s = s .* exp(-1j * 4 * pi * fc * R(:) / c); % 慢时间相位，速度信息所在这段代码里最关键的是那个慢时间相位 (\exp(-j4\pi f_c R(t_m)/c))。它包含了目标运动引起的多普勒信息，后面测速度就要靠它。
2.3 脉冲压缩与运动补偿前的基线结果
脉冲压缩用频域匹配滤波实现。构造匹配滤波器的频响时要注意，把发射信号补零到和回波矩阵相同的列数再取共轭。
Nfft = 2^nextpow2(Nfast * 2); S = fft(s, Nfft, 2); ref = exp(1j * pi * K * ((tau - t_ref).^2)); % 参考信号 H = fft(ref, Nfft, 2); Y = ifft(S .* conj(H), [], 2); Y = Y(:, 1:Nfast);这里参考信号取以参考时延为中心，是为了让匹配滤波后的峰值出现在参考距离附近。如果不做任何走动校正，直接查看二维图的幅度，会看到目标峰值包络随慢时间明显偏移，一条斜线贯穿多个距离门。这就是我们后面要校正的对象。
3. 核心实现：MATLAB里的三种Keystone变换写法
3.1 方法一：慢时间维Sinc插值
最符合原理、也最好理解的实现方式是 Sinc 插值。对每个距离频率点 (f_\tau)，计算缩放因子 (\alpha = f_c / (f_c + f_\tau))，然后对慢时间序列在非均匀的 (t_m' = \alpha t_m) 处重新采样。
Yf = fftshift(fft(Y, Nfft, 2), 2); freq = ((-Nfft/2 : Nfft/2-1) / Nfft) * fs; ValidF = abs(freq) <= B/2; Ykt = zeros(size(Yf)); for i = 1:Nfft if ValidF(i) alpha = fc / (fc + freq(i)); tm_new = tm * alpha; Ykt(:, i) = ks_sinc_interp(Yf(:, i), tm, tm_new, dt, 16); else Ykt(:, i) = Yf(:, i); end end y_kt = ifft(ifftshift(Ykt, 2), [], 2); y_kt = y_kt(:, 1:Nfast);自定义 Sinc 插值函数如下。这里的思路是对每个目标采样点，用带限 Sinc 核卷积估计新时刻的值，核长截断到 16 个采样点。
function y = ks_sinc_interp(x, t, t_new, dt, Ntap) % x: Mx1 原始慢时间序列 % t: Mx1 原始慢时间 % t_new: Lx1 新慢时间 % dt: 慢时间采样间隔 % Ntap: 核截断长度 L = length(t_new); M = length(x); y = zeros(L, 1); for n = 1:L tau = (t_new(n) - t) / dt; idx = find(abs(tau) < Ntap); if ~isempty(idx) y(n) = sum(x(idx) .* sinc(tau(idx))); end end end这种写法最贴近公式，适合讲原理和验证算法。缺点是慢时间维每个距离频率都需要做一次插值，256 个脉冲、1024 个距离采样点就要做 1024 次循环，速度比较慢。实际仿真可以先跑通，再考虑优化。
3.2 方法二：Chirp-Z变换实现
工程上更推荐的实现是利用 Chirp-Z 变换。Keystone 重采样本质上是一种尺度变换，而尺度变换可以由“频域相位调制 + 时域卷积”来实现，正好是 Chirp-Z 的典型应用。这样做的好处是不需要显式地构造新时刻网格，也不会有 Sinc 核截断造成的边缘误差，计算速度更快。
MATLAB 自带czt函数，可以对慢时间序列做任意起止频率和点数的 Z 变换。对每个距离频率点，把重采样看成在频域做一次变尺度 Z 变换即可。核心代码框架如下：
for i = 1:Nfft if ValidF(i) alpha = fc / (fc + freq(i)); % 对慢时间序列做 Chirp-Z 变换，等效于重采样到 alpha 倍 temp = czt(Yf(:, i), M, exp(-1j * 2 * pi * alpha / M), 1); Ykt(:, i) = temp(:); end end这里我只是给出最基本的思路。实际使用时要仔细确认czt的输入参数定义，尤其是起点和螺旋因子的设置，否则很容易把尺度关系搞反。由于czt在信号处理工具箱中可用，大部分机器上都能直接跑，速度比逐点 Sinc 插值快一个数量级。
3.3 方法三：基于FFT的等效实现
第三种思路是利用 FFT 在频域完成插值。对慢时间序列先做 FFT 到多普勒域，通过补零和相位修正来近似实现分数阶重采样。这个方法实现起来最快，但精度受补零方式和尺度因子整数化影响较大，适合快速验证或者对精度要求不高的场景。
for i = 1:Nfft if ValidF(i) alpha = fc / (fc + freq(i)); Nout = round(M * alpha); spec = fft(Yf(:, i), M); spec_pad = zeros(Nout, 1); if Nout >= M spec_pad(1:M/2) = spec(1:M/2); spec_pad(end-M/2+1:end) = spec(M/2+1:end); else spec_pad(1:Nout/2) = spec(1:Nout/2); spec_pad(end-Nout/2+1:end) = spec(M-Nout/2+1:end); end Ykt(:, i) = ifft(spec_pad, Nout); % 长度会变化 end end注意，这个方法会导致每个距离频率上的输出长度不一致，后续需要对齐到统一的慢时间网格。实际工程中通常不建议直接这样用，但用来理解“重采样即变尺度频谱搬移”的物理意义很不错。
三种方法的取舍我整理成了表格。
| 方法 | 精度 | 速度 | 实现复杂度 | 场景建议 |
|---|---|---|---|---|
| Sinc插值 | 高 | 慢 | 低 | 原理验证、教学 |
| Chirp-Z | 高 | 快 | 中 | 工程部署、批量处理 |
| FFT近似 | 中 | 最快 | 中 | 快速预览、实时性要求高 |
4. 仿真结果：校正前与校正后的对比
4.1 B-scan图与距离包络对比
把脉冲压缩后的幅度沿慢时间画出来，校正前的 B-scan 图上目标是一条倾斜的亮线，跨过的距离单元数大约为 (2vT_{obs}/c / (1/(2B)) \approx 17) 个。做完 Keystone 变换之后再画，斜线变成接近水平的直线，目标能量被“收拢”到同一距离单元。
这里有个实操细节：画图时最好把显示范围固定在同一动态范围内，比如统一用分贝值20*log10(abs(Y))，否则人眼很容易被颜色映射欺骗，觉得“好像没什么变化”。另外校正后如果 Sinc 插值用了 0 外插，图像边缘可能会出现暗区，这是正常的，不代表算法失效。
4.2 慢时间FFT与速度测量
距离走动校正的最终目的是为了相参积累。校正前，对每个距离门做慢时间 FFT，目标能量分散在多个距离门，多普勒谱被展宽。校正后，目标集中在同一个距离门，对这个距离门做慢时间 FFT 会得到一个尖锐的峰值，峰值频率对应的多普勒频率 (f_d = -2vf_c/c)。
具体到我这组参数：(f_d = -2 \times 200 \times 10^{10} / 3\times 10^8 \approx -13333) Hz。PRF 是 2000 Hz，所以这个频率已经远远超出了 ([-PRF/2, PRF/2]) 的范围，出现多普勒模糊。测出来的模糊频率是 (13333 \mod 2000 = 1333) Hz，对应速度约 (-20) m/s。这说明仿真里 PRF 选择并不适合这么高的速度，真实工程中要么提高 PRF，要么用多普勒模糊数解算。
4.3 信噪比改善和计算代价
评价 Keystone 效果不能只看图，更要看量化指标。简单做法是：在校正前的二维矩阵里，找到目标所在距离门和附近距离门的所有慢时间 FFT 幅度，计算峰值幅度和噪声底之间的比值；然后用同样方法计算校正后同一距离门上的比值。参考仿真条件下，校正后峰值幅度相比校正前通常能提升 10 dB 以上，具体数字取决于距离走动跨越的单元数。
代价方面，用 Sinc 插值跑 256×1024 的数据量，在普通笔记本上可能要几秒到十几秒，而用 Chirp-Z 可以压到 1 秒以内。做科研时建议先用 Sinc 验证正确性，再换成 Chirp-Z 跑蒙特卡洛实验。
5. 仿真过程中最容易踩的坑
5.1 插值核截断导致边缘失真
Sinc 插值理论上需要无限长的核，但实际只能截断。截断会带来两个问题：一是回波边缘产生振铃，二是当缩放因子偏离 1 较多时，新时刻可能落在原始慢时间范围之外，直接外插会把虚假能量带进来。我的处理办法是在插值函数里对超出范围的值强制置 0，只保留中间约 80% 的慢时间数据进行后续测速和多普勒分析。
5.2 快时间频率轴偏移导致重采样比例错误
频域做 Keystone 时，距离频率轴必须对应基带频率，也就是范围是 ([-f_s/2, f_s/2])。很多人直接用fft后的下标当频率，忘记做fftshift，结果频率正负颠倒，重采样比例全部反了，效果自然不对。建议在生成频率轴之后，先用一个直流信号自检，确保频率为 0 的点位于Nfft/2+1的位置。
5.3 多普勒模糊对测速结果的影响
前面已经提到，PRF 太低时 Keystone 变换本身仍然可以校正距离走动，但测速会模糊，因为慢时间采样不满足奈奎斯特条件。遇到这种情况，可以先盲速分割，或者用多 PRF 解模糊。如果只是在仿真里验证 Keystone 本身，建议把目标速度降到 30 m/s 以内，这样 PRF 为 2000 Hz 时不会模糊，结果看起来更干净。
5.4 数据量太大时如何加速
Keystone 最耗时的部分在慢时间维插值。除了 Chirp-Z 之外，还可以用 MATLAB 的interp1并行计算。如果处理的是采信机采回来的数据，可能还涉及 ADC 量化、十六进制转有符号数这些环节，建议先把数据切成长度合适的块，再逐块做 Keystone。实在需要极致性能，就把插值核写成 C 代码，用 MEX 编译调进 MATLAB，这样比纯脚本快很多。
个人体会是，Keystone 变换在 MATLAB 里跑通不难，难的是把边界条件和频率轴关系理清楚。刚开始做仿真时，我第一版代码跑出来的结果是一条完全乱掉的曲线，检查了很久才发现是频率轴没有做fftshift。所以建议你拿到代码后，先跑一组低速目标确认基线正确，再加大速度，这样能快速定位问题。这个系列后面的内容，我打算继续深入聊 Keystone 在多目标场景下的表现，以及和长时间相参积累算法的配合，如果大家在复现过程中遇到其他锯齿，欢迎一起交流。
