简介本资源是一份面向信号处理初学者与工程实践者的Python时域频域特征提取工具脚本适用于机械故障诊断、生物电信号分析、振动监测等需量化信号特性的实际场景。资源聚焦核心特征计算逻辑覆盖6类关键时域指标方差、标准差、峭度、裕度、峰值、斜度与4类典型频域特征功率谱密度、谐波成分、带宽、中心频率全部封装于单个Python脚本中依赖numpy和scipy库开箱即用。压缩包仅含1个.py文件大小2KB轻量简洁便于嵌入项目或教学演示。目前已有1799人学习下载读者可直接运行脚本理解特征物理意义、复现计算流程、调试参数并拓展至自定义信号分析是掌握信号特征工程基础方法的实用入门参考。1. 特征提取到底在做什么1.1 从一段信号说起做设备故障诊断、语音识别、振动分析或者脑电信号处理的人几乎都绕不开“特征提取”这四个字。我最早接触这个概念是很多年前做旋转机械振动监测当时领导丢给我一段轴承振动数据说“看看能不能区分正常和故障状态”。我打开软件一看波形长得差不多肉眼根本分不出来。那时候才意识到原始信号本身很难直接当输入用必须从中提炼出一些有代表性的数值这些数值就是特征。特征提取的本质是把一段高维的、冗余的、难以直接理解的原始信号压缩成一组低维的、有区分度的、能反映系统状态的指标。打个比方你要判断一个人身体好不好不会盯着他每一秒的心电图曲线看而是看心率、血压、体温这几个数字这几个数字就是从原始生理信号里提取的特征。而在信号处理领域特征主要从两个角度去挖时域和频域。时域看的是信号幅度随时间的变化规律频域看的是信号能量在不同频率上的分布。两者看的是同一个信号但视角完全不同缺一不可。1.2 时域和频域各自解决什么问题时域特征最直观。均值告诉你信号的直流分量峰值告诉你最大冲击强度均方根值告诉你整体能量水平峭度告诉你波形是否出现尖锐的脉冲。这些指标计算量小、实时性好非常适合在线监测场景。我做过一个产线设备状态监测项目采样率20kHz每个通道每秒钟产生2万个数据点如果用频域分析每段都要做FFT计算压力不小但时域特征直接递推计算单片机都能跑得动。频域特征解决的是时域看不到的问题。一个信号里如果混着50Hz的工频干扰、200Hz的齿轮啮合频率、800Hz的轴承故障频率在时域波形里这些成分完全叠在一起你根本看不出谁是谁。但一做FFT频谱上清清楚楚出现三个峰每个峰的幅值、频率、能量都能单独量化。频域特征的最大价值就是能把混叠在一起的不同频率成分分开找到周期性故障的“指纹频率”。所以我的经验是时域特征是“粗筛”适合快速判断有没有异常频域特征是“精定位”适合找出异常来自哪里。两者结合才能构成一个完整的诊断链路。2. 时域特征公式、物理意义与选型逻辑2.1 最常用的一组时域指标先列一下我实际项目里最常用的时域特征每个都给出公式和物理含义这些不是教科书上的空泛定义而是我在现场验证过有区分度的指标。均值Mean[ \bar{x} \frac{1}{N}\sum_{i1}^{N}x_i ]代表信号的直流分量。注意一个坑很多传感器出来的信号做完隔直处理后均值理论上接近0但如果均值明显漂移说明传感器零漂或者存在低频趋势项这对后续频域分析影响很大。均方根值RMS[ x_{rms} \sqrt{\frac{1}{N}\sum_{i1}^{N}x_i^2} ]这是最常用的幅值特征反映信号的有效能量。对于振动信号RMS值和设备振动烈度直接挂钩很多ISO标准比如ISO 10816就是按RMS值来划定设备状态等级的。RMS对持续性的振动水平敏感但对瞬态冲击不敏感所以需要配合其他指标使用。峰值Peak和峰峰值Peak-to-Peak[ x_{peak} \max|x_i|, \quad x_{p-p} \max(x_i) - \min(x_i) ]峰值反映信号的最大瞬时幅值对冲击类故障很敏感但容易受随机噪声干扰。峰峰值对整体摆动幅度的刻画更稳定。实际使用时我习惯把峰值和RMS同时提取然后算峰值因子。峰值因子Crest Factor[ CF \frac{x_{peak}}{x_{rms}} ]这个指标的妙处在于它是一个无量纲的比值和信号的绝对幅值无关。正常磨损的轴承峰值因子通常在3到5之间如果出现早期剥落频谱上会出现周期性冲击峰值增长比RMS快得多峰值因子可能飙升到10以上。我做过一个案例轴承早期故障时RMS只涨了15%但峰值因子涨了一倍多这就是无量纲指标的敏感度优势。峭度Kurtosis[ K \frac{\frac{1}{N}\sum_{i1}^{N}(x_i - \bar{x})^4}{x_{rms}^4} ]峭度反映波形的尖峭程度正态分布的峭度为3。实际工程中有人直接用峭度值比如峭度大于3.5判定异常也有人用峭度减去3得到超值峭度。我在风电齿轮箱监测项目里就是靠峭度突变锁定了一个齿轮齿面点蚀故障——当时RMS完全正常但峭度从3.2跳到了7.8拆开一看齿面果然有麻点。波形因子Shape Factor、脉冲因子Impulse Factor、裕度因子Margin Factor这几个指标物理含义类似都是某种形式的峰值与RMS或绝对平均值的比值我不逐一写公式了核心思路就是通过比值关系放大冲击成分的敏感度。如果做故障诊断建议特征集里至少包含均值、RMS、峰值、峰峰值、峰值因子、峭度、偏度、波形因子、脉冲因子、裕度因子一共10个。2.2 时域特征在Matlab里的快速实现很多人拿到数据第一个问题是这些特征怎么算如果用的Matlab一段代码就能搞定。假设你的信号存在变量x里采样率是fs% 时域特征提取 N length(x); mean_val mean(x); % 均值 rms_val rms(x); % 均方根值 peak_val max(abs(x)); % 峰值 p2p_val peak2peak(x); % 峰峰值 crest_val peak_val / rms_val; % 峰值因子 kurtosis_val kurtosis(x); % 峭度 skewness_val skewness(x); % 偏度 % 波形因子 RMS / 绝对平均幅值 shape_val rms_val / mean(abs(x)); % 脉冲因子 峰值 / 绝对平均幅值 impulse_val peak_val / mean(abs(x)); % 裕度因子 峰值 / (绝对平均幅值的平方根)^2这里用方根幅值 margin_val peak_val / (mean(sqrt(abs(x))))^2; % 打包成特征向量 features_time [mean_val, rms_val, peak_val, p2p_val, ... crest_val, kurtosis_val, skewness_val, shape_val, ... impulse_val, margin_val];注意方根幅值的实际公式是[ x_r \left(\frac{1}{N}\sum_{i1}^{N}\sqrt{|x_i|}\right)^2 ]Matlab里用mean(sqrt(abs(x)))再平方就对了。这个特征在实际轴承诊断里对早期故障很敏感别漏掉。如果用的是Python核心逻辑一模一样换成numpy就行。需要注意的是这两种语言里kurtosis的定义有差异Matlab的kurtosis默认返回的是超值峭度即减去3Python的scipy.stats.kurtosis默认也减了3但pandas的kurt方法也是减过3的。不要直接拿数值去对照教科书公式3.5这种先搞清楚你用的工具到底返回的是什么定义。2.3 时域特征选型的一个判断逻辑每次做项目客户都会问特征这么多到底该用哪些我的建议是三步走。第一步看应用场景。在线监测、边缘计算这种实时性要求高的场景优先选计算量小的特征——均值、RMS、峰值、峰值因子就够了不要一上来就上全套10个指标。离线故障诊断场景可以奢侈一点把10个时域特征全部算出来。第二步做敏感性分析。拿已知状态的样本数据比如正常100组、故障100组分别计算每个特征在两类样本上的分布。如果两个分布的均值差异大、重叠区域小这个特征的区分度就高。我在实际项目里常用箱线图来可视化这个过程一眼就能看出哪个特征分得开、哪个分不开。第三步处理量纲和尺度问题。RMS、峰值的量纲是信号幅值的单位可能是m/s²、可能是mV不同通道之间没法直接比较。如果用机器学习做分类建议把特征做标准化z-score或者优先筛选无量纲特征比如峰值因子、峭度、裕度因子。这里有个容易踩的坑不要盲目堆特征。特征越多模型越容易过拟合尤其是在样本量不大的情况下。我见过一个同事从时域和频域一共提取了50多个特征全扔进随机森林结果在验证集上表现反而不如只用10个精心挑选的特征。特征的“可解释性”在工程诊断里非常重要——如果模型告诉你峭度高所以判断为故障你要能向客户解释清楚“峭度高了说明波形出现尖锐脉冲这是轴承早期故障的典型特征”。但如果你用的是某个复杂的组合特征根本没法解释客户也不敢信你的诊断结果。3. 频域特征从傅里叶变换到故障频率3.1 为什么要做频域分析时域特征说到底都是统计值丢失了信号的顺序和周期性信息。两个完全不同的波形完全可能计算出相同的RMS和峰值。但它们的频率成分一定不同——一个可能是50Hz的正弦波叠加噪声另一个可能是200Hz的方波叠加噪声时域统计值可能很像但频谱上区别巨大。频域分析的核心工具是傅里叶变换。它把一个时域信号分解成无数个不同频率、不同幅值、不同相位的正弦波的叠加。以连续信号为例傅里叶变换的定义是[ X(f) \int_{-\infty}^{\infty} x(t)e^{-j2\pi ft}dt ]但在计算机里处理的是离散信号用的是离散傅里叶变换DFT。实际工程中最常用的是快速傅里叶变换FFT是DFT的高效算法实现。Matlab里一条命令就行X fft(x); % 得到复数频谱 f (0:N-1) * fs / N; % 频率轴这里有一个极其重要的细节FFT的结果是复数幅值要取模。而且X(1)对应的是直流分量X(2)对应的是频率为fs/N的分量依此类推。很多人第一次做FFT直接plot(real(X))看到的完全不是自己想要的频谱这就是对FFT输出格式不熟悉导致的。Python的话用numpyimport numpy as np X np.fft.fft(x) f np.fft.fftfreq(N, d1/fs)3.2 功率谱和幅值谱别搞混了做频域特征首先要区分“幅值谱”和“功率谱”这两个概念这是新手最容易混淆的地方。幅值谱Amplitude Spectrum直接把FFT结果的模除以N得到物理意义是原始信号中某个频率分量的实际幅值amplitude_spectrum abs(X) / N;如果原始信号是x(t) A*sin(2*pi*f0*t)那么幅值谱在f0处的峰值就是A/2单边谱的情况下是A。不要忘记正弦信号的幅值会均等地分到正负两个频率上所以单边谱要把除直流外的幅值乘2。功率谱密度PSD表示单位频带宽度内的信号功率单位是g²/Hz或者(m/s²)²/Hz这种形式。Matlab里可以直接用pwelch函数估算[pxx, f] pwelch(x, window, noverlap, nfft, fs);功率谱的好处是抗噪能力强、平滑度好特别适合提取宽带背景下的窄带特征。在实际故障诊断中我更倾向于用功率谱因为故障特征频率往往表现为功率谱上的一个凸起比幅值谱更明显。顺带提一个热点词里出现过的概念频域OCTOptical Coherence Tomography光学相干层析的频谱域实现。医学成像领域把OCT信号从时域扫描深度方向变换到频域通过分析干涉光谱的频域信息来重构样品内部结构。核心思路和我们做振动信号频域分析完全一样都是利用傅里叶变换把信号从时域映射到频域只是在那个场景下频域编码了空间位置信息。这就是频域分析的魅力同一个数学工具在不同领域能玩出完全不同的花样。3.3 从频域里能提取哪些特征频谱不是一个特征而是一堆特征。我整理一下我实际用过的频域特征按用途分类。频率中心重心频率FC[ FC \frac{\sum_{i1}^{K}f_i \cdot P_i}{\sum_{i1}^{K}P_i} ]其中P_i是第i个频率点的功率或幅值f_i是该频率点的频率。重心频率反映信号能量在频谱上的集中位置。轴承正常情况下能量集中在低频段故障时往往会在高频段出现新的能量聚集重心频率就会往后移。均方频率MSF和频率方差VF[ MSF \frac{\sum_{i1}^{K}f_i^2 \cdot P_i}{\sum_{i1}^{K}P_i}, \quad VF MSF - FC^2 ]这两个指标描述的是频谱的分散程度。频率方差越大说明信号的能量分散在更宽的频带里通常意味着信号越复杂、越“不健康”。谱峰数和峰值频率统计频谱中超过某个阈值的峰值的数量以及最大峰值对应的频率。在齿轮箱故障诊断中啮合频率及其谐波通常是最明显的峰值如果峰值频率偏离理论值说明转速可能不稳定或存在打滑。频带能量比把频谱按照感兴趣频段划分为若干个频带计算每个频带能量占总能量的比例。比如在电机轴承诊断里低频段0-1kHz主要是转频和齿轮啮合频率中频段1-5kHz是轴承故障特征频率所在区域高频段5-20kHz是早期故障激发的结构共振区域。如果高频段能量比例突然上升那大概率是有早期故障了。这个特征比单个频谱峰更稳健因为早期故障的特征频率往往偏离理论值存在滑移但能量往高频段转移这个趋势是稳定的。边频带特征在齿轮故障诊断中故障会导致啮合频率两侧出现边频带边频带间距等于故障轴的转频。提取边频带的间距和幅值可以直接定位到是哪根轴出了问题。这个特征属于比较进阶的玩法但区分度极高。我在实际项目中一般从频域提取以下特征组成向量重心频率、均方频率、频率方差、最大峰值频率、最大峰值幅值、峰值数超过平均幅值3倍的峰个数、3个典型频带的能量比一共9个特征。3.4 频域特征提取的Matlab代码示例假设已经计算出了功率谱pxx和频率轴f下面是提取频域特征的示例% 假设 pxx 是功率谱密度f 是对应的频率轴 K length(pxx); % 重心频率 FC sum(f .* pxx) / sum(pxx); % 均方频率 MSF sum(f.^2 .* pxx) / sum(pxx); % 频率方差 VF MSF - FC^2; % 最大峰值频率和幅值 [max_pxx, idx] max(pxx); peak_freq f(idx); % 频带能量比以1kHz和5kHz为边界举例 band1 f 0 f 1000; band2 f 1000 f 5000; band3 f 5000; E_total sum(pxx); E_ratio_1 sum(pxx(band1)) / E_total; E_ratio_2 sum(pxx(band2)) / E_total; E_ratio_3 sum(pxx(band3)) / E_total; % 组装频域特征向量 features_freq [FC, MSF, VF, peak_freq, max_pxx, ... E_ratio_1, E_ratio_2, E_ratio_3];这里需要特别提醒一个工程细节如果使用pwelch返回的功率谱密度单位是每赫兹的功率如果频率轴分辨率不够高即nfft太小低频段的细节会被抹平。经验法则是 nfft 至少取信号长度的两倍或者直接用信号本身长度同时用NFFT 2^nextpow2(length(x))这种取2的幂次的做法来加速运算。3.5 如何从频域选出一个特定频率分量的幅值这是网络热词里被问烂了的问题如何在Matlab中将一组时域数据转换为频域数据并选取出某一频率。很多人在这一步卡住我详细写一下步骤。整体思路是先对时域信号做FFT构造频率轴然后找到目标频率对应的索引读取该索引处的幅值。fs 1000; % 采样率单位Hz N length(x); % 信号长度 t (0:N-1) / fs; % 时间轴 X fft(x); % FFT X_mag abs(X) / N; % 幅值谱 X_mag(2:end-1) 2 * X_mag(2:end-1); % 单边幅值谱修正 f (0:N-1) * fs / N; % 频率轴 target_freq 50; % 要提取的目标频率 [~, idx] min(abs(f - target_freq)); amplitude_at_target X_mag(idx);这里的核心是通过min(abs(f - target_freq))找到频率轴上离目标频率最近的索引。之所以要加取模和单边修正是因为fft结果是复数且能量均分在正负频率上。需要说明的是这里的频率分辨率是fs/N如果你的信号长度不够长分辨率比目标频率的精度要求还粗那这个操作就不精确了。提高分辨率的方法是增加N也就是延长采样时间而不是提高采样率。这个细节80%的新手都会搞混。4. 从时域到频域的坑非整数数据和非等间隔采样4.1 非整数频率数据怎么处理网络热词里有个问题是“如何将一组时域下非整数数据转换为频域数据”。我理解这个“非整数”通常指两种情况。第一种时间点不是整数。比如采样时刻是0.001、0.0027、0.0051这种非等间隔的时间点。这在实际测量中很常见传感器时钟抖动或者数据采集卡触发不均都会导致。直接对非等间隔的数据做FFT是行不通的因为标准FFT要求等间隔采样。处理方法是先做插值重采样把数据调整到等间隔的时间轴上再用标准FFT。Matlab里用resample或者interp1都行。第二种频率值是整数但信号长度不是2的幂次或者频率分辨率不整除。这种其实不用太担心现代FFT库Matlab和numpy对任意长度的信号都能处理只是速度上2的幂次会更快。所谓“非整数数据”很多时候只是对FFT原理不清楚导致的误解FFT并不要求信号长度是2的幂只是Cooley-Tukey算法对2的幂次的分解最方便。4.2 非等间隔重采样的实操如果你遇到的是采样时刻不均的问题我分享一个简单可靠的做法% 假设 t_interp 是等间隔时间轴x_interp 是重采样后的信号 t_uniform linspace(t(1), t(end), N); x_uniform interp1(t, x, t_uniform, pchip); X fft(x_uniform);选pchip而不是默认的linear的原因pchip保持单调性且不会出现过冲对非周期性的信号插值更平滑。如果信号本身非常干净比如实验室里的正弦信号用spline精度更高但要小心过冲问题。另一个思路是用nufftNon-Uniform FFTMatlab从R2021a开始自带这个函数可以直接对非等间隔采样的信号做FFT不用先插值X nufft(x, t, f_query);这个方法频谱泄漏特性比插值后FFT要好一些但没有插值法那么直观。对大多数工程场景我建议先插值再FFT因为整个流程的可解释性好、参数调试方便。4.3 频谱泄漏和加窗做频域分析绕不开“频谱泄漏”这个坑。理想情况下FFT的前提是信号在截取区间内是周期的。但实际采样往往截取的是非整周期长度这会导致FFT结果中能量从真实频率分散到相邻频率上看起来像频谱“漏油”了一样。解决频谱泄漏的标准方案是加窗。窗函数的作用是让信号在两端逐渐衰减到0从而削弱截断处的非连续性。常用的窗有汉宁窗Hann通用分析首选主瓣适中、旁瓣衰减快汉明窗Hamming频率分辨率略好于汉宁但旁瓣衰减慢布莱克曼窗Blackman旁瓣衰减很大适合检测幅度差异悬殊的相邻频率分量平顶窗Flattop幅值测量精度最高适合精确测量单频幅值但主瓣很宽实际选择逻辑如果关心的是频率定位精度用汉宁窗如果关心的是幅值测量精度用平顶窗如果信号里有强弱差异悬殊的频率成分用布莱克曼窗。一个最基础的建议默认用汉宁窗遇到问题再调。加窗会导致幅值衰减所以做幅值谱时要做幅值恢复修正。Matlab里可以直接用自带函数pwelch内部已经处理了窗的幅值修正不需要手动校正。但如果你手动X fft(x .* hann(N))再幅值修正就需要乘以一个系数一般是除以窗函数的均值w hann(N); X fft(x .* w); X_mag abs(X) / sum(w) * 2; % 单边幅值修正4.4 为什么FFT结果和想象的不一样实操中每个人都会遇到“FFT结果看起来不对”的时刻。排查顺序我建议这样来第一步看基线。FFT结果的第一个点索引1是直流分量如果信号有直流偏置这里会有一个巨大的峰值把整个频谱的形状压扁。解决办法是先去均值x x - mean(x)。第二步看幅值。如果信号是一个幅值为1的50Hz正弦波FFT幅值谱在50Hz处的峰值应该是1单边修正后但如果忘了修正或者加窗了没修正结果可能是0.5或者更小。第三步看频率轴。检查频率轴的单位和取值范围fftfreq返回的是以循环/秒为单位的频率要确认是不是乘以了采样率。第四步看对称性。双边频谱关于fs/2对称只取前半段才是单边频谱。第五步看泄漏。如果频谱峰值旁边出现了很多小旁瓣说明存在频谱泄漏考虑加窗或者增加采样长度。5. 特征提取实战案例滚动轴承故障诊断5.1 案例背景与数据准备用一个我做过的最典型的例子来串起整个流程滚动轴承故障诊断。轴承是旋转机械里最容易损坏的部件之一故障信号早期特征非常微弱被淹没在齿轮噪声和背景振动里时域特征不明显必须靠频域特征和包络分析来捕捉。假设数据来自一个驱动端轴承座上的加速度传感器采样率 fs 12kHz这是凯斯西储大学轴承数据中心的标准配置也是很多人入门用的经典数据集内圈故障特征频率约为 157.94Hz外圈约 104.56Hz滚动体约 68.93Hz保持架约 9.94Hz。转频约29.5Hz。数据分成两类正常状态和内圈故障状态。每段信号取1024个点大概相当于0.085秒在轴承工频29.5Hz下能覆盖2.5个转周期足够反映一个周期的故障冲击。5.2 完整特征提取流程整个流程分成三步。第一步时域特征提取。对每段数据计算10个时域特征。重点关注的指标是RMS和峭度。正常状态RMS大约0.02到0.03内圈故障状态RMS可能到0.09到0.15峭度正常状态接近3故障状态往往在5到10之间。第二步频域分析。先对原始信号做FFT看频谱你会发现直接FFT的频谱中故障特征频率处的峰值并不明显因为故障冲击激发的共振频率通常很高几kHz到十几kHz但能量分散。这时候需要用包络分析先对信号做带通滤波比如8kHz到12kHz然后取希尔伯特变换的模得到包络信号再对包络信号做FFT故障特征频率就会在低频段清晰可见。这个技术叫包络解调分析是轴承诊断的黄金标准。包络分析的核心代码% 带通滤波 [b, a] butter(4, [8000 12000] / (fs/2), bandpass); x_filtered filtfilt(b, a, x); % 希尔伯特变换求包络 envelope abs(hilbert(x_filtered)); % 对包络做FFT N_env length(envelope); X_env fft(envelope); X_env_mag abs(X_env) / N_env; f_env (0:N_env-1) * fs / N_env; % 只看0-1000Hz段 idx f_env 0 f_env 1000; plot(f_env(idx), X_env_mag(idx));第三步提取频域特征。在包络谱上找理论故障频率附近的峰值幅值比如在157.94Hz附近找最大峰值这个幅值就是内圈故障的强度指标。同时计算重心频率、均方频率等全局特征。5.3 特征组合带来的分类效果把时域频域的特征合并成特征向量后我用随机森林做二分类正常/故障交叉验证准确率能达到99%以上。但如果只用时域特征准确率大约90%只用频域特征大约95%。两者的组合效果最好这就是时频特征结合的价值所在。这里顺便提一嘴MFCC特征。MFCC梅尔频率倒谱系数是语音识别里的经典特征流程是分帧→加窗→FFT→梅尔滤波器组→取对数→DCT→得到倒谱系数。本质上就是在频域特征上再做一次压缩变换。如果你做的是音频信号的特征提取MFCC基本是绕不开的方案但做机械振动诊断时MFCC用得少因为梅尔滤波器组的频率刻度是为听觉设计的不适合机械故障频率分析。这个例子说明特征提取方法必须和应用场景匹配不能盲目套用。6. 常见问题与排查技巧实录6.1 时域特征看起来没区别怎么办这是最常遇到的问题正常和故障状态下提取出来的特征值差异极小完全区分不开。我遇到这种问题第一反应是看数据本身有没有问题。先用时域波形画出来看看是否存在削波传感器量程不够导致波形顶部被截平、断线、非线性漂移。削波会让峰值失真峰值因子和峭度都会异常。数据没问题的话考虑是不是特征选的阶段不对。早期微弱故障RMS和峰值都变化很小但高频段能量和峭度可能已经有反应。这时不要死盯着时域特征直接上包络谱分析。还有一种情况是采样长度不够。特征提取对样本长度有要求特别是峭度这种高阶统计量数据太短比如只有几十个点估计值方差巨大根本不可靠。我建议单段特征提取的样本长度至少覆盖10个转频周期做轴承诊断的话每段至少1024个点最好4096个点。6.2 频谱峰值偏移对不上理论值故障特征频率应该在某一个频率处出现峰值实测发现峰值偏了这可能有两种情况。第一种转速波动。理论故障频率是假设转速恒定算出来的实际设备转速可能有波动负载变化、电源波动导致故障频率偏移。解决方法是同时采集转速信号用转速信号做等角度重采样角域平均或者把故障频率按照实际转速归一化到转频倍数。第二种频率分辨率不够。比如ls12kHzN1024频率分辨率就是fs/N ≈ 11.7Hz而内圈故障频率和滚动体故障频率的差可能不到10Hz频谱上两个峰值合成一个根本分辨不出来。解决方法是增加样本长度N分辨率变细但代价是计算时间变长。顺带提一下“方位向压缩是压缩频域还是时域”这个问题。这是合成孔径雷达SAR成像里的概念方位向压缩确实是在频域多普勒频域做的处理——通过匹配滤波对方位向信号进行压缩提高方位向分辨率。本质上和我们在信号处理里做FFT后滤波再IFFT回到时域是一回事。这说明了频域处理在很多领域都是绕不开的核心操作。6.3 FFT结果幅值来回跳不稳定同一段信号重复做FFT幅值应该是一样的因为FFT是确定性计算。如果你发现幅值不稳定大概率是信号本身不稳定不是FFT的问题。一个常见原因是数据采集的触发条件不一致导致每帧数据的起始相位不同。对非相参信号这会导致FFT幅值出现波动虽然功率谱密度会稳定一些。解决方案是用重叠率更高的窗函数增大overlap或者用Welch方法的平均功率谱来替代单帧FFT。另外检查一下是否有瞬态干扰混入。比如轴承故障信号的冲击性很强如果某帧数据恰好包含了一个大冲击RMS和峰值都会偏高。处理方法是做中值滤波或者异常值剔除。6.4 特征矩阵里出现NaN怎么办做特征批量提取时如果某个特征出现了NaN通常原因是分母为0。比如波形因子是RMS除以绝对平均幅值如果某帧数据全为0传感器没输出分母为0结果就是NaN或者Inf。这个一定要在代码里加保护if abs(mean(abs(x))) eps shape_val 0; % 或者标记为异常样本 else shape_val rms_val / mean(abs(x)); endMatlab里还可以用isnan检查最终特征矩阵把含NaN的行剔除掉。但要注意如果剔除的比例很高比如超过5%说明数据质量有问题先排查传感器和数据采集系统而不是简单粗暴地删除样本。6.5 时域掩蔽效应和滚动时域优化是什么网络热词里提到了时域掩蔽效应和滚动时域优化这两个词容易混在一起但完全是两个领域的概念。时域掩蔽效应是心理声学听觉感知里的概念当一个强的声音出现后它会对后面紧接着出现的弱声音产生掩蔽导致弱声音在一段时间内听不见。这在音频编码比如MP3、AAC中被利用通过心理声学模型把听不见的信息丢掉来压缩数据。如果你做音频处理这个效应会影响特征提取——有些频谱细节虽然真实存在但人耳听不见在特征中是否能被利用取决于你是做人耳感知相关的任务还是机器诊断任务。滚动时域优化Receding Horizon Optimization是控制理论里的概念核心思路是每个控制周期都基于当前状态做一个有限时域的最优控制问题执行第一步然后滚动到下一个周期重新优化。这个名字里虽然有“时域”两个字但它和时域特征提取完全不是一个层面的事。这两个词提醒我们搜索时域相关的内容时要注意区分“时域分析”在信号处理和在其他领域的含义避免概念混用。6.6 特征提取的速度优化如果特征是用于在线监测的性能优化很重要。讲几个实测有效的优化技巧。技巧一是避免在循环里逐点操作。用numpy或Matlab的向量化计算比在循环里逐点加要快几十倍。写特征提取脚本先考虑有没有现成的向量化函数比如rms、kurtosis、pwelch这些内部都是高度优化的比自己手写循环快得多。技巧二是合理选择FFT长度。不是越长越好长FFT计算慢而且对实时性要求高时没必要。如果只需要关心某个频段的特征可以先带通滤波再降采样然后再做FFT计算量可以大幅降低。技巧三是用缓存复用中间结果。比如特征提取里RMS、峰值、峭度都要用到信号幅值如果对数据做了归一化处理可以先缓存归一化后的数据避免重复计算。技巧四是用流式计算处理时域特征。均值、RMS这些统计量可以写成递推形式每次进来一个新数据点就更新一次不需要缓存整段数据。对于嵌入式设备这种流式特征提取方式几乎是唯一选择。7. 特征提取的工程落地经验分享7.1 一套完整的特征提取代码框架最后分享一个我实际在项目里用的特征提取框架把时域和频域特征合并成一个函数输入一段信号输出完整的特征向量。function feats extract_features(x, fs) % 输入: x - 一段信号, fs - 采样率 % 输出: feats - 特征向量依次为时域10个、频域9个 % 时域特征 N length(x); mean_val mean(x); rms_val rms(x); peak_val max(abs(x)); p2p_val peak2peak(x); crest_val peak_val / rms_val; kurtosis_val kurtosis(x); skewness_val skewness(x); shape_val rms_val / (sum(abs(x)) / N); impulse_val peak_val / (sum(abs(x)) / N); margin_val peak_val / (mean(sqrt(abs(x))))^2; feats_time [mean_val, rms_val, peak_val, p2p_val, ... crest_val, kurtosis_val, skewness_val, ... shape_val, impulse_val, margin_val]; % 频域特征 % 去均值避免直流干扰 x x - mean(x); % 功率谱估计用汉宁窗50%重叠 window hann(256); noverlap 128; nfft 512; [pxx, f] pwelch(x, window, noverlap, nfft, fs); % 只分析正频率范围 pos_idx f 0; f f(pos_idx); pxx pxx(pos_idx); K length(pxx); FC sum(f .* pxx) / sum(pxx); MSF sum(f.^2 .* pxx) / sum(pxx); VF MSF - FC^2; [max_pxx, max_idx] max(pxx); peak_freq f(max_idx); % 频带能量比这里用的6个等宽频带 band_ratios zeros(1, 6); band_edges linspace(0, max(f), 7); E_total sum(pxx); for k 1:6 band_mask f band_edges(k) f band_edges(k1); band_ratios(k) sum(pxx(band_mask)) / E_total; end feats_freq [FC, MSF, VF, peak_freq, max_pxx, band_ratios]; % 合并 feats [feats_time, feats_freq]; end这个函数可以直接用来批量处理数据。需要注意一点pwelch的窗口长度、nfft这些参数要根据你的采样率和目的频率范围调整如果故障特征频率在100Hz附近窗口别设太短否则频率分辨率不够pwelch会把低频细节抹掉。我这里的256窗口在12kHz采样率下对应约47Hz分辨率对轴承故障分析是够用的但如果做齿轮箱分析特征频率通常在1kHz以上窗口可以设更短来换取时间分辨率。7.2 特征提取和机器学习的衔接特征提取得再好最终还是要喂给分类器或者回归模型。这里有几个实操建议。特征矩阵的格式每一行是一个样本每一列是一个特征。拿正常样本30组、故障样本30组提取特征后得到60行×19列的特征矩阵加上标签列一共60×20。特征标准化很多机器学习算法对特征的尺度敏感比如SVM和KNN。建议在送入模型前做标准化在Python里用StandardScaler在Matlab里用zscore。标准化系数均值和标准差必须只用训练集计算然后应用到测试集否则会造成数据泄漏导致模型评估结果虚高。特征选择如果特征维度高可以用随机森林的特征重要性排序、递归特征消除RFE等方法筛掉不重要的特征。但对于故障诊断这种场景我更推荐保留结果可解释性强的特征哪怕模型准确率稍微低一点客户能看懂为什么报警价值更大。模型选择样本量不大几十到几百组的情况下随机森林和梯度提升树是默认选择不需要神经网络也能达到很好的效果。样本量上几千组之后可以尝试一维CNN直接从原始信号提取特征这就进入端到端学习的范畴了但可解释性会大幅下降。7.3 最后再提醒几个实践中的经验根据我的经验时域频域特征提取这个领域最核心的能力不是会套公式而是建立“信号→特征→状态”的完整链路。我见过太多人卡在中间某一步——或者FFT的幅值没搞清楚或者特征选了一堆但没法解释导致整个项目交付不了。几个最后的补充提醒关于采样率特征提取之前先确认采样率是否满足你的分析需求。采样定理要求采样率至少大于最高关注频率的两倍实际工程建议留3到5倍裕量。做轴承故障诊断如果故障特征频率在10kHz附近采样率至少要20kHz最好40kHz以上。关于样本长度时域统计特征需要足够长的数据才能稳定估计。峭度这种高阶统计量特别依赖于样本长度。经验法则峭度至少要1024个点才可靠最好4096以上。关于参考标准很多行业都有关于振动特征的评判标准比如ISO 10816用RMS做设备状态分级ISO 7919用轴振动位移做评判。做工程诊断时先查一下你所在行业有没有对应标准而不是自己定一个阈值。关于数据记录做特征提取前把传感器类型、采样率、量程、安装位置、设备转速、负载情况全部记录下来。同样的特征值在不同工况下含义完全不同。我做过一个项目现场工程师把传感器装反了方向导致采集的所有数据极性反转特征值全乱了——这种问题单纯靠算法永远发现不了必须有完整的现场记录来辅助判断。本文还有配套的精品资源点击获取 信号特征提取实战：时域频域指标详解与故障诊断应用
简介本资源是一份面向信号处理初学者与工程实践者的Python时域频域特征提取工具脚本适用于机械故障诊断、生物电信号分析、振动监测等需量化信号特性的实际场景。资源聚焦核心特征计算逻辑覆盖6类关键时域指标方差、标准差、峭度、裕度、峰值、斜度与4类典型频域特征功率谱密度、谐波成分、带宽、中心频率全部封装于单个Python脚本中依赖numpy和scipy库开箱即用。压缩包仅含1个.py文件大小2KB轻量简洁便于嵌入项目或教学演示。目前已有1799人学习下载读者可直接运行脚本理解特征物理意义、复现计算流程、调试参数并拓展至自定义信号分析是掌握信号特征工程基础方法的实用入门参考。1. 特征提取到底在做什么1.1 从一段信号说起做设备故障诊断、语音识别、振动分析或者脑电信号处理的人几乎都绕不开“特征提取”这四个字。我最早接触这个概念是很多年前做旋转机械振动监测当时领导丢给我一段轴承振动数据说“看看能不能区分正常和故障状态”。我打开软件一看波形长得差不多肉眼根本分不出来。那时候才意识到原始信号本身很难直接当输入用必须从中提炼出一些有代表性的数值这些数值就是特征。特征提取的本质是把一段高维的、冗余的、难以直接理解的原始信号压缩成一组低维的、有区分度的、能反映系统状态的指标。打个比方你要判断一个人身体好不好不会盯着他每一秒的心电图曲线看而是看心率、血压、体温这几个数字这几个数字就是从原始生理信号里提取的特征。而在信号处理领域特征主要从两个角度去挖时域和频域。时域看的是信号幅度随时间的变化规律频域看的是信号能量在不同频率上的分布。两者看的是同一个信号但视角完全不同缺一不可。1.2 时域和频域各自解决什么问题时域特征最直观。均值告诉你信号的直流分量峰值告诉你最大冲击强度均方根值告诉你整体能量水平峭度告诉你波形是否出现尖锐的脉冲。这些指标计算量小、实时性好非常适合在线监测场景。我做过一个产线设备状态监测项目采样率20kHz每个通道每秒钟产生2万个数据点如果用频域分析每段都要做FFT计算压力不小但时域特征直接递推计算单片机都能跑得动。频域特征解决的是时域看不到的问题。一个信号里如果混着50Hz的工频干扰、200Hz的齿轮啮合频率、800Hz的轴承故障频率在时域波形里这些成分完全叠在一起你根本看不出谁是谁。但一做FFT频谱上清清楚楚出现三个峰每个峰的幅值、频率、能量都能单独量化。频域特征的最大价值就是能把混叠在一起的不同频率成分分开找到周期性故障的“指纹频率”。所以我的经验是时域特征是“粗筛”适合快速判断有没有异常频域特征是“精定位”适合找出异常来自哪里。两者结合才能构成一个完整的诊断链路。2. 时域特征公式、物理意义与选型逻辑2.1 最常用的一组时域指标先列一下我实际项目里最常用的时域特征每个都给出公式和物理含义这些不是教科书上的空泛定义而是我在现场验证过有区分度的指标。均值Mean[ \bar{x} \frac{1}{N}\sum_{i1}^{N}x_i ]代表信号的直流分量。注意一个坑很多传感器出来的信号做完隔直处理后均值理论上接近0但如果均值明显漂移说明传感器零漂或者存在低频趋势项这对后续频域分析影响很大。均方根值RMS[ x_{rms} \sqrt{\frac{1}{N}\sum_{i1}^{N}x_i^2} ]这是最常用的幅值特征反映信号的有效能量。对于振动信号RMS值和设备振动烈度直接挂钩很多ISO标准比如ISO 10816就是按RMS值来划定设备状态等级的。RMS对持续性的振动水平敏感但对瞬态冲击不敏感所以需要配合其他指标使用。峰值Peak和峰峰值Peak-to-Peak[ x_{peak} \max|x_i|, \quad x_{p-p} \max(x_i) - \min(x_i) ]峰值反映信号的最大瞬时幅值对冲击类故障很敏感但容易受随机噪声干扰。峰峰值对整体摆动幅度的刻画更稳定。实际使用时我习惯把峰值和RMS同时提取然后算峰值因子。峰值因子Crest Factor[ CF \frac{x_{peak}}{x_{rms}} ]这个指标的妙处在于它是一个无量纲的比值和信号的绝对幅值无关。正常磨损的轴承峰值因子通常在3到5之间如果出现早期剥落频谱上会出现周期性冲击峰值增长比RMS快得多峰值因子可能飙升到10以上。我做过一个案例轴承早期故障时RMS只涨了15%但峰值因子涨了一倍多这就是无量纲指标的敏感度优势。峭度Kurtosis[ K \frac{\frac{1}{N}\sum_{i1}^{N}(x_i - \bar{x})^4}{x_{rms}^4} ]峭度反映波形的尖峭程度正态分布的峭度为3。实际工程中有人直接用峭度值比如峭度大于3.5判定异常也有人用峭度减去3得到超值峭度。我在风电齿轮箱监测项目里就是靠峭度突变锁定了一个齿轮齿面点蚀故障——当时RMS完全正常但峭度从3.2跳到了7.8拆开一看齿面果然有麻点。波形因子Shape Factor、脉冲因子Impulse Factor、裕度因子Margin Factor这几个指标物理含义类似都是某种形式的峰值与RMS或绝对平均值的比值我不逐一写公式了核心思路就是通过比值关系放大冲击成分的敏感度。如果做故障诊断建议特征集里至少包含均值、RMS、峰值、峰峰值、峰值因子、峭度、偏度、波形因子、脉冲因子、裕度因子一共10个。2.2 时域特征在Matlab里的快速实现很多人拿到数据第一个问题是这些特征怎么算如果用的Matlab一段代码就能搞定。假设你的信号存在变量x里采样率是fs% 时域特征提取 N length(x); mean_val mean(x); % 均值 rms_val rms(x); % 均方根值 peak_val max(abs(x)); % 峰值 p2p_val peak2peak(x); % 峰峰值 crest_val peak_val / rms_val; % 峰值因子 kurtosis_val kurtosis(x); % 峭度 skewness_val skewness(x); % 偏度 % 波形因子 RMS / 绝对平均幅值 shape_val rms_val / mean(abs(x)); % 脉冲因子 峰值 / 绝对平均幅值 impulse_val peak_val / mean(abs(x)); % 裕度因子 峰值 / (绝对平均幅值的平方根)^2这里用方根幅值 margin_val peak_val / (mean(sqrt(abs(x))))^2; % 打包成特征向量 features_time [mean_val, rms_val, peak_val, p2p_val, ... crest_val, kurtosis_val, skewness_val, shape_val, ... impulse_val, margin_val];注意方根幅值的实际公式是[ x_r \left(\frac{1}{N}\sum_{i1}^{N}\sqrt{|x_i|}\right)^2 ]Matlab里用mean(sqrt(abs(x)))再平方就对了。这个特征在实际轴承诊断里对早期故障很敏感别漏掉。如果用的是Python核心逻辑一模一样换成numpy就行。需要注意的是这两种语言里kurtosis的定义有差异Matlab的kurtosis默认返回的是超值峭度即减去3Python的scipy.stats.kurtosis默认也减了3但pandas的kurt方法也是减过3的。不要直接拿数值去对照教科书公式3.5这种先搞清楚你用的工具到底返回的是什么定义。2.3 时域特征选型的一个判断逻辑每次做项目客户都会问特征这么多到底该用哪些我的建议是三步走。第一步看应用场景。在线监测、边缘计算这种实时性要求高的场景优先选计算量小的特征——均值、RMS、峰值、峰值因子就够了不要一上来就上全套10个指标。离线故障诊断场景可以奢侈一点把10个时域特征全部算出来。第二步做敏感性分析。拿已知状态的样本数据比如正常100组、故障100组分别计算每个特征在两类样本上的分布。如果两个分布的均值差异大、重叠区域小这个特征的区分度就高。我在实际项目里常用箱线图来可视化这个过程一眼就能看出哪个特征分得开、哪个分不开。第三步处理量纲和尺度问题。RMS、峰值的量纲是信号幅值的单位可能是m/s²、可能是mV不同通道之间没法直接比较。如果用机器学习做分类建议把特征做标准化z-score或者优先筛选无量纲特征比如峰值因子、峭度、裕度因子。这里有个容易踩的坑不要盲目堆特征。特征越多模型越容易过拟合尤其是在样本量不大的情况下。我见过一个同事从时域和频域一共提取了50多个特征全扔进随机森林结果在验证集上表现反而不如只用10个精心挑选的特征。特征的“可解释性”在工程诊断里非常重要——如果模型告诉你峭度高所以判断为故障你要能向客户解释清楚“峭度高了说明波形出现尖锐脉冲这是轴承早期故障的典型特征”。但如果你用的是某个复杂的组合特征根本没法解释客户也不敢信你的诊断结果。3. 频域特征从傅里叶变换到故障频率3.1 为什么要做频域分析时域特征说到底都是统计值丢失了信号的顺序和周期性信息。两个完全不同的波形完全可能计算出相同的RMS和峰值。但它们的频率成分一定不同——一个可能是50Hz的正弦波叠加噪声另一个可能是200Hz的方波叠加噪声时域统计值可能很像但频谱上区别巨大。频域分析的核心工具是傅里叶变换。它把一个时域信号分解成无数个不同频率、不同幅值、不同相位的正弦波的叠加。以连续信号为例傅里叶变换的定义是[ X(f) \int_{-\infty}^{\infty} x(t)e^{-j2\pi ft}dt ]但在计算机里处理的是离散信号用的是离散傅里叶变换DFT。实际工程中最常用的是快速傅里叶变换FFT是DFT的高效算法实现。Matlab里一条命令就行X fft(x); % 得到复数频谱 f (0:N-1) * fs / N; % 频率轴这里有一个极其重要的细节FFT的结果是复数幅值要取模。而且X(1)对应的是直流分量X(2)对应的是频率为fs/N的分量依此类推。很多人第一次做FFT直接plot(real(X))看到的完全不是自己想要的频谱这就是对FFT输出格式不熟悉导致的。Python的话用numpyimport numpy as np X np.fft.fft(x) f np.fft.fftfreq(N, d1/fs)3.2 功率谱和幅值谱别搞混了做频域特征首先要区分“幅值谱”和“功率谱”这两个概念这是新手最容易混淆的地方。幅值谱Amplitude Spectrum直接把FFT结果的模除以N得到物理意义是原始信号中某个频率分量的实际幅值amplitude_spectrum abs(X) / N;如果原始信号是x(t) A*sin(2*pi*f0*t)那么幅值谱在f0处的峰值就是A/2单边谱的情况下是A。不要忘记正弦信号的幅值会均等地分到正负两个频率上所以单边谱要把除直流外的幅值乘2。功率谱密度PSD表示单位频带宽度内的信号功率单位是g²/Hz或者(m/s²)²/Hz这种形式。Matlab里可以直接用pwelch函数估算[pxx, f] pwelch(x, window, noverlap, nfft, fs);功率谱的好处是抗噪能力强、平滑度好特别适合提取宽带背景下的窄带特征。在实际故障诊断中我更倾向于用功率谱因为故障特征频率往往表现为功率谱上的一个凸起比幅值谱更明显。顺带提一个热点词里出现过的概念频域OCTOptical Coherence Tomography光学相干层析的频谱域实现。医学成像领域把OCT信号从时域扫描深度方向变换到频域通过分析干涉光谱的频域信息来重构样品内部结构。核心思路和我们做振动信号频域分析完全一样都是利用傅里叶变换把信号从时域映射到频域只是在那个场景下频域编码了空间位置信息。这就是频域分析的魅力同一个数学工具在不同领域能玩出完全不同的花样。3.3 从频域里能提取哪些特征频谱不是一个特征而是一堆特征。我整理一下我实际用过的频域特征按用途分类。频率中心重心频率FC[ FC \frac{\sum_{i1}^{K}f_i \cdot P_i}{\sum_{i1}^{K}P_i} ]其中P_i是第i个频率点的功率或幅值f_i是该频率点的频率。重心频率反映信号能量在频谱上的集中位置。轴承正常情况下能量集中在低频段故障时往往会在高频段出现新的能量聚集重心频率就会往后移。均方频率MSF和频率方差VF[ MSF \frac{\sum_{i1}^{K}f_i^2 \cdot P_i}{\sum_{i1}^{K}P_i}, \quad VF MSF - FC^2 ]这两个指标描述的是频谱的分散程度。频率方差越大说明信号的能量分散在更宽的频带里通常意味着信号越复杂、越“不健康”。谱峰数和峰值频率统计频谱中超过某个阈值的峰值的数量以及最大峰值对应的频率。在齿轮箱故障诊断中啮合频率及其谐波通常是最明显的峰值如果峰值频率偏离理论值说明转速可能不稳定或存在打滑。频带能量比把频谱按照感兴趣频段划分为若干个频带计算每个频带能量占总能量的比例。比如在电机轴承诊断里低频段0-1kHz主要是转频和齿轮啮合频率中频段1-5kHz是轴承故障特征频率所在区域高频段5-20kHz是早期故障激发的结构共振区域。如果高频段能量比例突然上升那大概率是有早期故障了。这个特征比单个频谱峰更稳健因为早期故障的特征频率往往偏离理论值存在滑移但能量往高频段转移这个趋势是稳定的。边频带特征在齿轮故障诊断中故障会导致啮合频率两侧出现边频带边频带间距等于故障轴的转频。提取边频带的间距和幅值可以直接定位到是哪根轴出了问题。这个特征属于比较进阶的玩法但区分度极高。我在实际项目中一般从频域提取以下特征组成向量重心频率、均方频率、频率方差、最大峰值频率、最大峰值幅值、峰值数超过平均幅值3倍的峰个数、3个典型频带的能量比一共9个特征。3.4 频域特征提取的Matlab代码示例假设已经计算出了功率谱pxx和频率轴f下面是提取频域特征的示例% 假设 pxx 是功率谱密度f 是对应的频率轴 K length(pxx); % 重心频率 FC sum(f .* pxx) / sum(pxx); % 均方频率 MSF sum(f.^2 .* pxx) / sum(pxx); % 频率方差 VF MSF - FC^2; % 最大峰值频率和幅值 [max_pxx, idx] max(pxx); peak_freq f(idx); % 频带能量比以1kHz和5kHz为边界举例 band1 f 0 f 1000; band2 f 1000 f 5000; band3 f 5000; E_total sum(pxx); E_ratio_1 sum(pxx(band1)) / E_total; E_ratio_2 sum(pxx(band2)) / E_total; E_ratio_3 sum(pxx(band3)) / E_total; % 组装频域特征向量 features_freq [FC, MSF, VF, peak_freq, max_pxx, ... E_ratio_1, E_ratio_2, E_ratio_3];这里需要特别提醒一个工程细节如果使用pwelch返回的功率谱密度单位是每赫兹的功率如果频率轴分辨率不够高即nfft太小低频段的细节会被抹平。经验法则是 nfft 至少取信号长度的两倍或者直接用信号本身长度同时用NFFT 2^nextpow2(length(x))这种取2的幂次的做法来加速运算。3.5 如何从频域选出一个特定频率分量的幅值这是网络热词里被问烂了的问题如何在Matlab中将一组时域数据转换为频域数据并选取出某一频率。很多人在这一步卡住我详细写一下步骤。整体思路是先对时域信号做FFT构造频率轴然后找到目标频率对应的索引读取该索引处的幅值。fs 1000; % 采样率单位Hz N length(x); % 信号长度 t (0:N-1) / fs; % 时间轴 X fft(x); % FFT X_mag abs(X) / N; % 幅值谱 X_mag(2:end-1) 2 * X_mag(2:end-1); % 单边幅值谱修正 f (0:N-1) * fs / N; % 频率轴 target_freq 50; % 要提取的目标频率 [~, idx] min(abs(f - target_freq)); amplitude_at_target X_mag(idx);这里的核心是通过min(abs(f - target_freq))找到频率轴上离目标频率最近的索引。之所以要加取模和单边修正是因为fft结果是复数且能量均分在正负频率上。需要说明的是这里的频率分辨率是fs/N如果你的信号长度不够长分辨率比目标频率的精度要求还粗那这个操作就不精确了。提高分辨率的方法是增加N也就是延长采样时间而不是提高采样率。这个细节80%的新手都会搞混。4. 从时域到频域的坑非整数数据和非等间隔采样4.1 非整数频率数据怎么处理网络热词里有个问题是“如何将一组时域下非整数数据转换为频域数据”。我理解这个“非整数”通常指两种情况。第一种时间点不是整数。比如采样时刻是0.001、0.0027、0.0051这种非等间隔的时间点。这在实际测量中很常见传感器时钟抖动或者数据采集卡触发不均都会导致。直接对非等间隔的数据做FFT是行不通的因为标准FFT要求等间隔采样。处理方法是先做插值重采样把数据调整到等间隔的时间轴上再用标准FFT。Matlab里用resample或者interp1都行。第二种频率值是整数但信号长度不是2的幂次或者频率分辨率不整除。这种其实不用太担心现代FFT库Matlab和numpy对任意长度的信号都能处理只是速度上2的幂次会更快。所谓“非整数数据”很多时候只是对FFT原理不清楚导致的误解FFT并不要求信号长度是2的幂只是Cooley-Tukey算法对2的幂次的分解最方便。4.2 非等间隔重采样的实操如果你遇到的是采样时刻不均的问题我分享一个简单可靠的做法% 假设 t_interp 是等间隔时间轴x_interp 是重采样后的信号 t_uniform linspace(t(1), t(end), N); x_uniform interp1(t, x, t_uniform, pchip); X fft(x_uniform);选pchip而不是默认的linear的原因pchip保持单调性且不会出现过冲对非周期性的信号插值更平滑。如果信号本身非常干净比如实验室里的正弦信号用spline精度更高但要小心过冲问题。另一个思路是用nufftNon-Uniform FFTMatlab从R2021a开始自带这个函数可以直接对非等间隔采样的信号做FFT不用先插值X nufft(x, t, f_query);这个方法频谱泄漏特性比插值后FFT要好一些但没有插值法那么直观。对大多数工程场景我建议先插值再FFT因为整个流程的可解释性好、参数调试方便。4.3 频谱泄漏和加窗做频域分析绕不开“频谱泄漏”这个坑。理想情况下FFT的前提是信号在截取区间内是周期的。但实际采样往往截取的是非整周期长度这会导致FFT结果中能量从真实频率分散到相邻频率上看起来像频谱“漏油”了一样。解决频谱泄漏的标准方案是加窗。窗函数的作用是让信号在两端逐渐衰减到0从而削弱截断处的非连续性。常用的窗有汉宁窗Hann通用分析首选主瓣适中、旁瓣衰减快汉明窗Hamming频率分辨率略好于汉宁但旁瓣衰减慢布莱克曼窗Blackman旁瓣衰减很大适合检测幅度差异悬殊的相邻频率分量平顶窗Flattop幅值测量精度最高适合精确测量单频幅值但主瓣很宽实际选择逻辑如果关心的是频率定位精度用汉宁窗如果关心的是幅值测量精度用平顶窗如果信号里有强弱差异悬殊的频率成分用布莱克曼窗。一个最基础的建议默认用汉宁窗遇到问题再调。加窗会导致幅值衰减所以做幅值谱时要做幅值恢复修正。Matlab里可以直接用自带函数pwelch内部已经处理了窗的幅值修正不需要手动校正。但如果你手动X fft(x .* hann(N))再幅值修正就需要乘以一个系数一般是除以窗函数的均值w hann(N); X fft(x .* w); X_mag abs(X) / sum(w) * 2; % 单边幅值修正4.4 为什么FFT结果和想象的不一样实操中每个人都会遇到“FFT结果看起来不对”的时刻。排查顺序我建议这样来第一步看基线。FFT结果的第一个点索引1是直流分量如果信号有直流偏置这里会有一个巨大的峰值把整个频谱的形状压扁。解决办法是先去均值x x - mean(x)。第二步看幅值。如果信号是一个幅值为1的50Hz正弦波FFT幅值谱在50Hz处的峰值应该是1单边修正后但如果忘了修正或者加窗了没修正结果可能是0.5或者更小。第三步看频率轴。检查频率轴的单位和取值范围fftfreq返回的是以循环/秒为单位的频率要确认是不是乘以了采样率。第四步看对称性。双边频谱关于fs/2对称只取前半段才是单边频谱。第五步看泄漏。如果频谱峰值旁边出现了很多小旁瓣说明存在频谱泄漏考虑加窗或者增加采样长度。5. 特征提取实战案例滚动轴承故障诊断5.1 案例背景与数据准备用一个我做过的最典型的例子来串起整个流程滚动轴承故障诊断。轴承是旋转机械里最容易损坏的部件之一故障信号早期特征非常微弱被淹没在齿轮噪声和背景振动里时域特征不明显必须靠频域特征和包络分析来捕捉。假设数据来自一个驱动端轴承座上的加速度传感器采样率 fs 12kHz这是凯斯西储大学轴承数据中心的标准配置也是很多人入门用的经典数据集内圈故障特征频率约为 157.94Hz外圈约 104.56Hz滚动体约 68.93Hz保持架约 9.94Hz。转频约29.5Hz。数据分成两类正常状态和内圈故障状态。每段信号取1024个点大概相当于0.085秒在轴承工频29.5Hz下能覆盖2.5个转周期足够反映一个周期的故障冲击。5.2 完整特征提取流程整个流程分成三步。第一步时域特征提取。对每段数据计算10个时域特征。重点关注的指标是RMS和峭度。正常状态RMS大约0.02到0.03内圈故障状态RMS可能到0.09到0.15峭度正常状态接近3故障状态往往在5到10之间。第二步频域分析。先对原始信号做FFT看频谱你会发现直接FFT的频谱中故障特征频率处的峰值并不明显因为故障冲击激发的共振频率通常很高几kHz到十几kHz但能量分散。这时候需要用包络分析先对信号做带通滤波比如8kHz到12kHz然后取希尔伯特变换的模得到包络信号再对包络信号做FFT故障特征频率就会在低频段清晰可见。这个技术叫包络解调分析是轴承诊断的黄金标准。包络分析的核心代码% 带通滤波 [b, a] butter(4, [8000 12000] / (fs/2), bandpass); x_filtered filtfilt(b, a, x); % 希尔伯特变换求包络 envelope abs(hilbert(x_filtered)); % 对包络做FFT N_env length(envelope); X_env fft(envelope); X_env_mag abs(X_env) / N_env; f_env (0:N_env-1) * fs / N_env; % 只看0-1000Hz段 idx f_env 0 f_env 1000; plot(f_env(idx), X_env_mag(idx));第三步提取频域特征。在包络谱上找理论故障频率附近的峰值幅值比如在157.94Hz附近找最大峰值这个幅值就是内圈故障的强度指标。同时计算重心频率、均方频率等全局特征。5.3 特征组合带来的分类效果把时域频域的特征合并成特征向量后我用随机森林做二分类正常/故障交叉验证准确率能达到99%以上。但如果只用时域特征准确率大约90%只用频域特征大约95%。两者的组合效果最好这就是时频特征结合的价值所在。这里顺便提一嘴MFCC特征。MFCC梅尔频率倒谱系数是语音识别里的经典特征流程是分帧→加窗→FFT→梅尔滤波器组→取对数→DCT→得到倒谱系数。本质上就是在频域特征上再做一次压缩变换。如果你做的是音频信号的特征提取MFCC基本是绕不开的方案但做机械振动诊断时MFCC用得少因为梅尔滤波器组的频率刻度是为听觉设计的不适合机械故障频率分析。这个例子说明特征提取方法必须和应用场景匹配不能盲目套用。6. 常见问题与排查技巧实录6.1 时域特征看起来没区别怎么办这是最常遇到的问题正常和故障状态下提取出来的特征值差异极小完全区分不开。我遇到这种问题第一反应是看数据本身有没有问题。先用时域波形画出来看看是否存在削波传感器量程不够导致波形顶部被截平、断线、非线性漂移。削波会让峰值失真峰值因子和峭度都会异常。数据没问题的话考虑是不是特征选的阶段不对。早期微弱故障RMS和峰值都变化很小但高频段能量和峭度可能已经有反应。这时不要死盯着时域特征直接上包络谱分析。还有一种情况是采样长度不够。特征提取对样本长度有要求特别是峭度这种高阶统计量数据太短比如只有几十个点估计值方差巨大根本不可靠。我建议单段特征提取的样本长度至少覆盖10个转频周期做轴承诊断的话每段至少1024个点最好4096个点。6.2 频谱峰值偏移对不上理论值故障特征频率应该在某一个频率处出现峰值实测发现峰值偏了这可能有两种情况。第一种转速波动。理论故障频率是假设转速恒定算出来的实际设备转速可能有波动负载变化、电源波动导致故障频率偏移。解决方法是同时采集转速信号用转速信号做等角度重采样角域平均或者把故障频率按照实际转速归一化到转频倍数。第二种频率分辨率不够。比如ls12kHzN1024频率分辨率就是fs/N ≈ 11.7Hz而内圈故障频率和滚动体故障频率的差可能不到10Hz频谱上两个峰值合成一个根本分辨不出来。解决方法是增加样本长度N分辨率变细但代价是计算时间变长。顺带提一下“方位向压缩是压缩频域还是时域”这个问题。这是合成孔径雷达SAR成像里的概念方位向压缩确实是在频域多普勒频域做的处理——通过匹配滤波对方位向信号进行压缩提高方位向分辨率。本质上和我们在信号处理里做FFT后滤波再IFFT回到时域是一回事。这说明了频域处理在很多领域都是绕不开的核心操作。6.3 FFT结果幅值来回跳不稳定同一段信号重复做FFT幅值应该是一样的因为FFT是确定性计算。如果你发现幅值不稳定大概率是信号本身不稳定不是FFT的问题。一个常见原因是数据采集的触发条件不一致导致每帧数据的起始相位不同。对非相参信号这会导致FFT幅值出现波动虽然功率谱密度会稳定一些。解决方案是用重叠率更高的窗函数增大overlap或者用Welch方法的平均功率谱来替代单帧FFT。另外检查一下是否有瞬态干扰混入。比如轴承故障信号的冲击性很强如果某帧数据恰好包含了一个大冲击RMS和峰值都会偏高。处理方法是做中值滤波或者异常值剔除。6.4 特征矩阵里出现NaN怎么办做特征批量提取时如果某个特征出现了NaN通常原因是分母为0。比如波形因子是RMS除以绝对平均幅值如果某帧数据全为0传感器没输出分母为0结果就是NaN或者Inf。这个一定要在代码里加保护if abs(mean(abs(x))) eps shape_val 0; % 或者标记为异常样本 else shape_val rms_val / mean(abs(x)); endMatlab里还可以用isnan检查最终特征矩阵把含NaN的行剔除掉。但要注意如果剔除的比例很高比如超过5%说明数据质量有问题先排查传感器和数据采集系统而不是简单粗暴地删除样本。6.5 时域掩蔽效应和滚动时域优化是什么网络热词里提到了时域掩蔽效应和滚动时域优化这两个词容易混在一起但完全是两个领域的概念。时域掩蔽效应是心理声学听觉感知里的概念当一个强的声音出现后它会对后面紧接着出现的弱声音产生掩蔽导致弱声音在一段时间内听不见。这在音频编码比如MP3、AAC中被利用通过心理声学模型把听不见的信息丢掉来压缩数据。如果你做音频处理这个效应会影响特征提取——有些频谱细节虽然真实存在但人耳听不见在特征中是否能被利用取决于你是做人耳感知相关的任务还是机器诊断任务。滚动时域优化Receding Horizon Optimization是控制理论里的概念核心思路是每个控制周期都基于当前状态做一个有限时域的最优控制问题执行第一步然后滚动到下一个周期重新优化。这个名字里虽然有“时域”两个字但它和时域特征提取完全不是一个层面的事。这两个词提醒我们搜索时域相关的内容时要注意区分“时域分析”在信号处理和在其他领域的含义避免概念混用。6.6 特征提取的速度优化如果特征是用于在线监测的性能优化很重要。讲几个实测有效的优化技巧。技巧一是避免在循环里逐点操作。用numpy或Matlab的向量化计算比在循环里逐点加要快几十倍。写特征提取脚本先考虑有没有现成的向量化函数比如rms、kurtosis、pwelch这些内部都是高度优化的比自己手写循环快得多。技巧二是合理选择FFT长度。不是越长越好长FFT计算慢而且对实时性要求高时没必要。如果只需要关心某个频段的特征可以先带通滤波再降采样然后再做FFT计算量可以大幅降低。技巧三是用缓存复用中间结果。比如特征提取里RMS、峰值、峭度都要用到信号幅值如果对数据做了归一化处理可以先缓存归一化后的数据避免重复计算。技巧四是用流式计算处理时域特征。均值、RMS这些统计量可以写成递推形式每次进来一个新数据点就更新一次不需要缓存整段数据。对于嵌入式设备这种流式特征提取方式几乎是唯一选择。7. 特征提取的工程落地经验分享7.1 一套完整的特征提取代码框架最后分享一个我实际在项目里用的特征提取框架把时域和频域特征合并成一个函数输入一段信号输出完整的特征向量。function feats extract_features(x, fs) % 输入: x - 一段信号, fs - 采样率 % 输出: feats - 特征向量依次为时域10个、频域9个 % 时域特征 N length(x); mean_val mean(x); rms_val rms(x); peak_val max(abs(x)); p2p_val peak2peak(x); crest_val peak_val / rms_val; kurtosis_val kurtosis(x); skewness_val skewness(x); shape_val rms_val / (sum(abs(x)) / N); impulse_val peak_val / (sum(abs(x)) / N); margin_val peak_val / (mean(sqrt(abs(x))))^2; feats_time [mean_val, rms_val, peak_val, p2p_val, ... crest_val, kurtosis_val, skewness_val, ... shape_val, impulse_val, margin_val]; % 频域特征 % 去均值避免直流干扰 x x - mean(x); % 功率谱估计用汉宁窗50%重叠 window hann(256); noverlap 128; nfft 512; [pxx, f] pwelch(x, window, noverlap, nfft, fs); % 只分析正频率范围 pos_idx f 0; f f(pos_idx); pxx pxx(pos_idx); K length(pxx); FC sum(f .* pxx) / sum(pxx); MSF sum(f.^2 .* pxx) / sum(pxx); VF MSF - FC^2; [max_pxx, max_idx] max(pxx); peak_freq f(max_idx); % 频带能量比这里用的6个等宽频带 band_ratios zeros(1, 6); band_edges linspace(0, max(f), 7); E_total sum(pxx); for k 1:6 band_mask f band_edges(k) f band_edges(k1); band_ratios(k) sum(pxx(band_mask)) / E_total; end feats_freq [FC, MSF, VF, peak_freq, max_pxx, band_ratios]; % 合并 feats [feats_time, feats_freq]; end这个函数可以直接用来批量处理数据。需要注意一点pwelch的窗口长度、nfft这些参数要根据你的采样率和目的频率范围调整如果故障特征频率在100Hz附近窗口别设太短否则频率分辨率不够pwelch会把低频细节抹掉。我这里的256窗口在12kHz采样率下对应约47Hz分辨率对轴承故障分析是够用的但如果做齿轮箱分析特征频率通常在1kHz以上窗口可以设更短来换取时间分辨率。7.2 特征提取和机器学习的衔接特征提取得再好最终还是要喂给分类器或者回归模型。这里有几个实操建议。特征矩阵的格式每一行是一个样本每一列是一个特征。拿正常样本30组、故障样本30组提取特征后得到60行×19列的特征矩阵加上标签列一共60×20。特征标准化很多机器学习算法对特征的尺度敏感比如SVM和KNN。建议在送入模型前做标准化在Python里用StandardScaler在Matlab里用zscore。标准化系数均值和标准差必须只用训练集计算然后应用到测试集否则会造成数据泄漏导致模型评估结果虚高。特征选择如果特征维度高可以用随机森林的特征重要性排序、递归特征消除RFE等方法筛掉不重要的特征。但对于故障诊断这种场景我更推荐保留结果可解释性强的特征哪怕模型准确率稍微低一点客户能看懂为什么报警价值更大。模型选择样本量不大几十到几百组的情况下随机森林和梯度提升树是默认选择不需要神经网络也能达到很好的效果。样本量上几千组之后可以尝试一维CNN直接从原始信号提取特征这就进入端到端学习的范畴了但可解释性会大幅下降。7.3 最后再提醒几个实践中的经验根据我的经验时域频域特征提取这个领域最核心的能力不是会套公式而是建立“信号→特征→状态”的完整链路。我见过太多人卡在中间某一步——或者FFT的幅值没搞清楚或者特征选了一堆但没法解释导致整个项目交付不了。几个最后的补充提醒关于采样率特征提取之前先确认采样率是否满足你的分析需求。采样定理要求采样率至少大于最高关注频率的两倍实际工程建议留3到5倍裕量。做轴承故障诊断如果故障特征频率在10kHz附近采样率至少要20kHz最好40kHz以上。关于样本长度时域统计特征需要足够长的数据才能稳定估计。峭度这种高阶统计量特别依赖于样本长度。经验法则峭度至少要1024个点才可靠最好4096以上。关于参考标准很多行业都有关于振动特征的评判标准比如ISO 10816用RMS做设备状态分级ISO 7919用轴振动位移做评判。做工程诊断时先查一下你所在行业有没有对应标准而不是自己定一个阈值。关于数据记录做特征提取前把传感器类型、采样率、量程、安装位置、设备转速、负载情况全部记录下来。同样的特征值在不同工况下含义完全不同。我做过一个项目现场工程师把传感器装反了方向导致采集的所有数据极性反转特征值全乱了——这种问题单纯靠算法永远发现不了必须有完整的现场记录来辅助判断。本文还有配套的精品资源点击获取