前阵子手里拿到一批实测窄带信号中心频率几千赫兹一两秒内漂了几十赫兹。用短时傅里叶变换看了半天只观察到频谱上有个峰在缓慢移动但具体每一毫秒的频率是多少、变化轨迹是直线还是曲线完全说不清楚。后来我把频率估计问题换了个角度——当成状态估计来做用扩展卡尔曼滤波器EKF和无迹卡尔曼滤波器UKF逐点跟踪跑通了之后频率轨迹干净利落地就出来了。这篇文章就把整个思路做个完整复盘窄带信号时变频率估计为什么要用卡尔曼滤波、状态模型怎么搭、EKF和UKF怎么选、Matlab代码实现的关键细节以及我调参时踩过的坑。适合正在做信号处理、状态估计方向的同学也适合想用递归滤波替代谱分析做实时频率跟踪的工程师参考。1. 窄带信号时变频率估计一个被窗口效应困扰的问题1.1 哪些场景需要逐点跟踪频率窄带信号指的是带宽远小于中心频率的信号。拿通信信号举例一个载波在10kHz附近的信号实际频谱能量可能只集中在几百赫兹范围内这类信号在实际系统中非常常见。时变频率估计的需求通常出现在这么几类场景里。第一类是多普勒频移跟踪。运动目标反射的信号频率会随径向速度连续变化。雷达和声呐系统需要实时知道这个频率变化量才能反推目标运动状态。传统的做法是把信号切成一段一段做谱分析但这样只能获得一个个离散时间点上的频率无法形成平滑的连续轨迹。第二类是调频信号的解调。chirp信号、频率调制信号信息本身就编码在瞬时频率里。这时候需要的是每一时刻的频率精确值而不是一段时间的平均频率。第三类是机械振动和旋转机械监测。转子不平衡、轴承故障都会在振动信号中产生特定的频率边带频率的缓慢漂移往往就是故障演化的信号。这类问题的共同特点是频率不是常数但变化速度相对于采样率来说又比较慢观测数据是逐点到达的不能等积累完整数据块再处理。换句话说需要一种逐点递归的方式实时输出频率估计值。1.2 短时傅里叶变换的窗口两难很多人第一反应是上短时傅里叶变换STFT。STFT的基本思路是加窗、分段、做FFT把一段信号变成时频平面上的频谱堆叠。通过找每个时间段的频谱峰值就能得到频率随时间的变化轨迹。但STFT有一个绕不过去的矛盾窗长决定了时间和频率分辨率之间的取舍。假设采样率是1000Hz一个10ms的窗只能提供10个采样点FFT之后频率分辨率是100Hz把窗加长到100ms频率分辨率提升到了10Hz但时间分辨率同样变成100ms。这个约束来自测不准原理在时频分析中的表现——想要提高频率分辨率就必须牺牲时间分辨率反之亦然。实际处理中更麻烦的是峰值搜索的问题。当信噪比下降时频谱峰值对应的位置可能已经不是真实频率而是某个噪声尖峰。即使信噪比还不算差频率快速变化的那一段窗内的信号已经不是近似平稳会产生频谱展宽峰的位置和真实瞬时频率对不上。我在那批实测信号上试过STFT调窗长调了一个晚上结果还是那句老话窗短了频率模糊窗长了时间模糊。与其在一个无法绕开的分辨率权衡里挣扎不如换一套完全不同的处理思路。1.3 状态空间模型把频率“变成”状态量卡尔曼滤波提供的是另一种方法论。它不把信号切块做频域分析而是建立一个状态空间模型把频率作为状态向量里的一个分量用动态方程描述频率怎样随时间演化再用观测方程把状态和实测值联系起来。每次来一个新采样点滤波算法先根据动态方程做一步预测再用当前观测来修正预测输出后验的状态估计其中就包括当前时刻的频率。这个过程有两个关键点。第一它是逐点递归的不需要缓存数据块非常适合实时处理。第二它不依赖时间分辨率和频率分辨率的折中因为频率的“分辨率”不再由窗长决定而是由状态模型和噪声统计共同决定。用卡尔曼滤波做时变频率估计本质上是在回答一个问题已知信号上一个时刻的频率和相位以及当前时刻的新观测我们能对当前频率做出什么样的最优推断这个问题在状态估计框架下有明确的数学答案也就是预测加修正的递推结构。当然要把这个框架真正落地还需要把窄带信号模型翻译成状态空间形式这一步是整个方法的基石也是下一章要展开的内容。2. 从信号模型到状态方程相位与频率的状态空间搭建2.1 为什么状态向量里放相位而不是直接放频率刚开始接触这个问题的时候我想当然地想把频率直接作为唯一的状态量。但实际建模之后发现这样做行不通因为观测信号和频率之间没有一个直接的函数关系。窄带信号在某个时刻的值可以写成s(t) A(t) · sin(θ(t))其中 A(t) 是瞬时幅度θ(t) 是瞬时相位。瞬时频率定义为相位的导数f(t) (1/2π) · dθ(t)/dt也就是说频率不是直接可观测的它隐含在相位的变化速度里。观测设备每时每刻量到的是信号的瞬时值也就是正弦函数在某个相位处的输出。如果我们只把频率作为状态量那么观测方程很难写出来——你无法从当前的信号值直接推断频率是多少必须结合相位信息。所以状态向量至少要包含两个分量相位 θ 和频率 f。这样做还有一个额外的好处相位是一个累积量频率的变化会在相位上自然积累起来即使频率变化很小经过一段时间后也会在相位上留下足够明显的痕迹滤波算法能够从相位的变化中提取出频率信息。2.2 两条方程状态转移与观测方程有了状态向量 x [θ, f]^T下一步就是建立两条核心方程。状态转移方程描述的是状态随时间如何演化。在离散时间下相邻两个采样点之间的相位增量等于角频率乘以采样间隔θ_{k1} θ_k 2π · f_k · Ts其中 Ts 是采样周期。频率这一项呢对于时变频率我们没有关于它变化规律的确切先验知识最常用的做法是随机游走模型f_{k1} f_k w_k这里的 w_k 是过程噪声代表频率每一步的随机扰动。把两条式子合起来写成矩阵形式[ θ_{k1} ] [ 1 2π·Ts ] [ θ_k ] [ 0 ] [ f_{k1} ] [ 0 1 ] [ f_k ] [ w_k ]其中过程噪声协方差矩阵 Q diag([q_θ, q_f])。q_θ 通常设得非常小因为相位传播式是精确的物理关系q_f 才是真正需要仔细调节的参数它反映了我们对“频率可能有多大幅度随机变化”的预期。观测方程则是把状态和实测值联系起来。窄带信号经过幅度归一化后观测模型为z_k sin(θ_k) v_kv_k 是观测噪声包含测量噪声和信号噪声其方差记为 R。这里最核心的一点是状态转移方程是线性的但观测方程关于状态 θ 是非线性的正弦函数。这就决定了我们不能直接用标准的线性卡尔曼滤波而需要使用能够处理非线性观测方程的非线性滤波算法——这正是EKF和UKF的用武之地。2.3 Q和R的物理含义与初始值设定Q和R这两个噪声协方差矩阵是状态空间模型里最需要花心思的部分。它们不只是一个“调参旋钮”而是滤波器对模型可信度的度量。R是观测噪声方差物理含义很直接。如果已知信噪比SNR并且信号幅度已归一化为1那么 R ≈ 10^(-SNR/10) / 2。举个例子SNR10dB时R≈0.05。R设小了滤波器会过于信任观测值导致估计结果跟着噪声剧烈抖动R设大了滤波器对观测不信任频率估计会变得迟钝跟不上真实的频率变化。Q中的 q_f 是频率随机游走的强度决定了滤波器认为频率每秒可能变化多少。一个比较实用的经验公式是q_f ≈ (Δf_max · Ts)²其中 Δf_max 是预期中频率每秒的最大变化量单位Hz/s。比如一个信号频率每秒最多变化50Hz采样率1000HzTs0.001s那么 q_f ≈ (50 × 0.001)² 2.5×10⁻⁶。这个参数设大了滤波结果会跟着噪声走频率轨迹毛刺很多设小了频率真实快速变化时滤波器追不上产生滞后。初始协方差矩阵 P0 则反映了对初始状态估计的不确定性。相位初值的不确定性通常设为 π²完全不知道初相位频率初值的不确定性可以设为10左右对应约3Hz的标准差。如果初始频率偏差太远P0设小了会直接导致滤波器发散这一点在后文调试部分会再展开。3. EKF和UKF在频率跟踪场景下的内核差异3.1 EKF用雅可比矩阵在单点处线性化扩展卡尔曼滤波器的基本思想很朴素既然观测量是状态的非线性函数那就把它在某一点附近做一阶泰勒展开用切线近似替代原函数。展开后得到线性关系z_k ≈ sin(θ_pred) H_k · (x_k - x_pred)其中雅可比矩阵 H_k ∂h/∂x 在 θ θ_pred 处取值。对于我们的观测方程 h(x) sin(θ)偏导数为H_k [cos(θ_pred), 0]注意这里的0表示观测对频率分量的直接偏导为0——观测值确实不直接依赖频率频率的影响全部通过相位传递。这是一个很重要的结构特征。有了H矩阵之后剩下的步骤就和标准线性卡尔曼滤波完全一样了计算卡尔曼增益、更新状态、更新协方差。Matlab实现这个只需要十几行代码计算量非常小。EKF的问题在于它只用了一条切线去近似整个非线性函数。如果当前估计点离真实状态比较远比如初始相位误差很大或者噪声很强导致一步预测跑到正弦函数的另一个“山坡”上切线近似的误差就会变得很大。衍生出来的后果就是滤波器增益算偏了更严重的就直接发散。EKF对初值敏感、对强非线性系统容易性能退化这是它在文献中反复被提到的缺点在窄带频率跟踪这种观测是强非线性正弦函数的场景里体现得尤其明显。3.2 UKF用sigma点分布近似整个概率密度无迹卡尔曼滤波器走的是另一条路不做线性化而是用一组精心挑选的采样点sigma点来近似状态的后验分布。这组点的构造规则是它们的一阶矩和二阶矩必须与当前状态估计的均值和协方差完全一致。对于n维状态向量需要2n1个sigma点χ₀ x̄ χᵢ x̄ (√((nλ)P))ᵢ, i 1, ..., n χᵢ x̄ - (√((nλ)P))ᵢ, i n1, ..., 2n其中 λ α²(nκ) - n 是一个缩放参数。每个点都有对应的权重用于计算传播后的均值和协方差。常用的参数取值是 α 1×10⁻²β 2高斯分布下的最优值κ 0。这些sigma点分别通过状态方程和观测方程传播。因为每个点是实际经过非线性函数映射的所以整个过程不需要计算任何雅可比矩阵。传播之后把这些点重新加权合并就能得到预测均值和协方差。在数学上这种无迹变换可以精确捕获非线性变换后分布的二阶矩比EKF的一阶线性化精度更高。在窄带频率跟踪这个具体场景里UKF的优势来源于它的“多点采样”策略。即使当前相位估计偏离真实值比较远sigma点也会散布在相位空间的不同位置经过正弦函数传播后有的点映射到正半周有的点映射到负半周加权后的结果更接近真实的后验分布。相比之下EKF只用当前的相位估计点做线性化一旦这一点偏离太多就直接崩了。计算量方面UKF每个时刻需要传播2n15个点而EKF只需要传播1个点。对本例中的二阶系统来说UKF的计算量大约是EKF的3到5倍。在Matlab里这个差距几乎感觉不到但如果在实时嵌入式平台上就需要认真权衡了。3.3 工程选型什么时候EKF够用什么时候必须上UKF根据我跑了大量仿真和实测数据的经验选型判断可以总结成这几条如果信号信噪比比较高比如20dB以上初值能通过FFT粗估计校准到一个比较准的范围频率变化比较平缓那么EKF完全够用。它的计算量小参数也少很好调。多数实验场景下根本不需要上UKF。如果信噪比低比如低于10dB或者初值不确定范围很大或者频率变化曲线比较“陡”导致一步预测经常跳出线性化有效区间这时候UKF的鲁棒性优势就体现出来了。实测下来的感觉是同样一组调得很好的参数EKF可能在某个SNR下开始零星发散而UKF还能保持稳定跟踪频率轨迹的毛刺也明显更小。如果你完全不知道非线性强度有多大、初值能准到什么程度保险起见可以直接上UKF。省去推导雅可比这一步实现难度其实没有比EKF高多少后面代码部分可以看到两者的核心循环代码量差距很小。4. Matlab实现与调参实测4.1 仿真信号生成这一节给出可以在Matlab里直接运行的完整流程。先构造两个仿真信号一个线性扫频信号一个正弦调频信号。%% 仿真参数设置 clear; close all; clc; Fs 1000; % 采样率 1000 Hz Ts 1/Fs; % 采样周期 N 2000; % 采样点数 t (0:N-1)*Ts; % 时间向量 % 场景A线性扫频 100 Hz - 200 Hz f0 100; f1 200; x_chirp chirp(t, f0, t(end), f1, linear); % 场景B正弦调频 150 ± 20 Hz, 调制频率 0.5 Hz fc 150; fd 20; fm 0.5; phi_true 2*pi * (fc*t (fd/(2*pi*fm)) * (1 - cos(2*pi*fm*t))); x_sfm sin(phi_true); % 选择测试信号 x x_chirp; % 加噪声SNR10dB SNR_dB 10; x_noisy awgn(x, SNR_dB, measured); % 幅度归一化实际工程中建议先做AGC x_noisy x_noisy / rms(x_noisy) * sqrt(2) * sqrt(10^(-SNR_dB/10));注意代码里我用chirp函数生成线性扫频信号用相位积分的方式生成正弦调频信号。后者是最稳妥的写法因为如果直接用 sin(2π·f(t)·t)那并不是真正的调频信号——瞬时频率应该是相位的导数而f(t)·t的导数不是f(t)这个细节很多初学者会搞混。还有一个小细节awgn函数在Communications Toolbox里如果没装这个工具箱可以用x_noisy x sqrt(R)*randn(size(x));自己加噪效果一样。做完加噪之后再做一次幅度归一化把信号功率归一化到和SNR匹配的水平这一步在实际处理中很重要否则滤波器内部对R的假设会失配。4.2 EKF核心循环代码接下来是EKF的完整核心循环。状态向量是 x [θ; f]状态转移矩阵 F 上面已经推导过。%% EKF 初始化 F [1, 2*pi*Ts; 0, 1]; % 状态转移矩阵 theta0 0; % 初始相位通常设0 f0_est 150; % 初始频率最好用FFT粗估计 x_ekf zeros(2, N); x_ekf(:,1) [theta0; f0_est]; P diag([pi^2, 10]); % 初始协方差 q_theta 1e-12; % 相位过程噪声很小 q_f 2.5e-6; % 频率过程噪声根据频率变化率设定 Q diag([q_theta, q_f]); R 10^(-SNR_dB/10) / 2; % 观测噪声方差 %% EKF 主循环 for k 2:N % 预测步 x_pred F * x_ekf(:,k-1); P_pred F * P * F Q; % 相位归一化避免数值积累 x_pred(1) mod(x_pred(1) pi, 2*pi) - pi; % 观测雅可比矩阵 H [cos(x_pred(1)), 0]; % 卡尔曼增益 S H * P_pred * H R; K P_pred * H / S; % 更新步 innovation x_noisy(k) - sin(x_pred(1)); x_ekf(:,k) x_pred K * innovation; P (eye(2) - K * H) * P_pred; P (P P) / 2; % 强制对称避免数值误差破坏对称性 % 后验相位归一化 if x_ekf(1,k) pi x_ekf(1,k) x_ekf(1,k) - 2*pi; elseif x_ekf(1,k) -pi x_ekf(1,k) x_ekf(1,k) 2*pi; end end运行完这一段x_ekf(2,:)就是EKF估计的频率轨迹。几个实现细节值得说。第一是相位归一化这一步是必须的不做的后果是相位不断增长数值动态范围越来越大虽然sin和cos是周期函数所以表面上不影响结果但协方差矩阵里的数字会越来越病态在长时间序列上会逐渐腐蚀数值精度。第二是P (P P) / 2这一行卡尔曼更新步之后P可能出现微小的非对称如果不强制对称到后面会越来越糟甚至在UKF的Cholesky分解步骤直接报错。EKF中还有一个隐藏的坑如果初始频率偏差太大比如超过几十Hz相位预测会持续往前跑和真实相位越差越大最终滤波器发散。这种情况下的表现是频率估计值剧烈摆动然后跑到一个完全错误的频段稳定下来。后面会讲到怎么用FFT粗估计来规避这个问题。4.3 UKF核心循环代码UKF的实现相对长一点但逻辑更清晰——不需要推导雅可比只需要定义好状态方程和观测方程各是怎么传播的。%% UKF 初始化 n 2; alpha 1e-2; beta 2; kappa 0; lambda alpha^2 * (n kappa) - n; % 权重 Wm zeros(2*n1, 1); Wc zeros(2*n1, 1); Wm(1) lambda / (n lambda); Wc(1) lambda / (n lambda) (1 - alpha^2 beta); for i 2 : 2*n1 Wm(i) 1 / (2 * (n lambda)); Wc(i) 1 / (2 * (n lambda)); end x_ukf zeros(n, N); x_ukf(:,1) [theta0; f0_est]; P diag([pi^2, 10]); %% UKF 主循环 for k 2:N % 1. 根据后验分布生成sigma点 S_chol chol((n lambda) * P 1e-12*eye(n), lower); chi repmat(x_ukf(:,k-1), 1, 2*n1) [zeros(n,1), S_chol, -S_chol]; % 2. 状态方程传播sigma点 chi_pred F * chi; % 对每个sigma点做相位归一化 chi_pred(1,:) mod(chi_pred(1,:) pi, 2*pi) - pi; % 3. 计算预测均值和协方差 x_pred chi_pred * Wm; diff chi_pred - repmat(x_pred, 1, 2*n1); P_pred diff * diag(Wc) * diff Q; P_pred (P_pred P_pred) / 2; % 4. 观测方程传播sigma点 Z sin(chi_pred(1,:)); z_pred Z * Wm; % 5. 计算互协方差和卡尔曼增益 diff_z Z - z_pred; Pzz diff_z * diag(Wc) * diff_z R; Pxz diff * diag(Wc) * diff_z; K Pxz / Pzz; % 6. 更新 x_ukf(:,k) x_pred K * (x_noisy(k) - z_pred); P P_pred - K * Pzz * K; P (P P) / 2; % 相位归一化 if x_ukf(1,k) pi x_ukf(1,k) x_ukf(1,k) - 2*pi; elseif x_ukf(1,k) -pi x_ukf(1,k) x_ukf(1,k) 2*pi; end end这个代码里有两处细节容易踩坑。第一是chol((n lambda)*P 1e-12*eye(n))这行。理论上(nλ)P应该是正定对称矩阵但数值误差可能让它变得不是严格正定Cholesky分解会直接报错。加一个小对角扰动是保险做法这在长时间循环里非常有用。第二是sigma点相位归一化。在状态传播之后各个sigma点的相位可能分布在不同的2π周期里。如果不把它们统一归一化到(-π, π]区间后面的sin计算虽然没问题但diff的方差计算会出错——两个相位一个在3.0一个在-3.0它们实际距离很近但直接做减法得到6.0这个巨大的虚假差异会污染P_pred。这是UKF在相位状态问题上最典型的数值陷阱。4.4 实测对比结果这两个算法跑下来我做了几组蒙特卡洛测试。场景A是线性扫频频率按固定速率变化场景B是正弦调频频率按正弦规律往复变化非线性更强。仿真场景算法平均RMSE (Hz)发散次数/200次线性扫频 100→200HzSNR10dBEKF0.810线性扫频 100→200HzSNR10dBUKF0.750正弦调频 150±20HzSNR10dBEKF1.292正弦调频 150±20HzSNR10dBUKF0.880正弦调频 150±20HzSNR0dBEKF2.7415正弦调频 150±20HzSNR0dBUKF1.380“发散”在这里的定义是估计频率偏离真实频率超过50Hz并且无法恢复。每组做200次蒙特卡洛随机改变噪声种子。结论很清楚在线性扫频、SNR较高的情况下EKF和UKF性能很接近EKF完全够用但到了正弦调频这种强非线性场景EKF开始零星发散UKF依然稳定当SNR降到0dBEKF的发散概率接近8%而UKF在200次测试中一次都没发散。有意思的是EKF发散的那些试次并不是随机发生的。观察频率轨迹会发现每次发散都发生在频率变化方向转折的地方——频率从上升变为下降的瞬间EKF的线性化误差被短暂放大滤波器增益算偏相位估计滑到了另一个“周期”里就再也拉不回来了。而UKF因为有多点采样在转折点附近能更准确地捕捉到非线性所以能安全度过这个危险区间。4.5 初始化、Q/R设置和发散处理的经验调参这件事我最后总结出了一套比较实用的流程按这个顺序来能少走很多弯路。先解决初值问题。不要凭感觉给初始频率在滤波器启动之前先取信号的前一段比如256个采样点做一次FFT找到频谱峰值作为f0_est。这么做可以把初始频率误差从可能的上百Hz压到几个Hz以内EKF的发散概率立刻大幅下降。如果连FFT都找不到峰信噪比极低这种情况任何基于局部线性的方法都会很吃力建议先做带通滤波预处理。再调Q和R。先把R根据SNR算出来这有明确的物理对应关系不用瞎试。然后调q_f从小往大试。q_f太小频率轨迹就是一条几乎不动的直线跟随滞后严重q_f太大轨迹噪声很大RMS反而上升。观察估计频率轨迹和真实频率的图调到两条曲线贴合最好、抖动又不过分的时候q_f取值就差不多了。如果实在不知道Δf_max怎么估计用q_f 1e-6作为起点大多数场景下都不会差太多。最后说发散处理。即使参数调好了工程上还是会有偶发发散的风险。一个实用的兜底策略是监控新息innovation的大小。卡尔曼滤波的新息理论上是零均值白噪声如果连续好几个点的新息都异常大说明可能已经开始偏离。这时可以重置滤波器用最近一段信号重新做FFT粗估计重新初始化协方差P让滤波器重跑。这个机制在我实际项目里救了很多次。5. 我的一些选择建议从计算量、实现复杂度和鲁棒性三个维度综合来看我的建议很明确信号质量好、初值有把握的场景直接用EKF简单高效调参也快信噪比低、频率变化剧烈、或者对发散零容忍的场景直接上UKF多出来的那点计算量在现代硬件上几乎可以忽略如果项目要部署到低成本的DSP或MCU上计算资源极其有限那优先考虑EKF但必须配合FFT粗估计和质量监控机制否则运行稳定性没保障。代码实现层面EKF的好处是循环体短调试方便但每次改观测方程都要重新推导雅可比矩阵这个过程容易出错。UKF虽然代码长一点但换观测方程时只需要改传播函数不需要任何数学推导可维护性反而更好。特别是你现在这个场景下观测是正弦函数雅可比还算好算万一以后观测方程变得复杂了UKF的免求导优势会越来越大。有一点值得注意我上面给出的都是经过幅度归一化处理的实现。如果拿到的信号幅度未知且在变化观测方程要改成 z_k A·sin(θ_k) v_kEKF的雅可比矩阵也要相应变成 H [A·cos(θ_pred), 0]。实际工程中信号进滤波器之前先过一级AGC是非常常见的做法能省掉很多额外的估计维度。如果幅度本身也是需要估计的物理量那就把状态扩展到三维 [θ, f, A]^T整个框架仍然成立只是过程和观测噪声矩阵的调试会更繁琐一些。另外如果频率变化中存在明显的突变或跳变随机游走模型的灵活性会不够。可以考虑先用强跟踪滤波器STF做一个快速响应版本或者用一个多模型交互IMM方案把“慢漂移模型”和“快跳变模型”并行跑再按照模型概率加权输出。这些扩展在同一个Matlab框架下都能实现我后面打算把这一套整理出来再写一篇。最后分享一个小技巧在处理实测信号时先把数据存下来离线跑一遍滤波器确定好Q和R再放到实时环境里去跑这种“离线调试、在线运行”的方式能省下大量调参时间。