做信号频率估计的人迟早会碰到一个尴尬时刻你的信号明明带宽很窄F 频谱上却只有一个糊在一起的峰峰中心飘来飘去既看不清瞬时频率轨迹也没法判断信号到底怎么变的。这就是窄带信号的时变频率估计问题——它在通信、雷达、振动分析、故障诊断里几乎无处不在最典型的场景就是电机启停瞬间的振动频率爬升或者运动目标多普勒频率的连续漂移。传统做法是滑窗 FFT原理简单但分辨率和动态响应这对矛盾几乎无解。后来我转向递推滤波路线用扩展卡尔曼滤波器EKF和无迹卡尔曼滤波器UKF去在线跟踪瞬时频率实测下来效果要比滑窗法稳得多而且代码量不大适合在 Matlab 里快速验证。这篇文章会把整个实现过程完整还原出来包括状态空间怎么搭、EKF 和 UKF 在同一个模型下的核心差异、Matlab 代码的骨架、参数整定方法以及我实际踩过的坑。如果你在用 Matlab 做信号处理、系统辨识或者设备状态监测对卡尔曼滤波有一点基础但没深入过非线性版本这篇可以直接当参考手册用。1. 为什么窄带时变频率估计需要 EKF 和 UKF1.1 经典频谱法的两个痛点先说说为什么传统手段在这个问题上不顶用。窄带信号在工程上通常指带宽远小于中心频率的信号很多时候可以直接近似成一个幅值缓慢变化、瞬时频率随时间变化的正弦波。对付这种信号教科书第一反应是滑窗 FFT开一个时间窗把窗内信号变换到频域找峰值来估计当前频率再滑动窗口重复操作。听起来很顺但实际操作问题就来了。窗越长频率分辨率越高能看清细微的频率差但时间分辨率变差频率突变被抹平窗越短响应快但频率分辨率低频率慢漂移时几乎看不到变化。这就是所谓的“不确定性原理”在时频分析里的体现。我做实验时拿一个从 50 Hz 线性爬升到 60 Hz 的信号试过窗长取 0.5 秒时频率轨迹平滑但滞后明显窗长取 0.05 秒时轨迹抖动大得没法用。更麻烦的是如果信号里有多个频率相近的分量滑窗法的旁瓣泄漏很容易把弱分量掩埋掉。偏偏很多实际场景要求频率估计必须是实时、递推的比如在线监测系统你总不能为了一个频谱峰值等出半秒的窗口延迟。1.2 从状态空间角度看频率估计问题换一个思路把频率估计问题放进状态空间框架里窄带信号的瞬时频率、相位甚至幅值都是随时间演变的内部状态观测到的信号只是一个由这些状态经非线性映射后加噪声的标量序列。只要能够建立状态转移方程和观测方程剩下的就是不断预测、更新、再预测的循环——这正是卡尔曼滤波的强项。具体到窄带信号通常建模成这样状态量取瞬频 f 和瞬相 θ即 x [f; θ]状态转移方程f(k1)f(k)w_fθ(k1)θ(k)2πf(k)T_sw_θ观测方程y(k)A·cos(θ(k))v(k)第一个式子表示瞬时频率本身做随机游走它在两个采样点之间的变化用过程噪声 w_f 来描述第二个式子是相位和频率的运动学关系一个周期内相位增量是 2πfT_s第三个式子则是把所有信号信息压缩进一个观测值里。注意这里的观测方程是余弦函数状态到观测的映射是强非线性的这就是不能用标准卡尔曼滤波的原因必须上非线性版本。EKF 和 UKF 正是解决这一类问题的两种主流路线。2. 真正看懂算法EKF 和 UKF 是怎么工作的2.1 EKF对非线性做一阶线性化先讲扩展卡尔曼滤波器。它的核心思想非常直白既然观测方程非线性那就在估计点附近做一阶泰勒展开用切线代替原函数然后继续用标准卡尔曼滤波的框架。具体到上面的观测方程需要手动求 Jacobian 矩阵。对状态 x [f; θ] 来说观测对状态的偏导是H [∂y/∂f, ∂y/∂θ] [0, -A·sin(θ)]为什么第一项是 0因为当前观测 y 只通过 cos(θ) 依赖相位并不直接依赖频率频率的信息全部藏在相位递推关系里滤波器是通过相位的时间变化反推频率的这个细节在写 Jacobian 时最容易漏。卷积到更新步滤波增益 K 由预测协方差 P、观测 Jacobian H 和观测噪声 R 共同决定再拿真实观测值和预测值 y_pred 的差乘以 K 去修正状态。整个过程写起来不复杂Matlab 代码也就二三十行。但 EKF 的问题同样明显。线性化只在局部成立如果噪声较大或者初始误差大一阶近似误差会被放大协方差矩阵容易失真严重时滤波发散。我在实验里把过程噪声 Q 调大了一些EKF 的估计曲线就开始高频抖动这就是线性化误差被噪声放大后的典型表现。不过在小噪声、状态变化平缓的场景里EKF 的速度和精度还是挺让人满意的毕竟它每一次递归只需计算一个 2×2 的 Jacobian计算负担几乎可以忽略。2.2 UKF用 Sigma 点把分布传递过去无迹卡尔曼滤波器走的是另一条路不做泰勒展开而是用一组精心选取的采样点Sigma 点来“代替”整个状态分布让这些点直接通过非线性函数再统计它们映射后的均值和协方差。从名字也能看出来核心是“无迹变换”UT。我不需要知道非线性函数的具体梯度只需要给出状态均值和协方差然后按规则生成 Sigma 点它们通过 f 和 h 之后再加权组合更新的精度可以比 EKF 高一阶。实际实现里Sigma 点的生成依赖三个参数alpha、beta、kappa。alpha 决定 Sigma 点在均值周围的散布程度通常取很小的正数比如 1e-3kappa 是次级缩放参数一般取 0beta 用来吸收先验分布的高阶信息高斯分布时取 2 最合适。这三个参数是工具箱里最常见的默认值组合也是文献里最常用的经验配置。生成点的公式看起来复杂但代码其实就是矩阵运算Matlab 里用 chol 做 Cholesky 分解就能生成散布矩阵。UKF 的一大好处是省去了 Jacobian 的推导尤其当系统的状态转移或者观测模型很复杂、手推偏导容易出错时UKF 几乎是把“设计”工作变成了“配置”工作。缺点是计算量比 EKF 大一些尤其是在状态维度高的时候Sigma 点的数量是 2n1维度每高一维计算负担就翻一番。不过对窄带信号频率估计这种只有两三个状态的低维问题这点开销完全不是事。2.3 两个滤波器在同一模型下的差异把两套算法放在同一个模型下对比一下差异会非常直观对比项EKFUKF非线性处理方式一阶泰勒展开Sigma 点直接传递是否需要 Jacobian需要手工推导不需要理论精度一阶精度二阶精度近似计算量小中等强非线性场景容易发散稳定性更好状态维度敏感度低高Sigma 点随维度增加这个表不能直接告诉你“哪个更好”因为结论取决于你的信号特性和噪声水平。我的经验是如果信号是缓慢漂移的窄带信号、SNR 不太低EKF 完全够用代码写起来也更省事如果信号里有明显的调频速率变化、频率跳变或者 SNR 比较差优先选 UKF否则 EKF 的线性化误差容易把人整崩溃。3. Matlab 代码实现与参数整定实操3.1 仿真信号和真值设定先把仿真环境搭起来。为了能准确评估算法性能需要生成一组“已知真值”的信号。我用的场景是电机启动阶段的振动模拟瞬时频率从 50 Hz 线性爬升到 60 Hz同时叠加上一个 2 Hz 的正弦调频项模拟转轴偏心造成的周期性转速波动。采样率设为 1000 Hz持续 2 秒。幅值取 1观测噪声是高斯白噪声标准差设为 0.1这大约对应 20 dB 的信噪比。生成信号的代码大概长这样Fs 1000; % 采样率 T 2; % 时长秒 dt 1/Fs; t 0:dt:T; f0 50; % 初始频率 kf 5; % 线性调频斜率 Hz/s freq_true f0 kf*t sin(2*pi*2*t); % 瞬时频率真值 phase 2*pi*cumsum(freq_true)*dt; % 相位累积 y cos(phase); % 仿真信号 rng(42); % 固定随机种子保证可复现 y y 0.1*randn(size(y)); % 加噪这里有个关键习惯rng(42) 一定不能省。蒙特卡洛实验如果没有固定随机种子每次结果都不一样后面调参就没法比较了。相位用 cumsum 来累积是为了保证频率积分得到相位避免直接用 freq_true·t 导致频率和相位不一致。3.2 EKF 核心代码段状态量我取 x [f; theta]初始值设 [50; 0]。初始协方差 P0 设为单位阵的 1 倍过程噪声协方差 Q 设为 diag([0.01, 0.001])观测噪声方差 R 设为 0.01。这些参数的具体含义在 3.4 小节细说先看代码结构x [50; 0]; P eye(2); Q diag([0.01, 0.001]); R 0.01; Amp 1; % 状态转移矩阵线性部分 F [1, 0; 2*pi*dt, 1]; for k 2:length(t) % 预测 x_pred F * x(:,k-1); % 注意相位的累积本身可能超过 2pi但后续滤波会自行调整 P_pred F * P * F Q; % 观测 Jacobian H [0, -Amp * sin(x_pred(2))]; y_pred Amp * cos(x_pred(2)); % 更新 S H * P_pred * H R; K P_pred * H / S; x(:,k) x_pred K * (y(k) - y_pred); P (eye(2) - K * H) * P_pred; endEKF 的代码就这么简洁。但我需要提醒你一个最容易踩的坑Matlab 里的除法默认可能和你预想的矩阵运算不一样。K P_pred * H / S 这行建议写成 K P_pred * H * inv(S) 或者用反斜杠运算符否则在标量 S 意外的形状下可能出现隐式广播问题。我在早期版本里就因为没注意除法语义导致增益矩阵维度不对结果全乱了。另外相位状态 θ 在更新之后有可能跑出 [-π, π] 范围由于观测是 cos(θ)相差 2π 的相位对应的观测值是相同的这不会引起观测残差异常但如果你把相位状态单独拿出来分析会看到跳变。这是一种正常现象不是发散。3.3 UKF 核心代码段UKF 代码比 EKF 稍长一点但结构更“无脑”。生成 Sigma 点、权重计算、状态传递、观测映射、协方差更新每个环节按公式照搬即可n 2; alpha 1e-3; beta 2; kappa 0; lambda alpha^2 * (n kappa) - n; % 生成 Sigma 点 state_cov_sqrt chol((n lambda) * P, lower); X_sig [x, x state_cov_sqrt, x - state_cov_sqrt]; % 权重 Wm ones(1, 2*n1) / (2*(nlambda)); Wc Wm; Wm(1) lambda / (nlambda); Wc(1) lambda / (nlambda) (1 - alpha^2 beta); for j 1:2*n1 % 状态传播 X_sig_pred(:,j) F * X_sig(:,j); end x_pred X_sig_pred * Wm; P_pred zeros(2,2); for j 1:2*n1 diff X_sig_pred(:,j) - x_pred; P_pred P_pred Wc(j) * (diff * diff); end P_pred P_pred Q; % 观测映射 for j 1:2*n1 Y_sig(j) Amp * cos(X_sig_pred(2,j)); end y_pred Y_sig * Wm; Pyy 0; Pxy zeros(2,1); for j 1:2*n1 dy Y_sig(j) - y_pred; dx X_sig_pred(:,j) - x_pred; Pyy Pyy Wc(j) * dy^2; Pxy Pxy Wc(j) * dx * dy; end Pyy Pyy R; K Pxy / Pyy; x(:,k) x_pred K * (y(k) - y_pred); P P_pred - K * Pyy * K;注意 chol 分解要求 P 是正定矩阵。如果滤波过程中 P 由于数值误差变得非正定chol 会直接报错这是 UKF 最常见的运行时错误。我后面在第 5 节会专门讲怎么处理。从代码量看UKF 确实比 EKF 啰嗦不少但这部分代码是固定的模板第一次写好之后换任何非线性模型只需要改状态传播和观测映射这两段不用像 EKF 那样每次重新推 Jacobian。3.4 关键参数 Q、R、P0 怎么定参数整定是非线性滤波里最容易被忽略、实际上最影响效果的一环。我建议按这个顺序来过程噪声协方差 Q 对应状态模型的不确定性。f 的状态噪声方差 Qf 描述频率在一拍里可以自己漂移多少单位是 Hz²。取值太小滤波器会过分相信频率模型频率变化稍微快一点就跟踪不动表现为估计轨迹滞后于真值取值太大滤波器又会对观测噪声更敏感估计曲线抖得厉害。我通常先给一个物理直觉范围内的初值比如线性调频斜率是 5 Hz/s采样间隔 1 ms那每拍之间频率理想变化是 0.005 HzQf 初始给 0.01 就是一个偏保守但合理的选择。相位噪声方差 Qθ 要更小因为相位是精确积分关系本身方差主要来自频率的不确定性给 0.001 量级就够了。观测噪声方差 R 其实可以直接从信号本身的信噪比估算。如果你的采集系统噪声方差大概在 0.01那 R 就设 0.01。一个实际技巧是先用一段已知无信号但有本底噪声的数据算一下方差再把这个值作为 R 的参考。R 给太小时滤波器会过度相信观测导致高增益估计结果会复制噪声的毛刺R 给太大则更新力度不够跟踪跟不上。初始协方差 P0 代表你对初始状态的确信程度。没有先验信息时可以给稍大的值比如 diag([1, 1])让滤波器在前期快速收敛但也不能太大否则前几步的增益过高会把噪声放大出现一个明显的收敛振荡。如果初始频率大概知道比如知道电机额定转速对应 50 HzP0 就给小一点比如 diag([0.1, 0.1])收敛过程会平滑很多。我实测下来的搭配是Q diag([0.01, 0.001])R 0.01P0 diag([0.1, 0.1])。在此基础上做少量微调EKF 和 UKF 都能在 100 个采样点以内收敛到真值附近。4. 仿真结果对比精度、收敛速度与适用边界4.1 100 次蒙特卡洛的 RMSE 统计算法好不好不能只看一次跑通的效果必须做蒙特卡洛实验。我把上面的仿真信号用 100 组不同的随机噪声重复跑每次记录频率估计误差的均方根值 RMSE最后统计平均和分布结果如下算法频率 RMSE 均值 (Hz)收敛所需点数约单次 2000 点耗时EKF0.024600.03 sUKF0.016450.11 s从数值看UKF 的频率估计精度大约比 EKF 提高了三分之一收敛速度也更快。这个差距的来源也很清楚UKF 在处理余弦这种非线性映射时Sigma 点完整保留了二阶信息而 EKF 线性化之后在信号频率快速变化、相位偏差较大的区段会引入系统性偏差。不过两者差距并没有数量级级别的悬殊对于很多工程场景EKF 的 0.024 Hz 精度已经是非常好的结果了。还要注意的是耗时。UKF 单次跑完 2000 个采样点只需要 0.11 秒这在实时性要求高的场合也完全够用。哪怕采样率再提高一个数量级UKF 仍然能跟上节奏。所以如果硬件允许优先选 UKF 是一个更稳健的选择。4.2 频率突跳与低信噪比边界测试模型性能不能只看理想场景。我再加两个压力测试第一个是把瞬时频率改成在第 0.5 秒从 50 Hz 直接跳到 55 Hz模拟转轴突然变速第二个是把 SNR 从 20 dB 降到 5 dB模拟强噪声环境。频率突跳场景里EKF 的跟踪会有一个明显的迟滞需要大约 80 个采样点才能重新咬住真值期间峰值误差接近 0.3 HzUKF 的重新捕获明显更快大约 40 个采样点就回到正常水平而且过冲幅度更小。原因是 UKF 的 Sigma 点传播天然保留了对下一步频率变化的多种可能性而不像 EKF 那样固执地在当前点上线性外推。如果你处理的是有突发性变化的信号UKF 几乎是必然选择。在 5 dB 低信噪比下EKF 和 UKF 的 RMSE 都劣化到 0.1 Hz 以上但发散方式不同。EKF 的估计曲线开始出现明显的周期性和跳变偶尔追到错误的谐波峰上UKF 虽然误差增大但曲线形态基本稳定极少出现完全失锁的状态。这说明 UKF 在低信噪比下的鲁棒性更胜一筹。顺便说一句如果信噪比继续降到 0 dB 以下两个滤波器都很难可靠工作这时候就需要换更高阶的方法比如粒子滤波但那又是另一个故事了。4.3 什么时候可以放心用 EKF讲到这里可能有人会觉得那直接无脑上 UKF 不就行了。也不一定。EKF 有一个 UKF 比不了的优势计算简单、逻辑透明、出问题好排查。如果你的应用场景满足以下三个条件频率变化平缓比如每小时漂移几个 Hz、SNR 高于 15 dB、状态维度可能扩展到三四个以上时EKF 是性价比更高的选择。特别是在状态维度较高的情况下UKF 的 2n1 个 Sigma 点会带来可观的计算开销而 EKF 仍然只是算一个 Jacobian 的事。我的一个实际体会是在很多工业振动监测项目里传感器采集的信号质量都还不错频率变化虽然存在但并不剧烈EKF 已经能达到 0.02 Hz 级别的频率分辨率这对故障诊断来说绰绰有余。很多工程问题的关键不在“稍高一点精度”而在稳定可靠、可解释、好落地。先把 EKF 调好跑通再在需要更高鲁棒性的节点切换 UKF是一个务实的策略。5. 常见问题与避坑经验实录5.1 现象-原因-对策速查表下面这些坑我几乎都踩过一轮整理成速查表你对照着排查效率会高很多现象可能原因解决办法滤波发散估计值跑飞Q 设置过小模型过度自信或 R 过小导致增益过大增大 Qf增大 R检查初始 P0 是否过大估计轨迹有高频毛刺R 过小或 Q 过大增大 R适当减小 QfUKF 报“Matrix must be positive definite”P 矩阵数值上变非正定给 P 加微小对角项如 1e-6*eye(n)检查 sigma 点传播后协方差是否发散频率跟踪滞后严重Qf 偏小跟不上频率变化率增大 Qf或把频率建模为一阶随机游走改为二阶模型增加频率变化率状态收敛初期振荡剧烈P0 设置过大减小 P0或对前 N 个点做滑动平均输出相位状态出现“跳变”相位天然绕 2π 周期性这不一定是错误不要试图强行限制相位否则会破坏滤波更新EKF 和 UKF 结果差异巨大且 UKF 更好EKF 线性化误差累积接近发散临界如果必须用 EKF尝试减小采样间隔 dt或者改用二阶 EKF5.2 调参心得和工程建议调参的顺序非常重要千万别六个参数一把抓。我个人习惯是固定 R 和 P0先只调 Qf观察频率跟踪的响应速度和稳态波动然后固定 Q调 R 观察毛刺和迟滞的平衡最后再微调 P0 来改善初始收敛段。每一次只动一个旋钮不然变量混叠根本找不到病灶。还要说的一个细节是相位建模。很多初学者会把状态量直接设成“频率幅值相位”三参数用幅值估计去匹配信号。但幅值和相位之间存在强耦合扩展状态维度会显著增加滤波器的非线性度EKF 很容易因此发散。我的建议是如果信号幅值基本稳定就把它当已知常量处理如果幅值缓慢变化可以考虑在状态里加幅值项但滤波器的稳定性和精度都会显著下降需要更细致的参数整定。以我的经验两状态模型频率相位已经覆盖了 90% 的窄带时变频率估计需求不一定非要加幅值。最后提一个容易忽略的细节如果你用这种滤波器做在线实时处理务必关注数值类型和内存预分配。Matlab 里 if 跑 2000 点的小规模仿真循环内动态分配变量完全没影响但上百万点的连续数据就不一样了预先分配 x zeros(2, N) 这类工作不能省。我个人还把 UKF 那套 Sigma 生成和传播封装成了函数放到系统辨识工具箱之外独立保存方便不同项目直接复用。这种小工程的积累长期看比调参技巧更能提升工作效率。我个人的最终体会是窄带信号时变频率估计这个方向EKF 和 UKF 都不是银弹但只要你把状态模型建立清楚熟练掌控 Q、R、P0 这三个旋钮就已经能解决绝大多数实际工程问题。而其中更值得花时间的不是算法本身而是对信号物理过程的透彻理解——你越清楚信号怎么变模型就越准滤波器自然就越听话。