简介：本资源是一份面向信号处理初学者与工程实践者的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_{i=1}^{N}x_i ]
代表信号的直流分量。注意一个坑：很多传感器出来的信号做完隔直处理后，均值理论上接近0，但如果均值明显漂移，说明传感器零漂或者存在低频趋势项，这对后续频域分析影响很大。
均方根值（RMS）：
[ x_{rms} = \sqrt{\frac{1}{N}\sum_{i=1}^{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_{i=1}^{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_{i=1}^{N}\sqrt{|x_i|}\right)^2 ]
Matlab里用mean(sqrt(abs(x)))再平方就对了。这个特征在实际轴承诊断里对早期故障很敏感，别漏掉。
如果用的是Python，核心逻辑一模一样，换成numpy就行。需要注意的是，这两种语言里kurtosis的定义有差异，Matlab的kurtosis默认返回的是超值峭度（即减去3），Python的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的话用numpy：
import numpy as np X = np.fft.fft(x) f = np.fft.fftfreq(N, d=1/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);功率谱的好处是抗噪能力强、平滑度好，特别适合提取宽带背景下的窄带特征。在实际故障诊断中，我更倾向于用功率谱，因为故障特征频率往往表现为功率谱上的一个凸起，比幅值谱更明显。
顺带提一个热点词里出现过的概念：频域OCT（Optical Coherence Tomography，光学相干层析）的频谱域实现。医学成像领域把OCT信号从时域（扫描深度方向）变换到频域，通过分析干涉光谱的频域信息来重构样品内部结构。核心思路和我们做振动信号频域分析完全一样，都是利用傅里叶变换把信号从时域映射到频域，只是在那个场景下频域编码了空间位置信息。这就是频域分析的魅力，同一个数学工具在不同领域能玩出完全不同的花样。
3.3 从频域里能提取哪些特征
频谱不是一个特征，而是一堆特征。我整理一下我实际用过的频域特征，按用途分类。
频率中心（重心频率，FC）：
[ FC = \frac{\sum_{i=1}^{K}f_i \cdot P_i}{\sum_{i=1}^{K}P_i} ]
其中P_i是第i个频率点的功率（或幅值），f_i是该频率点的频率。重心频率反映信号能量在频谱上的集中位置。轴承正常情况下能量集中在低频段，故障时往往会在高频段出现新的能量聚集，重心频率就会往后移。
均方频率（MSF）和频率方差（VF）：
[ MSF = \frac{\sum_{i=1}^{K}f_i^2 \cdot P_i}{\sum_{i=1}^{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精度更高，但要小心过冲问题。
另一个思路是用nufft（Non-Uniform FFT），Matlab从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 频谱峰值偏移，对不上理论值
故障特征频率应该在某一个频率处出现峰值，实测发现峰值偏了，这可能有两种情况。
第一种，转速波动。理论故障频率是假设转速恒定算出来的，实际设备转速可能有波动（负载变化、电源波动），导致故障频率偏移。解决方法是同时采集转速信号，用转速信号做等角度重采样（角域平均），或者把故障频率按照实际转速归一化到转频倍数。
第二种，频率分辨率不够。比如ls=12kHz，N=1024，频率分辨率就是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(k+1); 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用轴振动位移做评判。做工程诊断时，先查一下你所在行业有没有对应标准，而不是自己定一个阈值。
关于数据记录：做特征提取前，把传感器类型、采样率、量程、安装位置、设备转速、负载情况全部记录下来。同样的特征值在不同工况下含义完全不同。我做过一个项目，现场工程师把传感器装反了方向，导致采集的所有数据极性反转，特征值全乱了——这种问题单纯靠算法永远发现不了，必须有完整的现场记录来辅助判断。
本文还有配套的精品资源，点击获取
