拿到一段振动信号、声音信号或者任意一维时间序列之后第一步该做什么很多刚接触故障诊断和状态监测的同学都会愣住。采集仪给你一堆波形你盯着它知道这里面有信息但不知道怎样才能把它变成“正常的”“磨损的”“点蚀的”这种可以拿去训练的标签。我的习惯很直接先做特征提取一次性把时域、频域、时频域三个维度的统计指标全部算出来再用这些特征去观察数据、训练模型。这篇文章就把我在Matlab里常用的36种多征兆域特征提取方法完整梳理一遍包含可直接运行的代码、参数设置思路和踩过的坑。1. 36种特征指标的分域地图时域、频域、时频域各在解决什么问题很多教材会把特征提取封装成“提取特征”三个字但实际做项目时你会发现不同域的特征看的是信号的不同侧面。我在实际项目中习惯把36种特征分成三大类16个时域统计特征、13个频域谱特征、7个时频域特征。这样划分不是为了凑数是因为诊断一个故障往往需要多个域的信息互相印证。1.1 时域特征直接观察波形形态和能量分布时域特征直接对原始时间序列做统计计算它反映的是信号幅度、波动程度和分布形态。有量纲的特征比如均值、方差、标准差、均方根值、峰值、峰峰值、方根幅值、绝对平均值、能量、均方值这些跟信号的绝对幅度有关能反映信号的整体强度。无量纲的特征比如峭度、偏度、波形指标、峰值指标、脉冲指标、裕度指标这些剔除了幅度量纲的影响更多反映波形的形状特征。峭度是我做轴承故障诊断时最先看的指标。健康轴承振动信号的峭度接近3一旦出现早期点蚀或剥落冲击成分会让峭度迅速飙升到4.5以上甚至到6-8。波形指标和峰值指标在故障发展到中后期时变化更稳定可以作为峭度的补充。偏度则能提示波形的不对称性比如不对中或摩擦产生的削波往往让偏度发生变化。1.2 频域特征从频率成分分布中寻找故障指纹频域特征是对信号做傅里叶变换之后对频谱进行统计。故障信号不等于纯随机噪声它会在特定频率上出现能量聚集。比如齿轮啮合频率及其边频带、轴承的外圈或内圈故障特征频率这些都是故障的“指纹信息”。频域特征包括重心频率、均方频率、频率方差、均方根频率、谱幅值均值、谱幅值标准差、谱峰幅值、谱峰频率、频带总能量、谱熵、平均功率频率、边频带幅值均值、边频带幅值比这13个特征。谱熵反映频谱能量分布的均匀程度健康信号的谱熵通常较高因为能量分布相对分散而带故障的信号频谱能量集中在特定频率附近谱熵会下降。1.3 时频域特征应对非平稳信号的关键手段变速运行、启停机过程、冲击载荷作用下的信号频率成分随时间变化单独的时域或频域分析都不够。时频域特征需要同时保留时间和频率信息。我这里的7个时频域特征包括小波包总能量、小波包能量熵、特定频带能量占比、边际谱重心频率、瞬时频率均值、瞬时幅值均值、瞬时幅值标准差。小波包分解可以把信号按频带细分比普通小波更精细。瞬时频率和瞬时幅值通过希尔伯特变换提取能捕捉信号每个时刻的振荡快慢和强弱。这几个特征对变转速工况特别有效。特征类别数量代表指标主要解决的问题时域特征16个峭度、RMS、峰值指标、脉冲指标冲击、磨损、整体能量变化频域特征13个重心频率、谱熵、边频带幅值比频率成分变化、调制故障时频域特征7个小波包能量熵、瞬时频率均值非平稳信号、变速工况2. 为什么必须组合三个域才能完成可靠诊断一个轴承故障实例的完整推理只有时域特征很容易误判。我之前遇到过一组实测的轴承振动数据正常状态下RMS本身就波动很大仅靠RMS根本分不出正常和点蚀早期。原因是RMS同时受转速、负载、传感器灵敏度影响工况一变健康状态的RMS可能比故障状态还高。把频域特征加进来后情况立刻不同。正常轴承的振动频谱相对平坦故障轴承会在某一特征频率附近出现明显的谱峰。齿轮箱磨损时啮合频率的边频带幅值比会变大。这就是频域特征的价值—它能定位“故障发生在哪个频率上”。但单靠频域还不够。内圈故障特征频率附近会出现调制边带但边带的幅度会随转速波动而变化如果只看某一段频谱可能错过故障特征。这时时频域特征派上用场。小波包能量熵能把信号能量在各个频带中的分布情况提炼成一个标量故障发生后能量会向高频带转移能量熵数值随之变化。三个域的指标组合起来诊断逻辑就完整了时域确认冲击和能量突变频域锁定故障特征频率时频域验证非平稳特性和频带转移。我在实际项目里给出诊断结论时通常要求至少两个域的特征都指向同一类故障才会下判断。比如峭度升高同时对应谱峰频率出现在外圈特征频率处这样的组合结论基本跑不了。再补充一个数值经验。计算边频带特征时边频带的带宽不能死板地取固定值我会根据谱峰附近的能量分布自适应选取默认取谱峰频率左右3%作为边频带范围。如果已知故障特征频率比如轴承节径和接触角能算出来就把中心频率替换成理论值此时边频带特征更准确。3. Matlab环境下的数据清洗与预处理特征计算前的最后一道防线特征提取的结果再漂亮输入数据一塌糊涂也白搭。我在处理信号之前有一套固定的预处理流程不到二十行代码就能完成但能挡住大部分数据异常问题。3.1 输入数据的统一规范化Matlab里容易忽略的一件事是列向量和行向量的差异。很多采集软件导出的数据是一行有些是转置后的一列。我在特征提取函数开头会先强制把输入变成列向量function x ensure_column(x) % 强制转为列向量避免size不匹配问题 if isrow(x) x x.; end end这个细节看起来很基础但实际踩坑最多的就是这个。尤其是批量处理几十个文件时有的文件是行向量有的文件是列向量后面写循环索引时非常容易出错。3.2 NaN、Inf和直流分量处理直接对包含NaN的信号算均值和标准差值全是NaN后面特征筛选和分类模型全崩。我会先定位并替换非有限值。工程上最简单的做法是把NaN和Inf替换为0但如果一整段都是NaN替换成0反而会引入一段虚假的低幅信号。更稳妥的做法是在预处理前就检查异常段如果异常点占比超过5%直接删除整段数据。去趋势用的是Matlab自带的detrend函数。它默认把信号中的线性趋势去掉相当于移除传感器漂移和温度漂移带来的低频分量。做特征提取前不去趋势时域均值会偏移频域低频部分会被虚假抬高。function x_clean preprocess_signal(x) x ensure_column(x); x detrend(x, linear); % 替换非有限值 x(~isfinite(x)) 0; end3.3 分段策略窗口长度与重叠率的选择特征提取不是拿整段信号算一次就行那样只有一个样本根本不够训练模型。实际项目里要把长信号切成多个窗口样本。窗口长度怎么选我遵循一个经验公式窗口长度取信号主要周期的5到10倍。比如齿轮轴的转频是10Hz一个周期是0.1秒那么窗口至少取0.5到1秒。窗口太短频率分辨率不足频域特征会失真窗口太长时间分辨率下降会丢失瞬态信息。重叠率我一般取50%这样既增加了样本数量又不会让相邻窗口高度重复。重叠率超过75%时样本之间高度相关模型容易过拟合这个经验值不要轻易突破。3.4 滑动窗口批量生成样本的代码框架function segments segment_signal(x, win_len, overlap) step round(win_len * (1 - overlap)); n_seg floor((length(x) - win_len) / step) 1; segments zeros(n_seg, win_len); for i 1:n_seg start_idx (i - 1) * step 1; segments(i, :) x(start_idx : start_idx win_len - 1); end end拿到二维矩阵后每一行就是一个样本配合后续的特征提取函数一次性对所有窗口批量计算36个特征输出一张样本数乘以36的特征表。这比在循环里一个个重复写特征计算代码要省事得多。4. 核心代码36种特征批量提取函数的可直接落地版本我确定的36种特征时域16个、频域13个、时频域7个。下面这段主函数直接返回一个1×36的特征向量。代码里我用了表驱动的方式方便你后续修改或裁减。4.1 时域特征计算代码时域特征计算全部建立在基础统计量之上。要注意的是波形指标、峰值指标、脉冲指标、裕度指标的分子分母都涉及峰值或绝对均值当信号非常平稳且接近纯正弦时部分指标可能异常偏大这是正常现象。计算时我习惯加上一个极小值防除零。function tfeats extract_time_features(x) % 输入: 信号列向量 x % 输出: 16个时域特征的结构体 x ensure_column(x); N length(x); % 基础统计量 mx mean(x); vx var(x); sx std(x); armx rms(x); px max(abs(x)); ppx max(x) - min(x); ra mean(sqrt(abs(x)))^2; % 方根幅值 aav mean(abs(x)); % 绝对平均值 ku kurtosis(x); % 峭度 sk skewness(x); % 偏度 % 无量纲指标 wf armx / (aav eps); % 波形指标 cf px / (armx eps); % 峰值指标 if_ px / (aav eps); % 脉冲指标 clf px / (ra eps); % 裕度指标 % 能量类指标 en sum(x.^2); % 信号能量 smq mean(x.^2); % 均方值 tfeats struct(mean, mx, var, vx, std, sx, rms, armx, ... peak, px, ptp, ppx, root_amp, ra, aav, aav, ... kurt, ku, skew, sk, waveform, wf, crest, cf, ... impulse, if_, clearance, clf, energy, en, square_mean, smq); end峭度、偏度、标准差和均方根是最常用的四个。峭度对早期冲击敏感偏度能反映非对称性标准差和均方根反映整体波动。四个组合使用时我遇到过正常但偏度显著变化的情况排查后确认是传感器接线接触不良导致的单边削波——这类问题不做特征提取根本发现不了。4.2 频域特征计算代码频域特征先做FFT取单边频谱然后对频谱本身做统计。谱幅值均值、谱幅值标准差衡量的是频谱整体的平坦程度重心频率和均方根频率衡量频谱能量集中在高频还是低频谱熵则从信息量角度衡量频谱的分散程度。这段代码里有一个容易踩的坑FFT之后频谱是复数很多人直接用abs()之后不除以N结果谱幅值随信号长度变化导致同一个信号在不同长度下特征值不稳定。我统一除以N做幅值归一化。function ffeats extract_freq_features(x, fs) x ensure_column(x); N length(x); % 单边幅值谱归一化 Y abs(fft(x)) / N; P Y(1:floor(N/2)1); P(1) P(1) / 2; % 直流分量修正 f (0:floor(N/2)) * fs / N; % 谱统计量 sumP sum(P); fc sum(f .* P) / (sumP eps); % 重心频率 msf sum(f.^2 .* P) / (sumP eps); % 均方频率 vf sum((f - fc).^2 .* P) / (sumP eps); % 频率方差 rmsf sqrt(msf); % 均方根频率 smean mean(P); % 谱幅值均值 sstd std(P); % 谱幅值标准差 [maxP, maxIdx] max(P); peakf f(maxIdx); % 谱峰频率 bandE sum(P.^2); % 频带总能量(幅值平方和) % 谱熵: 概率分布形式 pP P / (sumP eps); pP(pP 0) 1e-12; specEntropy -sum(pP .* log(pP)); % 谱熵 % 平均功率频率 mpf sum(f .* P.^2) / (sum(P.^2) eps); % 边频带特征: 以谱峰为中心, 取±3%频率范围 df fs / N; bandWidth max(round(0.03 * peakf / df), 2); loIdx max(maxIdx - bandWidth, 1); hiIdx min(maxIdx bandWidth, length(P)); sidebandMean mean(P(loIdx:hiIdx)); % 边频带幅值均值 sidebandRatio sidebandMean / (sumP eps); % 边频带幅值比 ffeats struct(fc, fc, msf, msf, vf, vf, rmsf, rmsf, ... smean, smean, sstd, sstd, maxP, maxP, peakf, peakf, ... bandE, bandE, specEntropy, specEntropy, mpf, mpf, ... sidebandMean, sidebandMean, sidebandRatio, sidebandRatio); end如果你处理的是变转速信号做频谱之前建议先做角域重采样订单跟踪否则转频变化会把频谱峰值抹平边频带特征几乎失效。4.3 时频域特征计算代码时频域特征用到了小波包和希尔伯特变换。使用wprcoef获取小波包分解后各节点的重构系数再计算节点能量。db4小波是处理振动信号时比较折中的选择db2太短频率分辨率差db8以上计算量变大且易受噪声影响。三层小波包把频带分成8个节点对大多数机械信号足够。function tfteats extract_tf_features(x, fs) x ensure_column(x); N length(x); % 三层小波包分解 wpt wpdec(x, 3, db4); nodes [3 0; 3 1; 3 2; 3 3; 3 4; 3 5; 3 6; 3 7]; nodeE zeros(8, 1); for i 1:8 coefs wprcoef(wpt, nodes(i, :)); nodeE(i) sum(coefs.^2); end totalE sum(nodeE); % 小波包总能量 pNode nodeE / (totalE eps); % 各节点能量占比 % 小波包能量熵 pNode(pNode 0) 1e-12; wpEntropy -sum(pNode .* log(pNode)); % 边际谱重心频率(用小波包节点中心频率加权近似) deltaF fs / (2 * 8); % 每个节点的频带宽度 centers (nodes(:, 1) - 3) * 0 ... % 中心频率按节点序号和层数计算 (nodes(:, 2) 0.5) * deltaF; % 简化处理: 按节点索引编号0-7近似中心频率 centers (nodes(:, 2) 0.5) * deltaF; marginalFC sum(centers .* pNode); % 希尔伯特变换提取瞬时特征 z hilbert(x); instAmp abs(z); instFreq fs / (2 * pi) * diff(unwrap(angle(z))); instFreqMean mean(instFreq); tfteats struct(wpTotal, totalE, wpEntropy, wpEntropy, ... bandEnergyRatio, pNode(1), marginalFC, marginalFC, ... instFreqMean, instFreqMean, instAmpMean, mean(instAmp), ... instAmpStd, std(instAmp)); end这个函数里我为了保持输出是7个特征bandEnergyRatio取了第一个节点即最低频带的能量占比。如果你需要更细的频带区分度可以把pNode里其他几个也作为单独特征加到特征表里特征总数会超过36个这种扩展完全可行。4.4 主函数串联一次输出一张完整的特征表三个子函数写完后主函数负责调用并将所有特征拼接成一个1×36的行向量。为了让结果便于后续使用我建议输出Matlab table而不是普通数组因为列名能直接保留特征名称做特征筛选和可视化时一眼就能认出来。function featRow extract_all_features35(x, fs) % 完整36种特征提取函数 tX extract_time_features(x); fX extract_freq_features(x, fs); tfX extract_tf_features(x, fs); % 展平结构体为数值行向量 f1 struct2array(tX); f2 struct2array(fX); f3 struct2array(tfX); featRow [f1, f2, f3]; end批量处理时对segments矩阵的每一行调用这个函数就行。耗时取决于小波包分解和窗口数量10万点以下的信号单窗口在一秒内能算完。如果窗口数量特别大配合parfor并行循环可以明显提速。5. 特征筛选36个指标怎么浓缩成真正有区分度的特征集特征不是越多越好。36个特征之间往往存在高度相关性比如波形指标和峰值指标经常一起升高谱幅值均值和频带总能量也强相关。直接把36维特征扔给分类器不仅计算量变大还可能因为多重共线性让模型过拟合。我一般在特征提取之后做一道筛选。5.1 相关性矩阵找冗余先对全部样本的特征矩阵做corrcoef找出相关系数绝对值超过0.95的特征对。比如我多次看到波形指标和峰值指标的相关系数超过0.98此时保留其中一个即可。筛选时我倾向于保留物理意义更明确的那个。峰值指标对冲击更敏感波形指标则更偏波形形态两者在早期故障时表现差异不大后期故障时峰值指标更稳。5.2 对故障敏感性的经验排序针对机械设备诊断场景我总结了一个经验排序早期微弱故障峭度、脉冲指标、峰值指标最敏感其次是谱熵和小波包能量熵。中晚期故障边频带幅值比、重心频率、频带能量占比更可靠。齿轮类故障边频带幅值比、谱熵、边际谱重心频率优先。不平衡与不对中2倍频附近的谱能量特征、峰值指标优先。这个排序来自我实际做轴承、齿轮箱、电机诊断项目的统计结果不一定适用于所有场景但可以作为初筛时的起点。拿到新问题的数据后我会先用箱线图或者t-SNE看每个特征在正常与故障两类样本上的分布差异响应不够大的直接丢弃。5.3 工况-特征耦合陷阱有一个坑值得单独提醒。时域特征尤其是RMS、峰值、能量对转速和负载变化非常敏感。同样是正常状态800转时和1500转时的RMS差距明显而频域特征中的重心频率和边频带幅值比对工况变化的敏感度要低一些。如果你后续要用模型做跨工况诊断建议优先保留频域和时频域特征或者对时域特征单独做工况归一化比如除以当前转速下的均值。5.4 降维方法怎么选相关性阈值筛选之后如果特征维度仍然偏高可以用PCA再做一次压缩但我会控制特征数量不要过少。PCA得到的主成分物理意义不明确不利于后期诊断解释。实际项目里我通常先用相关性筛选加敏感度排序把特征压到8到12个如果分类精度仍然很理想就不再降维了。只有特征维度还是很高且样本量有限时才会考虑PCA或LDA。这样既保留了可解释性又控制了维度。6. 调试经验与常见问题实战记录前面讲的都是正常流程但真实操作中总会冒出各种问题。我把遇到过的几个高频问题列出来这些经验大多是用实际项目教训换来的。6.1 工具箱依赖检查kurtosis和skewness来自Statistics and Machine Learning Toolboxwpdec和wprcoef来自Wavelet Toolboxhilbert来自Signal Processing Toolbox。如果你用的是精简版Matlab环境可能这些函数缺失。我习惯在特征提取脚本开头加一段工具箱检查一旦缺失立即提示而不是等跑到那一行才报错assert(~isempty(ver(signal)), 缺少信号处理工具箱); assert(~isempty(ver(wavelet)), 缺少小波工具箱); assert(~isempty(ver(stats)), 缺少统计工具箱);6.2 内存不足与小波包计算缓慢小波包分解对信号长度敏感。一次处理100万点信号时wpdec会生成较大的树结构内存占用明显上升。我的处理方式有两种一是分段计算特征每个窗口单独做小波包分解避免一次性对整个长信号分解二是如果窗口数量很多改用parfor并行但也注意不能把parfor嵌套在大量结构体操作里否则通信开销反而更大。实测下来单机8核并行处理数十万窗口能提速3到5倍。6.3 频谱泄露导致边频带特征失真FFT之前不加窗函数时频谱会出现明显的泄漏现象谱峰周围假边带成分会干扰边频带幅值比计算。我解决后的问题是在做FFT之前加Hanning窗。这里有个小细节加窗会改变信号总能量所以计算频带总能量时需要对窗函数主瓣增益做补偿。简单做法是计算能量前除以窗函数的均方根值RMS这样加窗后能量近似不变。6.4 特征值出现Inf或NaN的常见来源除了原始信号里的NaN之外特征计算中Inf通常来自除零。波形指标、峰值指标、脉冲指标、裕度指标这四兄弟都有除零隐患当信号的绝对平均值或方根幅值非常小时计算结果会爆炸。我在代码里用eps做分母下限同时建议在特征计算完成后加一段数值检查把非有限值替换为该列的中位数防止个别异常样本破坏整个训练流程featTable array2table(featMat); for j 1:width(featTable) col featTable{:, j}; if any(~isfinite(col)) col(~isfinite(col)) nanmedian(col, omitnan); featTable{:, j} col; end end7. 经验收尾先全量计算再逐步筛选比一开始拍脑袋选特征更靠谱最后说一点我个人坚持了很久的工作习惯。拿到任何一段新的信号数据我都会先把36个特征全部批量算一遍存成一张完整的特征表格然后再做敏感度分析和筛选而不会一开始就凭经验只选三个特征。原因很简单指标之间的相关性经常反直觉。你以为峭度和峰值指标是一回事但实际数据里它们的相关系数可能只有0.6你以为谱熵和能量占比高度重合真正算出来却发现它们在区分早期故障时各有各的作用。只有把全量特征铺开看过、筛过才知道这个具体问题上哪些指标真正有区分度。这个习惯帮我避免了很多次重复返工。特征表一旦建立后续不管是训练分类模型、做聚类还是输出诊断报告都能直接复用。信号特征指标诊断这件事说到底是先让数据变得可比较、可解释坐上这个环节的功夫后面的模型和判断才有根基。