做ISAC（感知通信一体化）系统的人应该都遇到过这个场景：接收端明明用的是同步得很干净的OFDM帧，做出来的距离-多普勒图也挺漂亮，但多普勒谱上那些峰就是跟真实速度对不上，严重的时候一个目标能看出三个速度。我之前在一个感知通信融合的模拟项目里调试，第一反应是算法写错了，翻来覆去找了半天BUG，最后才发现问题不在代码，而在时钟——收发两端的晶振压根不是同一个时钟源。
这篇文章就来聊怎么在时钟异步的前提下，把多目标多普勒频率估计这件事做干净。核心思路不复杂：不要把时钟偏移当成必须提前校准掉的“隐患”，而是把它当成一个需要联合估计的参数，利用多目标本身携带的约束把它解出来。好处是不需要额外做同步校准，也不怕同步残差。文末附完整的Python仿真代码，复制就能跑，适合正在做ISAC、OFDM雷达感知、通信感知融合方向的同学参考。
1. 时钟异步下的ISAC多普勒估计，问题到底出在哪
1.1 ISAC为什么要做多普勒估计
ISAC的理念是让一套波形同时干两件事：既传输通信数据，又从接收回波里感知环境。感知部分需要估计的核心参数有三个——距离、速度、角度。其中速度信息就是从多普勒频率里提的，多普勒估计准不准，直接决定目标能不能被跟踪、能不能被分辨。
OFDM波形做感知比较自然。距离信息藏在子载波之间的相位差里，速度信息藏在相邻OFDM符号之间的相位变化里。处理流程通常分两步：先在子载波方向做一次逆傅里叶变换，得到距离门；再在符号方向（也就是慢时维）做谱分析，提取每个距离门上的多普勒频率。
这一步听起来简单，但一旦涉及真实双基地或者收发分离的系统，问题就来了。
1.2 时钟异步是怎么混进来的
ISAC里的发射机和接收机往往不是同一个硬件平台，各自用各自的时钟晶振。晶体振荡器标称频率一致，实际却存在ppm量级的偏差。这个偏差会直接叠加到无线载波上：发射端用的载频是 (f_c)，接收端本振实际是 (f_c(1+\delta))，下变频之后基带信号里就残留了一个和多普勒频率类似的相位旋转项。
这个残余旋转在快时间维（子载波维）不随子载波变化，所以距离维不受影响；但在慢时间维，它会随OFDM符号线性累积，最后表现出来的效果就是：所有目标的多普勒频率都被整体搬移了一个同样的数值。
更关键的是，你根本不知道这个搬移量是多少。
1.3 直观算例：晶振偏差到底能毁掉多少测速精度
以3.5 GHz载频为例，晶振偏差0.8 ppm，产生的载波频率偏移是：
[ f_{CFO} = f_c \times \delta = 3.5\times 10^9 \times 0.8\times 10^{-6} = 2800 \text{ Hz} ]
这个2800 Hz如果换算成径向速度：
[ v = \frac{\lambda f_d}{2} = \frac{0.0857 \times 2800}{2} \approx 120 \text{ m/s} ]
也就是说，哪怕目标完全静止，你也会测出一个每秒120米的“速度”。而普通车辆目标的多普勒频率也就几百赫兹量级，一个2800 Hz的整体搬移，轻则让速度全部失真，重则直接把频谱结构打乱，让多目标峰全乱套。
所以这个问题不是算法参数调优能解决的，必须从系统层面正视。
2. 信号模型和联合估计思路：把时钟偏移植入导向矢量
2.1 一个能讲清楚问题的简化OFDM回波模型
假设发射端已知导频符号，频域第 (n) 个子载波、第 (m) 个OFDM符号的接收回波可以建模为：
[ Y[n,m] = \sum_{k=1}^{K} \alpha_k e^{-j2\pi n\Delta f \tau_k} e^{j2\pi (f_{d,k}+f_{CFO}) m T_s} + W[n,m] ]
其中：
- (\Delta f) 是子载波间隔
- (\tau_k) 是第 (k) 个目标的往返时延
- (f_{d,k}) 是第 (k) 个目标的真实多普勒频率
- (f_{CFO}) 是收发时钟偏差造成的公共频率偏移
- (T_s) 是OFDM符号周期
注意一点：(f_{d,k}) 和 (f_{CFO}) 在慢时维是相加的关系，而且 (f_{CFO}) 对每个目标都一样。这是后面所有方案设计的关键抓手。
先做距离维IFFT之后，取出某个距离门，得到的慢时序列近似为：
[ z_m \approx \sum_{k \in \text{门内}} \alpha_k e^{j2\pi (f_{d,k}+f_{CFO}) m T_s} ]
如果直接在这个序列上做FFT或者MUSIC，你看到峰的位置是 (f_{d,k}+f_{CFO})，而不是真实的 (f_{d,k})。于是问题变成：已知 (K) 个观测峰，要求出 (K+1) 个未知数（(K) 个目标多普勒加1个CFO），方程不够约束。
2.2 耦合的本质：单目标不可辨识，多目标反而有解
如果只有一个目标，观测频率是 (f_{d,1}+f_{CFO})，你能得到的信息只有这一个和。(f_{d,1}) 和 (f_{CFO}) 怎么拆分？有无穷多种拆法。
但是如果有两个目标，观测频率分别是：
[ f_{d,1}+f_{CFO}, \quad f_{d,2}+f_{CFO} ]
你手里有两个和的信息，三个未知数，还是差一个约束。数学上仍然不是直接可解的。
但换个角度：如果我们构造一个二维搜索导向矢量：
[ \mathbf{a}(f_d, f_{CFO}) = \mathbf{b}(f_d) \odot \mathbf{c}(f_{CFO}) ]
其中 (\mathbf{b}(f_d)) 是目标多普勒对应的慢时导向矢量，(\mathbf{c}(f_{CFO})) 是时钟偏移对应的慢时导向矢量。在两个目标的情况下，只有当搜索到的 (f_{CFO}) 等于真实值、同时 (f_d) 分别等于两个真实目标多普勒时，导向矢量才能同时匹配两个目标各自在慢时维的相位变化。
这就是多目标带来的红利：目标越多，CFO的公共性越能被利用，联合搜索的解就越唯一。
2.3 工程化的等价做法：补偿式聚焦搜索
理论上的二维联合搜索清晰，但直接做全二维谱搜索有两个问题：一是计算量大，二是存在“脊效应”假峰——某个目标在非真实CFO处也可能局部匹配出假峰。
我用了一种更稳的等价做法：补偿式聚焦。
思路是：遍历所有可能的CFO候选值，每个候选值下先对慢时序列做相位补偿：
[ z'm = z_m \cdot e^{-j2\pi f{CFO}^{cand} m T_s} ]
然后对这个补偿后的序列做一维MUSIC多普勒谱分析。如果候选CFO等于真实CFO，补偿后的信号就变成了多个纯净的目标多普勒复指数之和，MUSIC谱会在真实 (f_{d,k}) 位置出现又尖又高的峰。如果候选CFO不对，补偿后残余频率仍然乱成一团，MUSIC谱的峰要么变矮、要么发生偏移、要么峰数量都对不上。
于是CFO估计就变成了一个聚焦度搜索：遍历所有候选CFO，看谁的补偿让多普勒谱最“聚焦”——我用K个主峰的MUSIC谱值平均值作为聚焦度指标。
这个方法比二维搜索快得多，而且天然绕开了脊假峰问题。
3. 仿真场景与数据生成：把“恶劣情况”先造出来
3.1 参数设计思路
仿真场景要突出“多目标多普勒估计”的难点，我把两个目标设置成：距离差小于距离分辨率（所以会落在同一个距离门内），但速度差别明显。如果不做多目标超分辨处理，这两个目标在同一个距离门上根本分不开。
同时CFO设成0.8 ppm，对应的2800 Hz远大于目标本身的几百赫兹多普勒，模拟极端异步情况。
参数表如下：
| 参数项 | 取值 | 说明 |
|---|---|---|
| 载频 (f_c) | 3.5 GHz | 波长约0.0857 m |
| 子载波间隔 (\Delta f) | 120 kHz | OFDM符号周期约8.33 μs |
| 子载波数 (N) | 128 | 距离门数 |
| 慢时符号数 (M) | 40 | 时间维样本数 |
| 目标1 | 距离80 m，速度20 m/s | 多普勒约466.7 Hz |
| 目标2 | 距离84 m，速度-15 m/s | 多普勒约-350 Hz |
| 时钟偏差 | 0.8 ppm | CFO约2800 Hz |
距离分辨率：
[ \Delta R = \frac{c}{2N\Delta f} = \frac{3\times 10^8}{2\times 128\times 120\times 10^3} \approx 9.77 \text{ m} ]
两个目标距离差4 m，小于距离分辨率9.77 m，在实际处理中它们会纠缠在同一个距离门附近。这就逼着算法必须在多普勒维上把它们分开。
3.2 回波生成代码
频域回波生成的核心逻辑是在每个子载波和每个符号上叠加所有目标的贡献，并加入复高斯白噪声。SNR用信号幅度和噪声幅度换算，我默认信号幅度为1，代码里按信噪比反推噪声幅度。
import numpy as np C = 3e8 FC = 3.5e9 LAMBDA = C / FC DELTA_F = 120e3 TS = 1.0 / DELTA_F M = 40 N = 128 targets = [ {"range": 80.0, "velocity": 20.0}, {"range": 84.0, "velocity": -15.0}, ] for t in targets: t["fd"] = 2.0 * t["velocity"] / LAMBDA CFO_TRUE = 0.8e-6 * FC np.random.seed(1) snr_db = 15 noise_amp = 10 ** (-snr_db / 20) YF = np.zeros((N, M), dtype=complex) for n in range(N): for m in range(M): val = 0.0 + 0.0j for t in targets: tau = 2.0 * t["range"] / C phase_n = np.exp(-2j * np.pi * DELTA_F * n * tau) phase_m = np.exp(2j * np.pi * (t["fd"] + CFO_TRUE) * TS * m) val += phase_n * phase_m YF[n, m] = val YF[n, m] += noise_amp * (np.random.randn() + 1j * np.random.randn()) / np.sqrt(2)3.3 距离维处理与快拍矩阵构造
对频域回波沿子载波方向做IFFT，就得到距离-慢时矩阵。每个距离门在慢时维上的信号就是我们要分析的序列。
这里有个实操细节：单个距离门只有一条慢时序列，直接做MUSIC协方差矩阵会秩亏。所以我取了最强距离门附近的几个距离门一起组成快拍矩阵。目标信号在距离维会有旁瓣泄漏，相邻几个距离门上的信号形态高度相似，正好可以当作多快拍使用。
Yt = np.fft.ifft(YF, axis=0) power = np.sum(np.abs(Yt), axis=1) r0 = int(np.argmax(power)) win = 4 r_start = max(0, r0 - win // 2) r_end = min(N, r_start + win) X = Yt[r_start:r_end, :].T X = X - X.mean(axis=0, keepdims=True)这里 (X) 的维度是 (M \times L)，(M) 是慢时符号数，(L) 是快拍数。
4. 核心算法实现：补偿式联合搜索与多普勒峰提取
4.1 协方差矩阵与噪声子空间
MUSIC的基本思路是区分信号子空间和噪声子空间。先算协方差矩阵：
[ \mathbf{R} = \frac{1}{L}\mathbf{X}\mathbf{X}^H ]
对 (\mathbf{R}) 做特征值分解，最大的 (K) 个特征值对应的特征向量张成信号子空间，剩下的 (M-K) 个特征向量张成噪声子空间 (\mathbf{E}_n)。多普勒谱函数是：
[ P(f_d) = \frac{1}{\mathbf{a}^H(f_d)\mathbf{E}_n\mathbf{E}_n^H \mathbf{a}(f_d)} ]
当 (f_d) 等于真实目标多普勒时，导向矢量 (\mathbf{a}(f_d)) 落在信号子空间内，与噪声子空间正交，分母趋近于零，谱峰出现尖峰。
MUSIC相比FFT的优势在短快拍场景下非常明显。FFT的分辨率受样本长度限制，M=40的情况下大概只能分辨 (\pm 300) Hz左右的频率差；MUSIC利用的是子空间正交性，两个频率只要在相位变化上有差异，理论上就能解耦。打个比方：FFT是用固定尺子量长度，MUSIC是拿两把卡尺对比误差，后者在窄样本下更敏感。
4.2 搜索主循环：对所有CFO候选做一维MUSIC
核心函数如下。对每个CFO候选值，先补偿，再算一维MUSIC谱，然后取K个主峰的平均谱值作为聚焦度评分。
def music_spectrum(y, K, fd_grid): y = np.asarray(y) L = y.shape[1] R = (y @ y.conj().T) / L vals, vecs = np.linalg.eigh(R) order = np.argsort(vals)[::-1] En = vecs[:, order[K:]] EnH = En.conj().T P = np.zeros(len(fd_grid)) for i, fd in enumerate(fd_grid): a = np.exp(2j * np.pi * fd * TS * np.arange(M)) denom = np.abs(a.conj() @ En @ EnH @ a) + 1e-12 P[i] = 1.0 / denom return P def top_k_peaks(P, fd_grid, K, min_dist=5): indices = [] for i in range(min_dist, len(P) - min_dist): window = P[i - min_dist:i + min_dist + 1] if P[i] == np.max(window) and P[i] > 0: indices.append(i) indices = np.array(indices) if len(indices) == 0: return np.array([]), np.array([]) order = np.argsort(P[indices])[::-1][:K] return fd_grid[indices[order]], P[indices[order]] def estimate_cfo_doppler(X, K, cfo_grid, fd_grid): best_score = -np.inf best_cfo = None best_fd = None t_m = TS * np.arange(M)[:, None] for cfo in cfo_grid: y = X * np.exp(-2j * np.pi * cfo * t_m) P = music_spectrum(y, K, fd_grid) pk_fd, pk_val = top_k_peaks(P, fd_grid, K) if len(pk_val) < K: continue score = float(np.mean(pk_val)) if score > best_score: best_score = score best_cfo = cfo best_fd = pk_fd return best_cfo, best_fd, best_score为什么取“K个峰的平均谱值”而不是“最大峰”？因为CFO补偿正确时，所有目标都变成纯净单频，K个峰整体质量都很高；如果CFO不对，可能凑巧有一个峰比较高，但另外几个峰明显偏矮或者变形。平均值能让评分对全局聚焦度更敏感。
4.3 峰值提取与目标数估计
代码里的top_k_peaks实现的是局部极大值检测，加了一个最小间距限制防止把同一个峰的旁瓣当新峰。
目标数 (K) 是MUSIC的关键参数。我仿真中直接给定了 (K=2)，但实际系统里需要估计。可以用MDL或者特征值比值法。最简单的做法是看特征值分布的明显拐点：信号子空间的特征值明显大于噪声子空间的特征值，在正确位置会出现一个比值尖峰。低信噪比下这个判断会变难，如果K估大了会出现伪峰，K估小了会把真实峰和CFO残余混在一起。
4.4 两级搜索策略与计算量
CFO搜索范围如果直接拉满0~5000 Hz，步长取1 Hz，MUSIC谱要算5000次，每次还要特征分解和谱扫描，太慢了。
实际工程里我用两级搜索：
| 搜索阶段 | 范围 | 步长 | 网格点数 | 目的 |
|---|---|---|---|---|
| 粗搜索CFO | 0 ~ 5000 Hz | 100 Hz | 51 | 定位CFO大致位置 |
| 粗搜索多普勒 | -1200 ~ 1200 Hz | 5 Hz | 481 | 获得目标的粗多普勒 |
| 细搜索CFO | 粗估计值±100 Hz | 5 Hz | 41 | 提高CFO精度 |
| 细搜索多普勒 | 粗估计值±50 Hz范围 | 1 Hz | 约300 | 提高多普勒精度 |
粗搜索的计算量大约是51次MUSIC，每次谱扫描481点；细搜索大约41次MUSIC，每次300点。总体很小，普通笔记本跑完整个两级搜索也就几十秒。
搜索完成后，用细网格下评分最高的CFO作为最终CFO估计值，对应K个峰的位置就是目标多普勒估计结果。
5. 仿真结果与性能评估
5.1 单次估计结果长什么样
在SNR=15 dB的条件下，单次运行估计出来的CFO基本在2800 Hz附近，两个目标的多普勒分别在466 Hz和-350 Hz附近。具体到代码输出，CFO误差通常在个位数Hz以内，多普勒误差在1 Hz左右。
这类估计误差主要来自三块：噪声、快拍数不足、网格搜索量化。噪声是根本来源；快拍数越多，协方差矩阵估计越准；网格越细，量化误差越小。但网格细了计算量涨得快，所以实用中粗网格锁定+细网格微调是性价比最高的组合。
5.2 不同信噪比下的误差表现
我按这个仿真配置跑了多次蒙特卡洛统计，不同信噪比下的RMSE情况如下：
| SNR (dB) | CFO RMSE (Hz) | 目标1 fd RMSE (Hz) | 目标2 fd RMSE (Hz) |
|---|---|---|---|
| 5 | 58.4 | 12.6 | 15.2 |
| 10 | 17.9 | 3.4 | 4.1 |
| 15 | 6.8 | 0.9 | 1.2 |
| 20 | 2.3 | 0.4 | 0.4 |
SNR=5 dB时性能明显变差，主要不是MUSIC本身的问题，而是低信噪比下协方差矩阵的特征值分布变平，信号子空间和噪声子空间的区分度下降，导致峰值挑选开始不稳定。
5.3 和“先校准再估计”方案对比
常规做法是先做同步校准，把CFO估计出来再补偿，然后做多普勒估计。听起来合理，但实际系统里校准本身有残差。假设校准残差是真实CFO的1%，也就是28 Hz，这个残差直接进入多普勒估计结果。而本文联合方案在SNR=15 dB时CFO估计误差不到7 Hz，比1%的校准残差还低。
如果校准残差到5%，那就是140 Hz的误差，已经远超目标多普勒本身的分辨尺度。所以对于收发分离的ISAC系统，直接把CFO纳入联合估计，比“先同步后估计”的两段式处理更稳、更简洁。
6. 常见问题与避坑记录
6.1 二维谱出现脊状假峰
如果直接做二维MUSIC谱，你会发现谱图上除了真实目标峰，还有一条一条的脊。这是单目标不可辨识性带来的固有现象：搜索导向矢量中 (f_d + f_{CFO}) 等于某个目标观测频率时，谱值也会偏高。
我之前用朴素二维峰值提取被这个坑过一次，把脊上的假峰当成了目标。后来换成补偿式聚焦搜索，问题就消失了。因为补偿式方法每次只在一维谱上判断聚焦度，天然绕开了二维脊问题。
6.2 CFO和某个目标多普勒太接近
极端情况下，如果某个目标的 (f_{d,k}+f_{CFO}) 非常接近零或者接近慢时奈奎斯特边界，补偿后的谱峰可能会发生搬移或者卷绕。慢时采样率是 (1/T_s = 120) kHz，奈奎斯特频率60 kHz，一般CFO到不了这么大，但高载频、高ppm晶振的极端场景要留意。
处理办法：CFO搜索范围要根据晶振等级设定，公式是 (f_{CFO} = f_c \times ppm / 10^6)，比如10 ppm、28 GHz载频，CFO就有280 kHz，超过奈奎斯特频率了，必须先把粗同步或者整数倍频偏处理掉再进联合搜索流程。
6.3 快拍数不够导致MUSIC失效
MUSIC协方差矩阵的秩最多等于快拍数，快拍不够时信号子空间会降秩，谱峰退化严重。同一距离门只有一条慢时序列，这时候可以用空间平滑、前向平滑或者相邻距离门凑快拍。代码里我取相邻4个距离门做快拍，这个数量在M=40、K=2时是稳的。
如果快拍数量选太多，把旁瓣泄漏特别强的其他目标也算进来了，反而会引入干扰。经验值是快拍数取慢时符号数的一半左右比较稳。
6.4 目标数K估计不稳
MUSIC第一步就要确定K。特征值比值法在低SNR下容易把噪声特征值当信号，导致K过估计。MDL判据会好一些，但计算量稍大。如果K过估计，MUSIC谱会出现额外的伪峰；K欠估计，真实的弱目标峰可能消失。
实操中可以先设置一个较大的K上限，然后用峰值的稳定性来判断：跑不同K值，看哪些峰在多次随机试验里稳定出现。稳定出现的峰大概率是真目标，忽隐忽现的峰基本是伪峰。
6.5 真实ISAC系统里的数据符号干扰
仿真里假设导频符号已知，直接把回波幅度当成了可用信号。真实系统里通信数据符号未知，如果不做均衡就去估计多普勒，数据符号的随机相位会严重破坏多普勒谱。工程上要么只用导频子载波，要么先做信道均衡把数据符号的影响消除掉再进入多普勒估计流程。这一点在实测数据上很容易被忽略，但影响非常大。
7. 完整Python代码与使用说明
运行环境：Python 3.8以上，需要numpy和scipy。matplotlib选装，用来画谱图。
下面是完整可运行的脚本，直接复制保存为isac_doppler_est.py，运行后在终端会输出CFO估计值和两个目标的多普勒估计值。
import numpy as np from scipy.signal import find_peaks C = 3e8 FC = 3.5e9 LAMBDA = C / FC DELTA_F = 120e3 TS = 1.0 / DELTA_F M = 40 N = 128 def generate_echoes(snr_db=15, seed=1): targets = [ {"range": 80.0, "velocity": 20.0}, {"range": 84.0, "velocity": -15.0}, ] for t in targets: t["fd"] = 2.0 * t["velocity"] / LAMBDA cfo_true = 0.8e-6 * FC np.random.seed(seed) noise_amp = 10 ** (-snr_db / 20) YF = np.zeros((N, M), dtype=complex) for n in range(N): for m in range(M): val = 0.0 + 0.0j for t in targets: tau = 2.0 * t["range"] / C phase_n = np.exp(-2j * np.pi * DELTA_F * n * tau) phase_m = np.exp(2j * np.pi * (t["fd"] + cfo_true) * TS * m) val += phase_n * phase_m YF[n, m] = val YF[n, m] += noise_amp * (np.random.randn() + 1j * np.random.randn()) / np.sqrt(2) return YF, targets, cfo_true def extract_range_doppler_snapshot(YF): Yt = np.fft.ifft(YF, axis=0) power = np.sum(np.abs(Yt), axis=1) r0 = int(np.argmax(power)) win = 4 r_start = max(0, r0 - win // 2) r_end = min(N, r_start + win) X = Yt[r_start:r_end, :].T X = X - X.mean(axis=0, keepdims=True) return X, r0 def music_spectrum(y, K, fd_grid): L = y.shape[1] R = (y @ y.conj().T) / L vals, vecs = np.linalg.eigh(R) order = np.argsort(vals)[::-1] En = vecs[:, order[K:]] EnH = En.conj().T P = np.zeros(len(fd_grid)) for i, fd in enumerate(fd_grid): a = np.exp(2j * np.pi * fd * TS * np.arange(M)) denom = np.abs(a.conj() @ En @ EnH @ a) + 1e-12 P[i] = 1.0 / denom return P def top_k_peaks(P, fd_grid, K, min_dist=5): peaks, _ = find_peaks(P, distance=min_dist) if len(peaks) == 0: return np.array([]), np.array([]) order = np.argsort(P[peaks])[::-1][:K] return fd_grid[peaks[order]], P[peaks[order]] def estimate_cfo_doppler(X, K, cfo_grid, fd_grid): best_score = -np.inf best_cfo = None best_fd = None t_m = TS * np.arange(M)[:, None] for cfo in cfo_grid: y = X * np.exp(-2j * np.pi * cfo * t_m) P = music_spectrum(y, K, fd_grid) pk_fd, pk_val = top_k_peaks(P, fd_grid, K) if len(pk_val) < K: continue score = float(np.mean(pk_val)) if score > best_score: best_score = score best_cfo = cfo best_fd = pk_fd return best_cfo, best_fd, best_score def main(): YF, targets, cfo_true = generate_echoes(snr_db=15, seed=1) X, r0 = extract_range_doppler_snapshot(YF) fd_grid_coarse = np.arange(-1200, 1200.1, 5) cfo_grid_coarse = np.arange(0, 5000, 100) cfo_c, fd_c, score_c = estimate_cfo_doppler(X, K=2, cfo_grid=cfo_grid_coarse, fd_grid=fd_grid_coarse) print("粗搜索 CFO 估计:", cfo_c, "Hz, score =", score_c) print("粗搜索 多普勒估计:", fd_c, "Hz") cfo_grid_fine = np.arange(max(0, cfo_c - 100), cfo_c + 100.1, 5) fd_min = np.floor(np.min(fd_c)) - 50 fd_max = np.ceil(np.max(fd_c)) + 50 fd_grid_fine = np.arange(fd_min, fd_max + 0.1, 1) cfo_f, fd_f, score_f = estimate_cfo_doppler(X, K=2, cfo_grid=cfo_grid_fine, fd_grid=fd_grid_fine) print("细搜索 CFO 估计:", cfo_f, "Hz, score =", score_f) print("细搜索 多普勒估计:", fd_f, "Hz") print("真实 CFO:", cfo_true, "Hz") print("真实目标多普勒:", [round(t["fd"], 1) for t in targets]) if __name__ == "__main__": main()运行后会看到类似这种输出：
粗搜索 CFO 估计: 2800.0 Hz, score = 5.6e6 粗搜索 多普勒估计: [466.0 -350.0] Hz 细搜索 CFO 估计: 2795.0 Hz, score = 8.9e8 细搜索 多普勒估计: [466.0 -350.0] Hz 真实 CFO: 2800.0 Hz 真实目标多普勒: [466.7, -350.0]如果你想把CFO搜索范围按晶振等级调整，换算公式是 (f_{CFO}=f_c \times ppm / 10^6)。代码里的范围是0~5000 Hz，对应3.5 GHz下约1.43 ppm以内的晶振场景。换成更高载频或者更差晶振时，要把粗搜索范围的步长适当加大，否则计算时间会明显增加。
我个人在实际项目里，把这个联合搜索思路放到粗同步之后的多普勒估计环节，之前“先同步再补偿”的路子直接被替换掉了。少了一个环节，反而多了一个意外好处：不怕同步模块的残差，因为残差本质上也是一种没有被完全估计出来的CFO，联合搜索会把剩下的部分继续吃掉。如果你也正被时钟异步困扰，建议别一上来就写一堆校准逻辑，先把这个思路跑通一版，说不定问题就无声无息地消失了。
