齿轮振动细化谱分析：Zoom-FFT原理与Matlab实现

做旋转机械故障诊断的同行应该都有过这种体验振动数据采回来在Matlab里直接做一次FFT整频谱一拉大峰值、小峰值、边频带全叠在一起眼睛盯得发酸也不敢确定某个峰值到底是不是故障特征。尤其在航空发动机齿轮箱这种多级齿轮啮合的场合齿轮啮合频率、各级边频带、轴承特征频率挤在同一个窄频段里常规FFT看全局够用想看细节就抓瞎。这次就来聊聊航空发动机齿轮振动信号细化谱分析把Zoom-FFT的原理、Matlab实现和实测经验一次讲透。文章适合正在写故障诊断程序的研究生、做状态监测的现场工程师也适合想搞懂细化谱和普通FFT差别的信号处理入门者。1. 为什么齿轮振动分析绕不开细化谱1.1 齿轮振动信号到底长什么样齿轮箱振动信号和轴承信号有个明显区别轴承故障特征频率通常出现在较低频段而齿轮啮合振动天然是“高频载波低频调制”的结构。一对齿轮啮合时轮齿交替进入和退出啮合刚度周期性变化产生以啮合频率为基波的振动啮合频率 fm Z × fr其中Z是齿轮齿数fr是转频。比如主动齿轮32齿、转频250Hz啮合频率就是8000Hz。这还没完齿轮副的制造误差、轴不对中、局部剥落都会让啮合振动的幅值发生周期变化相当于以转频fr为调制频率去调制啮合频率于是在频谱上以fm为中心、每隔fr出现一根边频带fm±fr、fm±2fr、fm±3fr……这就是齿轮故障诊断里最值钱的“边频带”。问题在于这些边频带幅度往往比中心啮合频率低一到两个数量级而且和高阶啮合谐波、轴承特征频率混在一起。想从频谱里把齿轮故障特征捞出来本质上就是解决两个问题频率分辨率够不够弱边频能不能从旁瓣和底噪里露出来。1.2 直接做FFT的坑分辨率不够、泄漏淹没普通FFT的频率分辨率很简单Δf fs / N也就是采样率除以FFT点数。现场采集的振动数据往往不是想要多长就有多长转速稳定段可能只有零点几秒或者为了控制数据量只存了有限长度。假设采样率51200Hz、分析窗口0.16秒直接用8192点FFT分辨率是6.25Hz看起来不差但把频谱图拉到8000Hz附近你会发现边频带和啮合频率之间的关系依然模糊。真正头疼的是频谱泄漏。齿轮振动信号里往往存在大幅值的低频分量、轴承冲击分量、其他轴系的啮合分量它们经过Hanning窗等有限长截断后旁瓣并不会收敛到零而是铺在整条频率轴上。弱边频的幅值一旦低于这些旁瓣底就直接被淹没。加大FFT点数确实能改善分辨率但把窗口拉长会引入转速波动问题而且带外强分量的旁瓣并不会因为点数增加就彻底消失。这个时候就需要细化谱分析。1.3 细化谱的思路把镜头拉近到窄带细化谱分析也叫Zoom-FFT思路其实特别朴素。普通FFT像拿广角镜头拍一张全景照整个0到fs/2频段一锅端细化谱则是先确定一个你关心的中心频率fc然后用“复调制低通滤波重采样”三步把fc附近一个很窄的频段单独取出来重新做频谱分析。打个不严谨但好记的比方普通FFT是站在广场中央拍全景细化谱是把长焦镜头架到三脚架上对准广场对面某一栋楼的某一扇窗户。窗户上的一根裂纹只有换长焦才能看见数码变焦是救不了的。细化谱就是振动分析里的“长焦镜头”。这样做的好处有三点一是窄带分析时可用更小的频率间隔把目标频段内的细节铺开二是带外强分量被低通滤波器滤掉不再通过窗函数旁瓣干扰目标频段三是抽取后数据量减少同样的FFT计算量可以覆盖更窄的频率范围适合做逐段精细扫描。2. 细化谱原理拆解Zoom-FFT怎么把镜头拉近2.1 第一步复调制把目标频段搬到零频假设原始信号的频谱以0Hz为中心正负频率各有一半。细化谱第一步是把fc处的频率成分搬到零频附近这是通过乘以一个复指数实现的x_shift(n) x(n) · e^(-j·2π·fc/fs·n)这里有一个关键细节为什么用复指数而不是用cos(2π·fc/fs·n)做实调制因为乘复指数做的是单边带搬移只把正频率方向上的fc附近分量往0频移动不会产生两个互相重叠的搬移副本。乘余弦会把频谱同时往正负方向各搬一次两部分叠在一起就乱套了。实际编程里这步就是一行复数乘法。唯一要注意的是数据量大时这步会产生一个和原数据等长的复数组内存占用翻倍老机器上做长数据得很注意。另一个容易忽略的点是fc不能超过fs/2也不能太靠近0否则搬移后的目标频段会不完整。2.2 第二步低通滤波加抽取把采样率降下来搬移完成后目标频段虽然到了零频附近但数据里还混着远处频段的成分这时候不能直接抽点必须先低通滤波。道理和采样定理一样后续要对信号做整数M倍抽取等效采样率变成fs/M新采样率能无混叠表示的最高频率只有fs/(2M)。如果不把高于这个频率的分量滤干净抽取后它们会折叠回目标频段形成虚假谱峰。滤波器截止频率的设置是个经验活。一般取抽取后新奈奎斯特频率的0.8倍左右也就是截止频率约为0.4×fs/M留出20%的过渡带余量。我用的是8阶巴特沃斯滤波器配合filtfilt做零相位滤波这样不会给谱线引入相位失真。用filtfilt的前提是数据长度要够至少是滤波器阶数的几倍以上否则首尾效果会很差。滤波完成后按M倍抽取每隔M个点取一个样本数据长度变为原来的1/M。抽完之后的信号采样率是fs/M目标频段完整保留频段外面的东西基本被滤掉了。这一步是把“镜头”对准窗户的关键也是细化谱计算量低于直接长数据FFT的原因。2.3 第三步重采样FFT与频率轴换算抽取后的数据做FFT得到的是零频附近的细化谱频率间隔为Δf_zoom (fs/M) / Nfft其中Nfft是细化后做FFT的点数。输出频率轴不能直接当作0Hz附近必须加回中心频率f fc (0 : Nfft/2) × Δf_zoom也就是说细化谱图横轴是fc左侧和右侧的窄带范围。很多新手写出来谱图位置不对十有八九是这步忘了加fc或者把原始采样率当成了抽取后的采样率去算频率间隔。幅值恢复也要注意。细化谱是单边幅值谱FFT结果要乘以2恢复实信号在该频率处的幅值同时用窗函数做幅值修正幅度修正系数是N/sum(win)而不是简单的2/N。如果不加窗才能直接用2/N一旦加了汉宁窗就必须除以窗函数系数之和否则幅值会系统性偏小。2.4 关于分辨率的一个常见误解有朋友问细化谱是不是能突破物理极限让相邻5Hz的两个频率分量也能分开别期望太高。理论上如果你有L点数据原始采样率fs普通全带FFT的频率间隔是fs/L。细化谱把数据每M点抽一个样本采样率变成fs/M有效点数变成L/M算下来频率间隔还是fs/L。也就是说独立频率信息的数量没有增加细化谱并没有“无中生有”。那为什么实战中细化谱就是比直接FFT清楚原因有两个。第一普通FFT如果只用1024点或2048点去覆盖全带分辨率确实有限细化谱把全部数据的分析能力集中到一个窄带里显示间隔小得多。第二也是更关键的带外强分量在细化前已经被低通滤波器干掉了窗函数旁瓣不再掩盖弱边频信噪比大幅提升。说得直白点细化谱提升的是“工程分辨率”和“谱图可读性”不是物理分辨率。真想突破物理分辨率得靠AR模型、MUSIC这类参数化谱估计方法那是另一套玩法。3. Matlab实现可直接套用的细化谱函数与参数调试3.1 核心函数zoom_fft的完整代码下面这个函数我把步骤拆得很开每一行都能对照上面的原理去理解。依赖的是Matlab信号处理工具箱里的butter、filtfilt、hanning函数没有这些函数的话用fir1设计FIR滤波器代替也能实现幅值修正的逻辑不变。function [f, Xz] zoom_fft(x, fs, fc, M) % 细化谱分析Zoom-FFT % 输入 % x - 时域振动信号行或列向量 % fs - 原始采样率Hz % fc - 细化中心频率Hz % M - 细化倍数通常取4~32 % 输出 % f - 细化谱频率轴Hz % Xz - 细化谱幅值单边谱 x x(:).; % 统一转成行向量 L length(x); % 1) 复调制频移把fc附近的频谱搬到零频附近 n 0:L-1; x_shift x .* exp(-1j*2*pi*fc/fs*n); % 2) 低通滤波截止频率取抽取后奈奎斯特频率的0.8倍 % 对应截止频率 fc_cut 0.8 * (fs/M) / 2 0.4*fs/M % 归一化角频率 0.8/M [b, a] butter(8, 0.8/M, low); x_filt filtfilt(b, a, x_shift); % 3) 整数抽取每隔M个点取一个样本 idx 1:M:L; x_dec x_filt(idx); N length(x_dec); % 4) 加汉宁窗补零到2的幂次后做FFT % 注意补零只做插值平滑不提高物理分辨率 Nfft 2^nextpow2(N); win hanning(N); x_win x_dec .* win(:); X fft(x_win, Nfft); % 5) 单边幅值谱按窗函数能量修正幅值 half floor(Nfft/2) 1; Xz 2 * abs(X(1:half)) / sum(win); % 6) 频率轴抽取后采样率变为fs/M df (fs / M) / Nfft; f fc (0:half-1) * df; end这个函数返回的f是从fc开始的单边细化频带plot(f, Xz)就能直接画图。想看fc左侧的谱画图时用xlim设置范围就行因为f覆盖从fc到fc带宽xlim可以往左拉比如xlim([fc-800, fc800])谱图会完整显示中心频率两侧。3.2 用仿真信号验证程序空跑函数不放心我写一个仿真信号来验证模拟“啮合频率8000Hz 转频250Hz调制 12kHz强干扰 随机噪声”的齿轮箱场景。复制到脚本里直接跑。fs 51200; % 采样率 N 8192; % 0.16秒数据 t (0:N-1)/fs; fr 250; % 转频 fm 8000; % 啮合频率 % 幅值调制产生对称边频带±250Hz一阶边频±500Hz二阶边频 AM 1 0.15*cos(2*pi*fr*t) 0.05*cos(2*pi*2*fr*t); x AM .* cos(2*pi*fm*t); % 加入带外强干扰和随机噪声模拟现场复杂成分 x x 0.8*cos(2*pi*12000*t) 0.02*randn(size(t)); % ---- 普通FFT用1024点看窄带效果 ---- figure subplot(2,1,1) Nfft1 1024; win1 hanning(Nfft1); X1 2*abs(fft(x(1:Nfft1).*win1, Nfft1))/sum(win1); f1 (0:Nfft1/2)*fs/Nfft1; plot(f1, X1(1:Nfft1/21)) xlim([6000 10000]) title(普通FFT1024点) xlabel(频率 / Hz), ylabel(幅值 / g) % ---- 细化谱M16 ---- subplot(2,1,2) [fz, Xz] zoom_fft(x, fs, 8000, 16); plot(fz, Xz) xlim([7500 8500]) title(细化谱M16) xlabel(频率 / Hz), ylabel(幅值 / g)普通FFT里8000Hz主峰附近几乎看不到边频12kHz的强干扰旁瓣把靠近主峰的弱分量都盖住了。细化谱一出来8000Hz两侧±250Hz和±500Hz的边频带清晰可见虽然二阶边频幅度只有0.025左右但已经能从底噪里分辨出来。这个对比就是细化谱“带外滤除窄带聚焦”带来的实际效果。3.3 细化倍数M和滤波器截止怎么定细化倍数M是整个程序里最需要花心思的参数。M越大等效采样率fs/M越低显示带宽越窄频率间隔越细。但M不是越大越好的无脑参数。我的选择思路是先定需求待分析中心频率fc通常是啮合频率或啮合频率的高次谐波。需要覆盖的边频带宽度B一般取±(3~5)倍转频。比如转频250Hz想看三阶边频B至少1500Hz。理论最小细化倍数M_min ≈ fs / B实际留2倍余量M取2M_min附近。比如fs51200Hz想看fc8000Hz附近±750Hz的边频带宽B1500HzM_min≈34取M48或M32都行。取32时显示带宽是fs/M1600Hz前文代码里低通截止约0.4×1600640Hz其实只能覆盖±640Hz三阶边频750Hz会被衰减。想看到±750Hz就得把滤波器截止从0.8/M改成0.95/M或者把M降到16并接受更粗略的间隔。这个权衡要具体情况具体分析。还有一个硬约束抽取后有效点数NL/M要足够做窗函数和FFT。工程上N最好在500点以上否则窗函数主瓣太宽细化谱的谱峰本身就糊掉了。如果数据长度只有4096点M取32抽取后只剩128点这时候的细化谱不用看也知道是一堆大馒头毫无细节。数据不够长就别硬上大M要么加长采样时间要么考虑用Chirp-Z变换做非整数倍细化。4. 工程实战案例把齿轮故障边频从底噪里捞出来4.1 工况与数据背景我曾经处理过一组航空发动机附件齿轮箱的振动数据。输入轴转速15000rpm对应转频250Hz主动齿轮32齿一级啮合频率8000Hz。传感器安装在齿轮箱壳体径向采样率51200Hz每段数据0.64秒共32768点。现场频谱最闹心的是8000Hz附近。壳体上除了齿轮啮合振动还有大量别的振动成分传来轴承外圈特征频率、不平衡引起的低频分量、以及下一级齿轮的啮合谐波全都掺在一起。看普通FFT全带谱8000Hz就是一根很宽的“小土包”旁边矮矮地堆着几个峰值根本分不清哪些是齿轮边频哪些是带外分量泄漏过来的假峰。4.2 普通FFT与细化谱的对比分析先用普通FFT处理取2048点一看频率分辨率25Hz。8000Hz主峰两侧确实有东西但边频幅度太低谱线底噪又高鼠标在图上拉来拉去也不敢判断。换成细化谱M取16目标频段变成8000Hz两侧各约640Hz细化后频率间隔约1.56Hz而且12kHz以外的强分量全被滤掉谱底一下干净了很多。谱图上的8000Hz主峰两侧能清楚看到三根边频7750Hz、8250Hz以及7500Hz和8500Hz处的第二阶边频。第一阶边频幅值约为主峰的14%第二阶约6%比值低但很稳定。再往外的第三阶边频若隐若现和噪声底接近不敢下结论。这种“一阶明显、二阶可见、三阶存疑”的谱形结合时域波形里没有明显冲击判断是齿轮早期均匀磨损或轻度点蚀还没到剥落阶段。4.3 谱图判读哪些边频值得警惕诊断齿轮故障不能只看有没有边频。健康齿轮也会有边频因为制造误差和装配误差不可避免地产生调制但健康状态下边频带窄、幅度低、两侧对称。真正要警惕的是三类变化第一高阶边频出现并逐渐超过一阶边频的幅度。一般局部剥落、齿根裂纹会产生很宽的边频家族3阶、4阶甚至更多阶都能看到这是故障能量往高频扩展的表现。第二边频带不对称。如果左右两侧边频幅度差很明显往往和轴不对中、齿轮偏心这类问题相关因为转频的径向分量让调制过程不再对称。第三细化谱里单根谱线变宽。正常边频是细细一根如果细化谱上变成一簇密集的细碎谱线说明转频存在波动齿轮啮合状态不稳定。这时候别急着判故障先看看是不是转速不稳条件允许的话做角域重采样才能进一步确认。5. 常见问题与避坑实录5.1 六个高频问题的快速排查表这段时间帮几个朋友看过细化谱程序踩来踩去基本都是下面这些坑现象可能原因排查与对策细化谱频率范围和预期不符频率轴没加fc或用原始fs算了间隔检查输出ffc(0:half-1)*dfdf要用fs/M谱图边缘幅值明显跌落滤波器截止留余量太大目标频段被压把0.8/M调到0.9/M或减小显示范围幅值比普通FFT小了一大截窗函数修正没做或漏了单边乘2用sum(win)修正单边谱乘2谱底噪声很高谱峰毛糙抽取后点数太少随机误差大加长数据或用多段平均谱出现镜像假峰抽取前滤波不够抽点后混叠提高滤波阶数检查M是否过大目标频段里还是看不清边频M取太大把边频滤掉了按B倒推M确保覆盖目标边频带5.2 细化倍数越大越好别再踩这个坑了我刚写Zoom-FFT函数的时候特别迷信大M。第一次调参数直接上M128觉得倍数越大肯定越精细。结果谱图出来8000Hz主峰孤零零立着边频全没了。我把滤波器、窗函数、幅值修正查了个遍最后才发现是M128时显示带宽只有400Hz边频在7500Hz和8500Hz早就被低通滤波器削掉了。这就像拿着长焦镜头结果焦段太长想拍的目标窗户已经出了画面。从那以后我定了个规矩先画普通FFT全局谱标出关心的频率范围再算目标带宽B反推M最后看一眼抽取后剩余点数够不够。三个条件都满足才动手写细化谱参数。另外细化谱要配合普通FFT一起看先全局后局部就像先看地图再放大街道两步走才能保证不迷路。最后分享一个我自己的操作习惯。拿到一段齿轮箱振动数据我不会一上来就做细化谱而是先看时域波形、看全局FFT、定位啮合频率和几个明显的边频段心里有了大致判断再用zoom_fft逐段放大。判断故障时我还会把细化谱和包络谱配合着看细化谱看啮合频率附近的边频结构包络谱看低频段的故障特征频率两个角度互相印证。Matlab代码也别急着封装成黑盒子第一次用的时候把复调制、滤波、抽取的中间结果逐个画出来看一遍搞清楚每一步数据发生了什么变化后面换工况、调参数才不慌。这行当的经验就是这么一行一行调出来的。