企业官网建设流程全解析

拿到一段振动信号、声音信号或者任意一维时间序列之后，第一步该做什么？很多刚接触故障诊断和状态监测的同学都会愣住。采集仪给你一堆波形，你盯着它，知道这里面有信息，但不知道怎样才能把它变成“正常的”“磨损的”“点蚀的”这种可以拿去训练的标签。我的习惯很直接：先做特征提取，一次性把时域、频域、时频域三个维度的统计指标全部算出来，再用这些特征去观察数据、训练模型。这篇文章就把我在Matlab里常用的36种多征兆域特征提取方法完整梳理一遍，包含可直接运行的代码、参数设置思路和踩过的坑。

1. 36种特征指标的分域地图：时域、频域、时频域各在解决什么问题

很多教材会把特征提取封装成“提取特征”三个字，但实际做项目时你会发现，不同域的特征看的是信号的不同侧面。我在实际项目中习惯把36种特征分成三大类：16个时域统计特征、13个频域谱特征、7个时频域特征。这样划分不是为了凑数，是因为诊断一个故障往往需要多个域的信息互相印证。

1.1 时域特征：直接观察波形形态和能量分布

时域特征直接对原始时间序列做统计计算，它反映的是信号幅度、波动程度和分布形态。有量纲的特征，比如均值、方差、标准差、均方根值、峰值、峰峰值、方根幅值、绝对平均值、能量、均方值，这些跟信号的绝对幅度有关，能反映信号的整体强度。无量纲的特征，比如峭度、偏度、波形指标、峰值指标、脉冲指标、裕度指标，这些剔除了幅度量纲的影响，更多反映波形的形状特征。

峭度是我做轴承故障诊断时最先看的指标。健康轴承振动信号的峭度接近3，一旦出现早期点蚀或剥落，冲击成分会让峭度迅速飙升到4.5以上，甚至到6-8。波形指标和峰值指标在故障发展到中后期时变化更稳定，可以作为峭度的补充。偏度则能提示波形的不对称性，比如不对中或摩擦产生的削波往往让偏度发生变化。

1.2 频域特征：从频率成分分布中寻找故障指纹

频域特征是对信号做傅里叶变换之后，对频谱进行统计。故障信号不等于纯随机噪声，它会在特定频率上出现能量聚集。比如齿轮啮合频率及其边频带、轴承的外圈或内圈故障特征频率，这些都是故障的“指纹信息”。

频域特征包括重心频率、均方频率、频率方差、均方根频率、谱幅值均值、谱幅值标准差、谱峰幅值、谱峰频率、频带总能量、谱熵、平均功率频率、边频带幅值均值、边频带幅值比这13个特征。谱熵反映频谱能量分布的均匀程度，健康信号的谱熵通常较高，因为能量分布相对分散；而带故障的信号频谱能量集中在特定频率附近，谱熵会下降。

1.3 时频域特征：应对非平稳信号的关键手段

变速运行、启停机过程、冲击载荷作用下的信号，频率成分随时间变化，单独的时域或频域分析都不够。时频域特征需要同时保留时间和频率信息。我这里的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; end

3.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 Toolbox，wpdec和wprcoef来自Wavelet Toolbox，hilbert来自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 end

7. 经验收尾：先全量计算再逐步筛选，比一开始拍脑袋选特征更靠谱

最后说一点我个人坚持了很久的工作习惯。拿到任何一段新的信号数据，我都会先把36个特征全部批量算一遍，存成一张完整的特征表格，然后再做敏感度分析和筛选，而不会一开始就凭经验只选三个特征。原因很简单：指标之间的相关性经常反直觉。你以为峭度和峰值指标是一回事，但实际数据里它们的相关系数可能只有0.6；你以为谱熵和能量占比高度重合，真正算出来却发现它们在区分早期故障时各有各的作用。只有把全量特征铺开看过、筛过，才知道这个具体问题上哪些指标真正有区分度。

这个习惯帮我避免了很多次重复返工。特征表一旦建立，后续不管是训练分类模型、做聚类还是输出诊断报告，都能直接复用。信号特征指标诊断这件事，说到底是先让数据变得可比较、可解释，坐上这个环节的功夫，后面的模型和判断才有根基。