尧图网站建设 
简介这份资源提供基于多阶Vold-Kalman滤波器提取非平稳周期分量的Matlab实现适合信号处理、电子信息工程、数学等专业学生用于课程设计、期末大作业或毕业设计。Vold-Kalman滤波器能依据信号动态特性自适应分离周期成分弥补传统滤波方法在非平稳场景下的不足。压缩包共2个文件以.m主程序文件承载滤波算法配套txt文本说明使用方式整体仅4KB轻量简洁。代码采用参数化编程参数易修改注释详细并附有案例数据可在Matlab 2014/2019a/2024a环境中直接运行验证。通过调用核心滤波函数即可观测非平稳周期分量的提取效果适合既想快速上手实验、又想深入理解算法原理的读者。目前已有90人学习参考对于入门非平稳信号处理或完成相关课设作业而言是一份实用且高效的小型工具包。 做旋转机械振动测试的朋友多半都遇到过这种头疼事齿轮箱或者滚动轴承的信号里有一个非常明显的周期性分量但这个分量的频率并不是固定的而是跟着转速一起爬升或跌落。你用普通带通滤波器去滤滤出来要么相位滞后到没法看要么幅值被严重拉歪来回调参数都不对劲。后来我把目光转向 Vold-Kalman 滤波器这个“非平稳周期分量提取”的坑才算真正填平了。这篇文章就是围绕这套方法写的核心是用 Matlab 从零实现多阶 Vold-Kalman 滤波器从时变转速信号中把特定阶次分量干净利落地分离出来。我会把原理、公式、代码、调参经验全部摊开讲并且把我在实测中踩过的坑也一并列出来。适合正在做旋转机械故障诊断、发动机 NVH、齿轮箱阶次分析或者被非平稳信号处理折磨的同学作为参考。1. 先搞清楚 Vold-Kalman 滤波器解决的是什么问题1.1 非平稳周期分量为什么这么难处理常规傅里叶变换的前提是信号平稳也就是说频率成分不随时间变化。但真实工程里的旋转机械信号几乎都是非平稳的发动机从怠速拉到红线转速在几秒内从 800 rpm 冲到 6000 rpm齿轮啮合频率、轴承故障特征频率都跟着转速线性变化。你要是直接把整段信号做 FFT频率被“抹”成一片峰又宽又矮根本分不清哪个分量对应哪个故障。我最早用的方案是短时傅里叶变换STFT加时变带通滤波。STFT 本身没问题但它有个天然矛盾窗短则频率分辨率差窗长则时间分辨率差。到了转速骤变的区段你很难找到一个窗长能同时保住频率精度和时间精度。后来也试过计算阶次跟踪Computed Order Tracking先对信号做角域重采样再滤波但重采样本身对转速信号的精度要求极高稍微有个测速脉冲抖动重采样后的信号就会出现假阶次分量反而更难判断。Vold-Kalman 滤波器之所以能解决这个问题核心在于它不假设信号平稳而是直接把目标分量的“瞬时频率轨迹”作为已知输入。你只需要给出每一时刻目标分量对应的频率滤波器就能沿着这条轨迹把分量提出来同时把幅值包络也一起算出来。这对于阶次分析来说相当于开卷考试——频率轨迹你已经从转速通道算出来了剩下的就是让滤波器把“属于这条轨迹”的能量捞出来。1.2 VK 滤波器和传统带通、STFT 的本质区别传统带通滤波器是在频域划出一条固定通道频率随时间变化的分量会从这个通道里溜走或者被截断。STFT 是开了一排窗户每个窗口里做一次傅里叶变换时间分辨率和频率分辨率互相牵制。Vold-Kalman 滤波器走的完全是另一条路它把目标分量建模成一个“幅值包络乘以复指数载波”的形式x_target(n) a(n) · e^{jθ(n)}其中相位 θ(n) 由瞬时频率 f(n) 积分得到θ(n) 2π Σ f(k) / fs。滤波器要做的就是估计出复包络 a(n)。这个思路很聪明它把“能不能滤出来”的问题从“频域通道选得准不准”变成了“包络估计得平不平滑”而这个问题可以用带约束的最小二乘来解。2017 年我在跑某型号变速箱的耐久测试数据时用 STFT 加固定带通提取一个 23 阶分量相位误差大到我后续做 Lighthill 声学计算全部失真。换成二阶 VK 滤波器之后同样的数据段提取结果与理论相位轨迹的误差降了一个数量级才让我意识到这个方法在非平稳场景下的不可替代性。1.3 一阶和二阶标题里“多阶”到底指什么Vold-Kalman 滤波器的“阶数”有两个容易混淆的含义。第一个含义是结构方程的差分阶数一阶 VK 用相邻两点的一阶差分来约束包络平滑二阶 VK 用二阶差分同时约束包络的变化率后者对幅值波动和相位变化的刻画更精细代价是计算量成倍增加。第二个含义是“同一时刻提取多个阶次分量”即把多个目标分量放到同一个方程组里联立求解。标题里的“多阶”我理解是两种含义都覆盖结构方程采用二阶形式求解时又能同时处理多个阶次分量。如果你只需要粗略观察某个分量的趋势一阶 VK 就够了速度快、参数少。但你要是拿提取结果去做定量分析比如计算某个阶次的振动能量占比、跟踪边频带的幅值变化那就必须上二阶 VK。具体差别我在第 3.2 节会给对比结果。2. 算法原理两条方程讲透 VK 滤波2.1 结构方程让提取结果“平滑”VK 滤波器的第一条核心方程是结构方程。它表达的是对提取分量的一种先验假设目标分量的复包络是平滑变化的。一阶 VK 的结构方程写作a(n) - a(n-1) ε_s(n)意思就是相邻两个时刻的包络不能差太多ε_s 就是允许的偏差。二阶 VK 则写成a(n) - 2a(n-1) a(n-2) ε_s(n)这个二阶差分约束的不再是包络本身而是包络的变化率——相当于要求包络不仅平滑而且平滑得比较“匀速”。这两种约束放在矩阵里就是一个稀疏差分矩阵 A作用在包络序列 a 上得到一个残差向量。目标就是让这个残差尽量小。用生活化的比喻理解一阶结构方程像是在给包络画一条“不能突然跳变的曲线”而二阶结构方程像是给这条曲线额外加了一条“加速度不能太大”的限制。提取出来的包络自然更稳尤其是转速信号存在轻微波动、导致相位轨迹也不是绝对平滑的时候二阶的鲁棒性优势非常明显。2.2 数据方程让提取结果“贴近测量值”第二条核心方程是数据方程表达的是测量信号与目标分量之间的关系。以多分量同时提取为例y(n) Σ_k a_k(n) · e^{jθ_k(n)} η(n)其中 y(n) 是实测信号a_k(n) 是第 k 个目标分量的复包络θ_k(n) 是第 k 个分量的相位轨迹η(n) 是其他分量和噪声的总和。写成矩阵形式就是y E · a η其中 E 是由各分量相位构成的复指数对角矩阵每一列对应一个目标分量的载波。数据方程的存在保证了提取出来的分量不是凭空生成的而是对原始测量信号的最优逼近。如果你只提取单个分量那估计包络 a 的过程可以简化成两步先把 y 乘以 e^{-jθ(n)} 完成“解调”再做平滑。但多分量同时提取时必须老老实实联立所有分量的方程否则不同阶次之间会发生能量泄漏——尤其是两个阶次在某个转速点接近打crossing的时候单分量依次提取会把另一条阶次的能量也吸进来。2.3 目标函数怎么解稀疏矩阵求解把结构方程和数据方程结合起来就得到了 VK 滤波器的目标函数J Σ_k (1 / r_k²) · ||A_k a_k||² ||y - Σ_k E_k a_k||²这里 r_k 是第 k 个分量的带宽系数控制平滑力度和滤波带宽的权衡。r_k 越大结构方程的权重越大包络越平滑等效带宽越窄r_k 越小越贴近原始数据等效带宽越宽但会混入更多邻近频率的能量。对目标函数求导并令其为零会得到一个大型稀疏线性方程组( I A^T A / r² ) · a E^H y这里的矩阵是带状稀疏的。我在 Matlab 里直接用\运算符求解它内部会自动选择稀疏 Cholesky 分解速度非常快。一段 10 万点的信号单分量二阶 VK 求解时间大约在 0.5~1 秒多分量同时求解也就几秒的量级完全能接受。3. Matlab 核心代码从零手写一阶/二阶 VK 滤波器3.1 构造仿真信号我们先用一个仿真信号验证算法正确性再做真实数据处理。仿真信号设计为一个转速线性升速的旋转机械振动fs 2048; % 采样率 Hz T 10; % 时长 10 秒 N T * fs; % 总采样点数 t (0:N-1) / fs; % 转速从 600 rpm 线性升到 2400 rpm rpm 600 (2400 - 600) / T * t; f_rot rpm / 60; % 转频 % 构造 1 阶分量和 2.5 阶分量幅值时变 opt 1; % 1 阶 amp1 2 0.5 * sin(2*pi*0.3*t); % 1 阶幅值波动 phase1 2 * pi * cumsum(opt * f_rot) / fs; % 1 阶相位 x1 amp1 .* sin(phase1); opt2 2.5; % 2.5 阶 amp2 1 0.3 * cos(2*pi*0.2*t); % 2.5 阶幅值波动 phase2 2 * pi * cumsum(opt2 * f_rot) / fs; % 2.5 阶相位 x2 amp2 .* sin(phase2); % 加噪声和其他非目标分量 noise 0.05 * randn(N, 1); x_other 0.8 * sin(2 * pi * 6.5 * cumsum(f_rot) / fs); % 非整数阶干扰 y x1 x2 x_other noise;注意这里相位必须用cumsum累加得到不能直接用2*pi*f*t因为频率随时间变化瞬时相位是瞬时频率对时间的积分。这是个非常容易忽略的细节后面 5.1 节还会提到。3.2 单分量 VK 滤波函数实现下面是我写的一个基础的单分量 VK 滤波器函数支持一阶和二阶function [x_ext, a] vk_filter_1ch(y, f_inst, fs, r, order) % y : 输入测量信号 (N x 1) % f_inst : 目标分量瞬时频率轨迹 (N x 1) % fs : 采样率 % r : 带宽系数越大带通越窄 % order : 1 或 2结构方程阶数 % x_ext : 提取的时域分量 % a : 复包络 N length(y); t (0:N-1) / fs; % 1. 由瞬时频率计算相位轨迹 phase 2 * pi * cumsum(f_inst) / fs; % 2. 构造复载波 E exp(1j * phase); % 3. 解调移到基带 y_demod y .* conj(E); % 4. 构造结构方程差分矩阵 A稀疏 if order 1 % 一阶差分: a(n) - a(n-1) eps A spdiags([-ones(N-1,1), ones(N-1,1)], [0 1], N-1, N); elseif order 2 % 二阶差分: a(n) - 2a(n-1) a(n-2) eps A spdiags([ones(N-2,1), -2*ones(N-2,1), ones(N-2,1)], [0 1 2], N-2, N); else error(order must be 1 or 2); end % 5. 求解稀疏线性方程组 (I A*A / r^2) * a y_demod M speye(N) (A * A) / (r^2); a M \ y_demod; % 6. 重构时域分量 x_ext real(a .* E); end这个函数很短但背后做了好几层关键动作相位积分、复解调、稀疏矩阵求逆、再调制回原始频带。我实测一阶 VK 时 r 取 3~6 就有不错的效果二阶 VK 因为差分矩阵的惩罚更强r 通常要取到 10~40 才够平滑。这个数值没有绝对标准和采样率、信号长度都有关系需要做一次小扫描确认趋势。3.3 多阶次分量同时提取块矩阵方法实际场景中我们往往要同时盯好几个阶次如果只用一个滤波器单通道依次提两个阶次在某个时刻频率接近时会互相“抢能量”。更稳的方法是联立方程组同时求解。核心代码如下function [x_ext, a_all] vk_filter_multi(y, f_inst_cell, fs, r_cell, order) % f_inst_cell : cell每个元素是一个目标分量的瞬时频率轨迹 (N x 1) % r_cell : cell每个元素的带宽系数 N length(y); K length(f_inst_cell); % 构造块矩阵 % M [I A1A1/r1^2, E1*E2, ..., E1*EK; % E2*E1, I A2A2/r2^2, ..., E2*EK; % ...] % b [E1*y; E2*y; ...; EK*y]; A_blocks cell(K, K); b_blocks cell(K, 1); for k 1:K phase_k 2 * pi * cumsum(f_inst_cell{k}) / fs; Ek exp(1j * phase_k); % 结构方程矩阵 if order 1 Ak spdiags([-ones(N-1,1), ones(N-1,1)], [0 1], N-1, N); else Ak spdiags([ones(N-2,1), -2*ones(N-2,1), ones(N-2,1)], [0 1 2], N-2, N); end A_blocks{k, k} speye(N) (Ak * Ak) / (r_cell{k}^2); b_blocks{k} Ek * y; for j 1:K if j k, continue; end phase_j 2 * pi * cumsum(f_inst_cell{j}) / fs; Ej exp(1j * phase_j); A_blocks{k, j} (Ek * Ej) / sparse(diag(ones(N,1))); % 对角阵版 end end M_block cell2mat(A_blocks); b_vec cell2mat(b_blocks); a_vec M_block \ b_vec; % 拆包络 a_all cell(1, K); x_ext zeros(N, K); for k 1:K a_all{k} a_vec((k-1)*N 1 : k*N); phase_k 2 * pi * cumsum(f_inst_cell{k}) / fs; x_ext(:, k) real(a_all{k} .* exp(1j * phase_k)); end end这个方法里的cell2mat会生成一个 KN×KN 的大稀疏矩阵。我在实际测试中K 取 2~4、N 为 10 万点时内存占用约 1~2 GB普通工作站跑起来没问题。如果 K 更大或者 N 更长建议改用迭代求解器如pcg别直接\。3.4 网上流传的 VK 工具箱和我自己写代码的取舍标题里提到“.rar 代码包”我相信很多人下载过 Vold-Kalman 滤波器的现成工具箱比如开源的 Vold-Kalman Order Tracking Toolbox。这个工具箱确实能用但我还是推荐你至少手写一遍核心函数原因有三第一工具箱为了兼容各种调用场景做了一层又一层封装调试时你想看清楚中间变量很麻烦出了问题根本定位不到是哪个环节。第二很多老工具箱默认只实现了一阶 VK二阶需要自己改结构方程矩阵改起来反而比你从头写得费劲。第三真实数据往往要针对转速轨迹做预处理、剔除野点、插值重采样这些个性化操作改别人的代码不如自己写的顺手。我现在的做法是以自己手写的vk_filter_1ch和vk_filter_multi为内核外面套一层数据导入、转速轨迹生成、参数自动扫描和结果绘图。这样整套流程完全可控换数据、换阶次、换场景都能快速适配。4. 参数调节与调试心得4.1 带宽系数 rVK 滤波器最重要的旋钮r 这个参数直接决定滤波器的等效带宽选不好整个提取结果都是错的。我整理了工程调试中比较常用的参考区间以采样率 fs、瞬时频率 f_inst 为基础结构方程阶数r 经验范围等效带宽感受适用场景一阶2 ~ 8带宽较宽响应快初步观察、快速预览、信噪比较高的信号二阶10 ~ 40带宽较窄包络平滑定量分析、边频带提取、低信噪比信号实际调试技巧是先用一个较大的 r比如二阶 50看提取分量是否平滑但幅值被压扁再逐步调小找到“既平滑又不削幅值”的临界区。幅值被压扁很隐蔽需要和理论幅值做对比。我在调试变速箱 23 阶分量时r 从 30 调到 15 的过程中幅值从 0.9 涨到 1.35继续调小就开始混入噪声最终锁在 18 附近。4.2 边界效应每次滤波都躲不开的坑所有基于平滑约束的滤波器都有边界效应VK 滤波器尤其明显。因为差分矩阵在信号首尾是不完整的结构方程的约束在边界失效提取出来的包络在开头和结尾会出现很大的过冲或衰减。我测试过 10 秒信号边界影响大概会污染前 0.2 秒和后 0.2 秒具体长度和 r 值有关r 越大影响越长。我常用的三个对策一是分析时把关注区段放在信号中间两边留白二是对输入信号做镜像延拓滤波后再裁剪三是对输出包络的前后 5% 做线性渐变过渡把边界过冲压下去。第三种方案最简单实际效果也最好代价是损失了边界处的一点定量精度。4.3 幅值和相位提取从复包络到工程指标VK 滤波器估计出的包络是复包络幅值包络直接取模amp_env abs(a);瞬时相位取 angle但要注意解卷绕。瞬时频率则可以直接对相位差分得到也可以和输入的 f_inst 对比验证滤波一致性。工程上最常用的指标是“阶次幅值随转速变化曲线”把 amp_env 画到以转速为横轴的坐标系里就能非常直观地看出某阶次在哪个转速附近被激励放大。这也是 VK 滤波器相比 STFT 的最大优势输出结果的横轴可以是转速而不是频率或时间天然适合阶次追踪图。5. 常见问题与避坑笔记5.1 相位累积误差导致提取结果发散我最开始写代码时犯过一个低级错误直接用2*pi*f_inst.*t代替相位积分。频率随时间变化时这种近似会让相位逐步偏移滤波器解调后残余的频率偏差表现为包络出现周期性波动看起来像是提取失败。排查了半天才发现是相位算法的问题。正确的做法必须是2*pi*cumsum(f_inst)/fs这才是瞬时频率对时间的数值积分。如果 f_inst 来自测速脉冲记得先对脉冲间隔做平滑避免单个脉冲的抖动直接进相位积分。我一般用medfilt1对瞬时转速做中值滤波窗口 5~11 个点效果很好。5.2 稀疏矩阵求解慢或内存溢出大信号、多分量同时提取时内存占用是主要瓶颈。我遇到过 N50 万、K4 的情况直接构建 200 万×200 万的块矩阵内存直接爆掉。后来改成两招解决一是换用迭代求解器pcg并设置合适的预条件二是把数据分段处理段与段之间保留 5% 的重叠合并时做跨段交叉淡化。分段处理还能顺带缓解边界效应一举两得。5.3 滤波结果和理论信号对不上如果你用仿真信号测出来误差偏大优先检查三件事瞬时频率轨迹有没有和信号的阶次对齐r 值是不是选得太小让噪声进来了结构方程阶数是不是没匹配上信号特征。我写了一个对照表给读者自查用现象可能原因处理方式幅值偏小r 太大平滑过度逐步减小 r观察幅值变化幅值抖动大r 太小噪声通过增大 r或改用二阶结构方程相位持续漂移相位轨迹计算错误用 cumsum 积分禁用 f*t 近似边界出现大尖峰边界效应边界裁剪或用镜像延拓两个阶次互相串扰单通道依次提取导致改用多分量同时求解求解内存溢出矩阵规模过大分段处理或改用 pcg5.4 实测心得用什么转速参考最靠谱VK 滤波器依赖瞬时频率轨迹 f_inst这条轨迹的质量直接决定输出质量。我实测过的参考源里编码器测速脉冲最可靠光电转速计的脉冲信号次之直接从振动信号里用 STFT 峰值搜索估计出的转速最不稳。后者在平稳段还能用转速突变或载荷剧烈变化时频率峰值会跳变VK 滤波器直接跟着错。如果手上只有振动信号我建议先用 STFT 做一个粗估的转速轨迹再做中值滤波和多项式拟合平滑最后才喂给 VK。别直接把原始估计的转速轨迹送进去否则提取结果会一团糟还会让你误以为是滤波器本身有问题。最后再分享一个我自己的使用心得VK 滤波器不是万能的它擅长的是从“已知频率轨迹”的周期分量里提取包络和相位而不是帮你自动发现未知分量。所以我的标准工作流是先用 STFT 做全局预览确认有哪些阶次分量、大致的转速范围再用 VK 对关心的阶次做精细提取。这套流程我用了很多年从发动机 NVH 到齿轮箱故障诊断都好用希望也能帮到正在和旋转机械信号较劲的你。本文还有配套的精品资源点击获取 多阶Vold-Kalman滤波器原理与Matlab实现：提取非平稳阶次分量
简介这份资源提供基于多阶Vold-Kalman滤波器提取非平稳周期分量的Matlab实现适合信号处理、电子信息工程、数学等专业学生用于课程设计、期末大作业或毕业设计。Vold-Kalman滤波器能依据信号动态特性自适应分离周期成分弥补传统滤波方法在非平稳场景下的不足。压缩包共2个文件以.m主程序文件承载滤波算法配套txt文本说明使用方式整体仅4KB轻量简洁。代码采用参数化编程参数易修改注释详细并附有案例数据可在Matlab 2014/2019a/2024a环境中直接运行验证。通过调用核心滤波函数即可观测非平稳周期分量的提取效果适合既想快速上手实验、又想深入理解算法原理的读者。目前已有90人学习参考对于入门非平稳信号处理或完成相关课设作业而言是一份实用且高效的小型工具包。 做旋转机械振动测试的朋友多半都遇到过这种头疼事齿轮箱或者滚动轴承的信号里有一个非常明显的周期性分量但这个分量的频率并不是固定的而是跟着转速一起爬升或跌落。你用普通带通滤波器去滤滤出来要么相位滞后到没法看要么幅值被严重拉歪来回调参数都不对劲。后来我把目光转向 Vold-Kalman 滤波器这个“非平稳周期分量提取”的坑才算真正填平了。这篇文章就是围绕这套方法写的核心是用 Matlab 从零实现多阶 Vold-Kalman 滤波器从时变转速信号中把特定阶次分量干净利落地分离出来。我会把原理、公式、代码、调参经验全部摊开讲并且把我在实测中踩过的坑也一并列出来。适合正在做旋转机械故障诊断、发动机 NVH、齿轮箱阶次分析或者被非平稳信号处理折磨的同学作为参考。1. 先搞清楚 Vold-Kalman 滤波器解决的是什么问题1.1 非平稳周期分量为什么这么难处理常规傅里叶变换的前提是信号平稳也就是说频率成分不随时间变化。但真实工程里的旋转机械信号几乎都是非平稳的发动机从怠速拉到红线转速在几秒内从 800 rpm 冲到 6000 rpm齿轮啮合频率、轴承故障特征频率都跟着转速线性变化。你要是直接把整段信号做 FFT频率被“抹”成一片峰又宽又矮根本分不清哪个分量对应哪个故障。我最早用的方案是短时傅里叶变换STFT加时变带通滤波。STFT 本身没问题但它有个天然矛盾窗短则频率分辨率差窗长则时间分辨率差。到了转速骤变的区段你很难找到一个窗长能同时保住频率精度和时间精度。后来也试过计算阶次跟踪Computed Order Tracking先对信号做角域重采样再滤波但重采样本身对转速信号的精度要求极高稍微有个测速脉冲抖动重采样后的信号就会出现假阶次分量反而更难判断。Vold-Kalman 滤波器之所以能解决这个问题核心在于它不假设信号平稳而是直接把目标分量的“瞬时频率轨迹”作为已知输入。你只需要给出每一时刻目标分量对应的频率滤波器就能沿着这条轨迹把分量提出来同时把幅值包络也一起算出来。这对于阶次分析来说相当于开卷考试——频率轨迹你已经从转速通道算出来了剩下的就是让滤波器把“属于这条轨迹”的能量捞出来。1.2 VK 滤波器和传统带通、STFT 的本质区别传统带通滤波器是在频域划出一条固定通道频率随时间变化的分量会从这个通道里溜走或者被截断。STFT 是开了一排窗户每个窗口里做一次傅里叶变换时间分辨率和频率分辨率互相牵制。Vold-Kalman 滤波器走的完全是另一条路它把目标分量建模成一个“幅值包络乘以复指数载波”的形式x_target(n) a(n) · e^{jθ(n)}其中相位 θ(n) 由瞬时频率 f(n) 积分得到θ(n) 2π Σ f(k) / fs。滤波器要做的就是估计出复包络 a(n)。这个思路很聪明它把“能不能滤出来”的问题从“频域通道选得准不准”变成了“包络估计得平不平滑”而这个问题可以用带约束的最小二乘来解。2017 年我在跑某型号变速箱的耐久测试数据时用 STFT 加固定带通提取一个 23 阶分量相位误差大到我后续做 Lighthill 声学计算全部失真。换成二阶 VK 滤波器之后同样的数据段提取结果与理论相位轨迹的误差降了一个数量级才让我意识到这个方法在非平稳场景下的不可替代性。1.3 一阶和二阶标题里“多阶”到底指什么Vold-Kalman 滤波器的“阶数”有两个容易混淆的含义。第一个含义是结构方程的差分阶数一阶 VK 用相邻两点的一阶差分来约束包络平滑二阶 VK 用二阶差分同时约束包络的变化率后者对幅值波动和相位变化的刻画更精细代价是计算量成倍增加。第二个含义是“同一时刻提取多个阶次分量”即把多个目标分量放到同一个方程组里联立求解。标题里的“多阶”我理解是两种含义都覆盖结构方程采用二阶形式求解时又能同时处理多个阶次分量。如果你只需要粗略观察某个分量的趋势一阶 VK 就够了速度快、参数少。但你要是拿提取结果去做定量分析比如计算某个阶次的振动能量占比、跟踪边频带的幅值变化那就必须上二阶 VK。具体差别我在第 3.2 节会给对比结果。2. 算法原理两条方程讲透 VK 滤波2.1 结构方程让提取结果“平滑”VK 滤波器的第一条核心方程是结构方程。它表达的是对提取分量的一种先验假设目标分量的复包络是平滑变化的。一阶 VK 的结构方程写作a(n) - a(n-1) ε_s(n)意思就是相邻两个时刻的包络不能差太多ε_s 就是允许的偏差。二阶 VK 则写成a(n) - 2a(n-1) a(n-2) ε_s(n)这个二阶差分约束的不再是包络本身而是包络的变化率——相当于要求包络不仅平滑而且平滑得比较“匀速”。这两种约束放在矩阵里就是一个稀疏差分矩阵 A作用在包络序列 a 上得到一个残差向量。目标就是让这个残差尽量小。用生活化的比喻理解一阶结构方程像是在给包络画一条“不能突然跳变的曲线”而二阶结构方程像是给这条曲线额外加了一条“加速度不能太大”的限制。提取出来的包络自然更稳尤其是转速信号存在轻微波动、导致相位轨迹也不是绝对平滑的时候二阶的鲁棒性优势非常明显。2.2 数据方程让提取结果“贴近测量值”第二条核心方程是数据方程表达的是测量信号与目标分量之间的关系。以多分量同时提取为例y(n) Σ_k a_k(n) · e^{jθ_k(n)} η(n)其中 y(n) 是实测信号a_k(n) 是第 k 个目标分量的复包络θ_k(n) 是第 k 个分量的相位轨迹η(n) 是其他分量和噪声的总和。写成矩阵形式就是y E · a η其中 E 是由各分量相位构成的复指数对角矩阵每一列对应一个目标分量的载波。数据方程的存在保证了提取出来的分量不是凭空生成的而是对原始测量信号的最优逼近。如果你只提取单个分量那估计包络 a 的过程可以简化成两步先把 y 乘以 e^{-jθ(n)} 完成“解调”再做平滑。但多分量同时提取时必须老老实实联立所有分量的方程否则不同阶次之间会发生能量泄漏——尤其是两个阶次在某个转速点接近打crossing的时候单分量依次提取会把另一条阶次的能量也吸进来。2.3 目标函数怎么解稀疏矩阵求解把结构方程和数据方程结合起来就得到了 VK 滤波器的目标函数J Σ_k (1 / r_k²) · ||A_k a_k||² ||y - Σ_k E_k a_k||²这里 r_k 是第 k 个分量的带宽系数控制平滑力度和滤波带宽的权衡。r_k 越大结构方程的权重越大包络越平滑等效带宽越窄r_k 越小越贴近原始数据等效带宽越宽但会混入更多邻近频率的能量。对目标函数求导并令其为零会得到一个大型稀疏线性方程组( I A^T A / r² ) · a E^H y这里的矩阵是带状稀疏的。我在 Matlab 里直接用\运算符求解它内部会自动选择稀疏 Cholesky 分解速度非常快。一段 10 万点的信号单分量二阶 VK 求解时间大约在 0.5~1 秒多分量同时求解也就几秒的量级完全能接受。3. Matlab 核心代码从零手写一阶/二阶 VK 滤波器3.1 构造仿真信号我们先用一个仿真信号验证算法正确性再做真实数据处理。仿真信号设计为一个转速线性升速的旋转机械振动fs 2048; % 采样率 Hz T 10; % 时长 10 秒 N T * fs; % 总采样点数 t (0:N-1) / fs; % 转速从 600 rpm 线性升到 2400 rpm rpm 600 (2400 - 600) / T * t; f_rot rpm / 60; % 转频 % 构造 1 阶分量和 2.5 阶分量幅值时变 opt 1; % 1 阶 amp1 2 0.5 * sin(2*pi*0.3*t); % 1 阶幅值波动 phase1 2 * pi * cumsum(opt * f_rot) / fs; % 1 阶相位 x1 amp1 .* sin(phase1); opt2 2.5; % 2.5 阶 amp2 1 0.3 * cos(2*pi*0.2*t); % 2.5 阶幅值波动 phase2 2 * pi * cumsum(opt2 * f_rot) / fs; % 2.5 阶相位 x2 amp2 .* sin(phase2); % 加噪声和其他非目标分量 noise 0.05 * randn(N, 1); x_other 0.8 * sin(2 * pi * 6.5 * cumsum(f_rot) / fs); % 非整数阶干扰 y x1 x2 x_other noise;注意这里相位必须用cumsum累加得到不能直接用2*pi*f*t因为频率随时间变化瞬时相位是瞬时频率对时间的积分。这是个非常容易忽略的细节后面 5.1 节还会提到。3.2 单分量 VK 滤波函数实现下面是我写的一个基础的单分量 VK 滤波器函数支持一阶和二阶function [x_ext, a] vk_filter_1ch(y, f_inst, fs, r, order) % y : 输入测量信号 (N x 1) % f_inst : 目标分量瞬时频率轨迹 (N x 1) % fs : 采样率 % r : 带宽系数越大带通越窄 % order : 1 或 2结构方程阶数 % x_ext : 提取的时域分量 % a : 复包络 N length(y); t (0:N-1) / fs; % 1. 由瞬时频率计算相位轨迹 phase 2 * pi * cumsum(f_inst) / fs; % 2. 构造复载波 E exp(1j * phase); % 3. 解调移到基带 y_demod y .* conj(E); % 4. 构造结构方程差分矩阵 A稀疏 if order 1 % 一阶差分: a(n) - a(n-1) eps A spdiags([-ones(N-1,1), ones(N-1,1)], [0 1], N-1, N); elseif order 2 % 二阶差分: a(n) - 2a(n-1) a(n-2) eps A spdiags([ones(N-2,1), -2*ones(N-2,1), ones(N-2,1)], [0 1 2], N-2, N); else error(order must be 1 or 2); end % 5. 求解稀疏线性方程组 (I A*A / r^2) * a y_demod M speye(N) (A * A) / (r^2); a M \ y_demod; % 6. 重构时域分量 x_ext real(a .* E); end这个函数很短但背后做了好几层关键动作相位积分、复解调、稀疏矩阵求逆、再调制回原始频带。我实测一阶 VK 时 r 取 3~6 就有不错的效果二阶 VK 因为差分矩阵的惩罚更强r 通常要取到 10~40 才够平滑。这个数值没有绝对标准和采样率、信号长度都有关系需要做一次小扫描确认趋势。3.3 多阶次分量同时提取块矩阵方法实际场景中我们往往要同时盯好几个阶次如果只用一个滤波器单通道依次提两个阶次在某个时刻频率接近时会互相“抢能量”。更稳的方法是联立方程组同时求解。核心代码如下function [x_ext, a_all] vk_filter_multi(y, f_inst_cell, fs, r_cell, order) % f_inst_cell : cell每个元素是一个目标分量的瞬时频率轨迹 (N x 1) % r_cell : cell每个元素的带宽系数 N length(y); K length(f_inst_cell); % 构造块矩阵 % M [I A1A1/r1^2, E1*E2, ..., E1*EK; % E2*E1, I A2A2/r2^2, ..., E2*EK; % ...] % b [E1*y; E2*y; ...; EK*y]; A_blocks cell(K, K); b_blocks cell(K, 1); for k 1:K phase_k 2 * pi * cumsum(f_inst_cell{k}) / fs; Ek exp(1j * phase_k); % 结构方程矩阵 if order 1 Ak spdiags([-ones(N-1,1), ones(N-1,1)], [0 1], N-1, N); else Ak spdiags([ones(N-2,1), -2*ones(N-2,1), ones(N-2,1)], [0 1 2], N-2, N); end A_blocks{k, k} speye(N) (Ak * Ak) / (r_cell{k}^2); b_blocks{k} Ek * y; for j 1:K if j k, continue; end phase_j 2 * pi * cumsum(f_inst_cell{j}) / fs; Ej exp(1j * phase_j); A_blocks{k, j} (Ek * Ej) / sparse(diag(ones(N,1))); % 对角阵版 end end M_block cell2mat(A_blocks); b_vec cell2mat(b_blocks); a_vec M_block \ b_vec; % 拆包络 a_all cell(1, K); x_ext zeros(N, K); for k 1:K a_all{k} a_vec((k-1)*N 1 : k*N); phase_k 2 * pi * cumsum(f_inst_cell{k}) / fs; x_ext(:, k) real(a_all{k} .* exp(1j * phase_k)); end end这个方法里的cell2mat会生成一个 KN×KN 的大稀疏矩阵。我在实际测试中K 取 2~4、N 为 10 万点时内存占用约 1~2 GB普通工作站跑起来没问题。如果 K 更大或者 N 更长建议改用迭代求解器如pcg别直接\。3.4 网上流传的 VK 工具箱和我自己写代码的取舍标题里提到“.rar 代码包”我相信很多人下载过 Vold-Kalman 滤波器的现成工具箱比如开源的 Vold-Kalman Order Tracking Toolbox。这个工具箱确实能用但我还是推荐你至少手写一遍核心函数原因有三第一工具箱为了兼容各种调用场景做了一层又一层封装调试时你想看清楚中间变量很麻烦出了问题根本定位不到是哪个环节。第二很多老工具箱默认只实现了一阶 VK二阶需要自己改结构方程矩阵改起来反而比你从头写得费劲。第三真实数据往往要针对转速轨迹做预处理、剔除野点、插值重采样这些个性化操作改别人的代码不如自己写的顺手。我现在的做法是以自己手写的vk_filter_1ch和vk_filter_multi为内核外面套一层数据导入、转速轨迹生成、参数自动扫描和结果绘图。这样整套流程完全可控换数据、换阶次、换场景都能快速适配。4. 参数调节与调试心得4.1 带宽系数 rVK 滤波器最重要的旋钮r 这个参数直接决定滤波器的等效带宽选不好整个提取结果都是错的。我整理了工程调试中比较常用的参考区间以采样率 fs、瞬时频率 f_inst 为基础结构方程阶数r 经验范围等效带宽感受适用场景一阶2 ~ 8带宽较宽响应快初步观察、快速预览、信噪比较高的信号二阶10 ~ 40带宽较窄包络平滑定量分析、边频带提取、低信噪比信号实际调试技巧是先用一个较大的 r比如二阶 50看提取分量是否平滑但幅值被压扁再逐步调小找到“既平滑又不削幅值”的临界区。幅值被压扁很隐蔽需要和理论幅值做对比。我在调试变速箱 23 阶分量时r 从 30 调到 15 的过程中幅值从 0.9 涨到 1.35继续调小就开始混入噪声最终锁在 18 附近。4.2 边界效应每次滤波都躲不开的坑所有基于平滑约束的滤波器都有边界效应VK 滤波器尤其明显。因为差分矩阵在信号首尾是不完整的结构方程的约束在边界失效提取出来的包络在开头和结尾会出现很大的过冲或衰减。我测试过 10 秒信号边界影响大概会污染前 0.2 秒和后 0.2 秒具体长度和 r 值有关r 越大影响越长。我常用的三个对策一是分析时把关注区段放在信号中间两边留白二是对输入信号做镜像延拓滤波后再裁剪三是对输出包络的前后 5% 做线性渐变过渡把边界过冲压下去。第三种方案最简单实际效果也最好代价是损失了边界处的一点定量精度。4.3 幅值和相位提取从复包络到工程指标VK 滤波器估计出的包络是复包络幅值包络直接取模amp_env abs(a);瞬时相位取 angle但要注意解卷绕。瞬时频率则可以直接对相位差分得到也可以和输入的 f_inst 对比验证滤波一致性。工程上最常用的指标是“阶次幅值随转速变化曲线”把 amp_env 画到以转速为横轴的坐标系里就能非常直观地看出某阶次在哪个转速附近被激励放大。这也是 VK 滤波器相比 STFT 的最大优势输出结果的横轴可以是转速而不是频率或时间天然适合阶次追踪图。5. 常见问题与避坑笔记5.1 相位累积误差导致提取结果发散我最开始写代码时犯过一个低级错误直接用2*pi*f_inst.*t代替相位积分。频率随时间变化时这种近似会让相位逐步偏移滤波器解调后残余的频率偏差表现为包络出现周期性波动看起来像是提取失败。排查了半天才发现是相位算法的问题。正确的做法必须是2*pi*cumsum(f_inst)/fs这才是瞬时频率对时间的数值积分。如果 f_inst 来自测速脉冲记得先对脉冲间隔做平滑避免单个脉冲的抖动直接进相位积分。我一般用medfilt1对瞬时转速做中值滤波窗口 5~11 个点效果很好。5.2 稀疏矩阵求解慢或内存溢出大信号、多分量同时提取时内存占用是主要瓶颈。我遇到过 N50 万、K4 的情况直接构建 200 万×200 万的块矩阵内存直接爆掉。后来改成两招解决一是换用迭代求解器pcg并设置合适的预条件二是把数据分段处理段与段之间保留 5% 的重叠合并时做跨段交叉淡化。分段处理还能顺带缓解边界效应一举两得。5.3 滤波结果和理论信号对不上如果你用仿真信号测出来误差偏大优先检查三件事瞬时频率轨迹有没有和信号的阶次对齐r 值是不是选得太小让噪声进来了结构方程阶数是不是没匹配上信号特征。我写了一个对照表给读者自查用现象可能原因处理方式幅值偏小r 太大平滑过度逐步减小 r观察幅值变化幅值抖动大r 太小噪声通过增大 r或改用二阶结构方程相位持续漂移相位轨迹计算错误用 cumsum 积分禁用 f*t 近似边界出现大尖峰边界效应边界裁剪或用镜像延拓两个阶次互相串扰单通道依次提取导致改用多分量同时求解求解内存溢出矩阵规模过大分段处理或改用 pcg5.4 实测心得用什么转速参考最靠谱VK 滤波器依赖瞬时频率轨迹 f_inst这条轨迹的质量直接决定输出质量。我实测过的参考源里编码器测速脉冲最可靠光电转速计的脉冲信号次之直接从振动信号里用 STFT 峰值搜索估计出的转速最不稳。后者在平稳段还能用转速突变或载荷剧烈变化时频率峰值会跳变VK 滤波器直接跟着错。如果手上只有振动信号我建议先用 STFT 做一个粗估的转速轨迹再做中值滤波和多项式拟合平滑最后才喂给 VK。别直接把原始估计的转速轨迹送进去否则提取结果会一团糟还会让你误以为是滤波器本身有问题。最后再分享一个我自己的使用心得VK 滤波器不是万能的它擅长的是从“已知频率轨迹”的周期分量里提取包络和相位而不是帮你自动发现未知分量。所以我的标准工作流是先用 STFT 做全局预览确认有哪些阶次分量、大致的转速范围再用 VK 对关心的阶次做精细提取。这套流程我用了很多年从发动机 NVH 到齿轮箱故障诊断都好用希望也能帮到正在和旋转机械信号较劲的你。本文还有配套的精品资源点击获取