简介面向信号处理学习者和MATLAB用户的时域频域特征提取代码包围绕从原始信号中提取统计特征与频谱特征的核心任务适用于机械故障诊断、语音识别等场景。压缩包共10个文件含6个mat数据文件和4个m脚本其中mat文件为示例信号或预处理后的特征矩阵m脚本则实现时域统计计算、频域统计计算及综合特征提取等流程整体体积20.68MB。已有2024人浏览学习。资料重点演示29维特征构建过程涵盖均值、峰峰值、均方根、峭度、过零率等时域指标以及频谱中心频率、带宽、功率谱密度等频域指标通过运行脚本可快速掌握特征提取从数据读取到结果输出的完整链路。相比孤立的理论公式这一代码包提供了可直接验证的MATLAB实现适合作为课程设计、项目预研或入门实践的参考。1. 从一段振动信号里找出设备故障matlab时域频域特征提取到底在提取什么做设备状态监测的工程师都有这种体验手里拿到一段加速度传感器数据几百兆的CSV波形看着乱糟糟FFT出来一堆峰值却说不清设备到底有没有故障。原因很简单原始信号是“黑匣子”你必须先从时域和频域两个角度把它压缩成一组有物理含义的数字这组数字就是特征。这套matlab时域频域特征提取资源包做的就是这件事把原始信号转成可供后续诊断或机器学习使用的特征向量。适合刚接触故障诊断、状态监测或振动分析的人也适合已经在跑深度学习但觉得“输入太长、模型太黑”想退回传统特征的工程师。它能用一段可运行的MATLAB代码帮你把峭度、峰值因子、重心频率、频带能量这些指标一次性算齐。2. 时域特征提取的代码骨架统计量只是起点形态指标才是关键2.1 时域特征为什么不能只靠mean和std很多人拿到信号第一时间算均值、方差、标准差然后发现这些特征对故障完全不敏感。原因不复杂均值反映的是信号中的直流分量标准差反映整体波动幅度但早期故障——比如轴承外圈点蚀——在波形上表现为周期性的窄脉冲它对均值和标准差的影响微乎其微却会让峭度和峰值因子发生剧烈变化。我在实际项目中一般把时域特征分成两组来看。第一组是有量纲特征包括均值、峰值、绝对平均值、方根幅值、标准差第二组是无量纲特征包括峭度、波形因子、峰值因子、脉冲因子、裕度因子。无量纲特征的好处是受工况变化影响小载荷变了、转速变了这些比值仍然比较稳定这是很多诊断标准选它们做报警参数的原因。资源包里的time_domain_features.m把这两组特征封装在一个函数里输入是一段信号向量输出是一个1×N的特征行向量。它的核心计算逻辑是先用abs取绝对值算绝对平均值用rms算有效值用max(abs(x))找峰值然后把它们两两相除得到无量纲指标。function features time_domain_features(x) % x: 单通道振动信号列向量 % features: 1x13 行向量依次为均值、峰值、绝对平均值、方根幅值、 % 标准差、有效值、峭度、波形因子、峰值因子、脉冲因子、裕度因子 x x(:); % 强制转为列向量 N length(x); mean_val mean(x); % 均值 peak_val max(abs(x)); % 峰值 abs_mean mean(abs(x)); % 绝对平均值 root_amp (mean(sqrt(abs(x))))^2; % 方根幅值 std_val std(x); % 标准差注意默认是N-1归一化 rms_val rms(x); % 有效值 kurtosis_val kurtosis(x); % 峭度MATLAB默认是三阶中心矩归一化 waveform_factor rms_val / abs_mean; % 波形因子 peak_factor peak_val / rms_val; % 峰值因子 pulse_factor peak_val / abs_mean; % 脉冲因子 margin_factor peak_val / root_amp; % 裕度因子 features [mean_val, peak_val, abs_mean, root_amp, std_val, ... rms_val, kurtosis_val, waveform_factor, peak_factor, ... pulse_factor, margin_factor]; end这段代码里最容易翻车的是向量的方向问题。如果x是行向量mean(x)得到的是一个行向量而不是标量最后拼出来的features维度会完全错乱。所以第一行x x(:)是必须的防御性写法。另一个细节是kurtosis的默认行为——MATLAB的kurtosis返回的是超额峭度也就是以3为基准做了平移。对于标准正态分布它返回0而不是3如果你后面要跟其他工具的数据对比一定要先搞清楚对方用的是不是同一套定义。2.2 特征排序和拼接给机器学习准备的标准输入格式单独算特征没有意义关键是要批量处理多段信号把每段信号的特征拼成一个大矩阵。资源包里提供了一个批处理脚本batch_extract_time.m核心思路是遍历一个文件夹下的所有mat文件对每个文件提取特征最后用vertcat拼接。% batch_extract_time.m % 批量提取时域特征输入文件夹路径输出特征矩阵和文件名单 folder D:\vibration_data\train; % 改成自己的数据路径 files dir(fullfile(folder, *.mat)); % 只读mat文件 num_files length(files); feature_matrix []; file_names {}; for i 1:num_files data load(fullfile(folder, files(i).name)); % 加载mat文件 % 假设mat文件里存的是变量 sig fea time_domain_features(data.sig); feature_matrix [feature_matrix; fea]; % 纵向拼接 file_names{end1} files(i).name; % 记录来源 end % 保存结果 save(time_features_result.mat, feature_matrix, file_names);这个脚本里我用feature_matrix [feature_matrix; fea]而不是预分配内存只适用于样本量不大于几百的情况。如果数据量大比如几千个mat文件最好改成预分配zeros(num_files, 13)再按行填入不然循环到后面会明显卡顿。但我也不建议一上来就做预分配先用这个简单版本跑通流程、确认特征维度正确再回头优化性能。data.sig这个变量名需要和你的mat文件内部变量名一致常常有人在这里踩坑——下载的数据集每个mat文件里变量名不一样有的叫vibration有的叫data。检查方式很简单MATLAB里用whos -file 文件名.mat就能查看变量名。2.3 时域特征怎么验证有效性特征算出来了怎么知道有用没用我常用的方法是做一次简单的“对比实验”取正常状态下的10段信号和故障状态下的10段信号分别提取特征然后看某个特征在两类样本之间的分布是否有明显分层。比如正常样本的峭度都在2.5到3.5之间故障样本的峭度在5到15之间那这个特征就基本可用。如果两个分布完全重叠直接丢弃。这一步不需要什么高级工具把特征矩阵按类别分组用histogram或boxplot画出来肉眼就能判断。这个习惯建议保留——我见过很多人辛辛苦苦算了几十个特征一股脑丢给分类器结果精度反而比只用两三个特征还差原因就是冗余特征太多把信噪比拉低了。3. 频域特征提取的两种路径FFT直接计算法和功率谱密度法3.1 从频谱到特征重心频率、均方根频率和频带能量频域特征的核心思想是把时域信号变换到频率轴然后从频谱形状里提取指标。最常见的三个指标是重心频率、均方根频率和频率标准差。重心频率表示频谱能量的重心位置设备出现故障时能量会向高频段移动重心频率随之变大。均方根频率和频率标准差则描述频谱分布的离散程度。资源包里的spectral_features.m用fft计算出单边幅值谱然后把频率轴和幅值谱代入公式计算。这里的关键不是fft本身——MATLAB里一行代码就能算——而是频率轴的构造和有效频段的选择。function [fc, rmsf, stdf, band_energy] spectral_features(x, fs) % x: 时域信号fs: 采样率 % fc: 重心频率, rmsf: 均方根频率, stdf: 频率标准差, band_energy: 频带能量 N length(x); X fft(x); % 复频谱 X X(1:floor(N/2)1); % 取单边谱 X_abs abs(X); % 幅值谱 freq (0:floor(N/2)) * fs / N; % 构造频率轴 fc sum(freq .* X_abs) / sum(X_abs); % 重心频率 rmsf sqrt(sum(freq.^2 .* X_abs) / sum(X_abs)); % 均方根频率 stdf sqrt(sum((freq - fc).^2 .* X_abs) / sum(X_abs)); % 频率标准差 % 频带能量把频谱分成10个等宽频带每个频带内求能量 num_bands 10; band_energy zeros(1, num_bands); max_freq fs / 2; band_width max_freq / num_bands; for k 1:num_bands idx freq (k-1)*band_width freq k*band_width; band_energy(k) sum(X_abs(idx).^2); end end这段代码有两个容易出错的点。第一fft的结果是复数取绝对值才是幅值谱但如果你要算能量应该用幅值的平方而不是幅值本身。我在band_energy里用的是X_abs(idx).^2这是功率的语义。第二单边谱长度是floor(N/2)1频率轴的间隔是fs/N如果N是奇数会有点绕但实际工程里采样点数基本都是2的幂次这个问题不常见。这个函数返回的band_energy是一个10维向量相当于把整个频谱等分成10个频段每个频段的能量构成一个子特征。这种做法在故障诊断里很实用——轴承故障的特征频率往往集中在某个特定频段你把全谱压缩成10个数字既降低了维度又保住了“哪个频段能量高”这个关键信息。3.2 功率谱密度估计为什么推荐pwelch而不是直接fft平方很多人算功率谱就是abs(fft(x)).^2 / N这在理论上是周期图法但它的问题是方差大谱线毛刺多。工程上更推荐Welch法也就是pwelch函数。它把信号分段、加窗、做FFT、再平均以方差为代价换取平滑的谱估计。资源包里psd_feature_extract.m用的是pwelch并在此基础上进一步提取谱特征。这部分的思路是先做PSD估计再从PSD曲线上找峰值频率、峰值的幅值以及在特征频率附近的边带能量。% psd_feature_extract.m % 使用Welch法估计功率谱密度提取谱峰特征 [x, fs] audioread(sample_vibration.wav); % 读音频文件实际项目中可能是其他格式 x x(:, 1); % 取单通道 [pxx, f] pwelch(x, hann(1024), 512, 4096, fs); % 窗长1024重叠512NFFT 4096 % 查找主峰频率 [peak_val, peak_idx] max(pxx); peak_freq f(peak_idx); % 在0.5~2倍主峰频率范围内查找次峰 band_mask f 0.5*peak_freq f 2*peak_freq; band_pxx pxx(band_mask); band_f f(band_mask); [second_val, second_idx] max(band_pxx); second_freq band_f(second_idx); % 输出特征向量 psd_features [peak_freq, peak_val, second_freq, second_val, ... sum(pxx), sum(pxx .* f) / sum(pxx)]; % 最后一个是重心频率(基于PSD)这段代码里我故意用一个.wav文件做输入示例因为audioread可以直接返回采样率省掉手动指定的麻烦。pwelch的参数设置是个门道窗长决定频率分辨率窗越长频率分辨率越高但时间分辨率越低重叠量一般取窗长的50%~75%NFFT是FFT点数建议取大于窗长的2的幂次这样可以获得更细的插值频率网格但注意它并不提高真实分辨率。我一般习惯hann(1024)窗长配512重叠这种组合在大多数振动信号上表现比较稳。如果你处理的是低频信号比如转速只有几百转的设备需要把窗长加大到4096甚至更长否则低频段的频率分辨率不够特征频率附近的边带根本看不出来。3.3 频段能量的工程解释边带分析频段能量这个特征在轴承故障诊断里特别有价值。滚动轴承出现局部损伤时振动信号会以某个特征频率为中心、两侧出现间隔为转频的边带。这些边带能量占总能量的比例和故障严重程度高度相关。所以我在实际处理中会把频段划分得不是等宽的而是围绕轴承故障特征频率的窄带。比如转频是30Hz外圈故障特征频率是210Hz那我就取200Hz~220Hz这个窄带的能量作为特征。这需要在知道轴承参数和转速的前提下才能做但如果数据里带转速信息这个特征往往比盲目的等宽频带能量有效得多。资源包里提供的band_energy是等宽频段划分这是一个通用方案适合不知道具体轴承包参数的情况。如果你的应用场景知道特征频率完全可以修改程序里的band_width和num_bands参数改成窄带能量提取。4. 把时域和频域特征合成一个向量联合提取脚本与归一化的时机4.1 联合提取怎么组织代码在实际项目里我从来不会只单独用时域特征或频域特征而是把两类特征拼成一个向量后续既可以喂给SVM做分类也可以作为统计模型的特征输入。资源包的extract_all_features.m就是干这个事的在一个脚本里依次调用时域特征和频域特征函数把结果拼接成行向量。% extract_all_features.m function feature_vector extract_all_features(signal, fs) % signal: 原始振动信号, fs: 采样率 % feature_vector: 时域13维 频域3维 频带能量10维 26维特征向量 td time_domain_features(signal); % 时域 1x13 [fc, rmsf, stdf, band_energy] spectral_features(signal, fs); % 频域特征 % 逐个拼起来 feature_vector [td, fc, rmsf, stdf, band_energy]; % 1x26 end这种写法最简单直接但有一个隐患如果signal是一个矩阵而不是向量——某些采集卡导出数据时每列是一个通道——time_domain_features会返回1xN*13维的向量和后面的频域特征拼接时会报维度不匹配错误。我在调用前都会加一句signal signal(:, 1);来强制取单通道。另外要注意的是这段代码用[td, fc, rmsf, stdf, band_energy]拼接时td是1×13band_energy是1×10其他三个是标量拼接结果维度是1×26。如果你后续想把这个特征向量存成表格方便查看或导入其他工具建议用writematrix直接写成CSV或者用table带上列名不然输出就是一排数字根本分不清哪个是哪个。% 导出特征表格 feature_vector extract_all_features(signal, fs); T array2table(feature_vector, VariableNames, ... {mean, peak, abs_mean, root_amp, std, rms, kurtosis, ... waveform_factor, peak_factor, pulse_factor, margin_factor, ... fc, rmsf, stdf, band1, band2, band3, band4, band5, ... band6, band7, band8, band9, band10}); writetable(T, features_output.csv);给特征命名这件事看起来是个小动作但后面排查特征重要性时帮大忙。我踩过的坑是一开始输出特征不带名称用机器学习模型算特征重要性只能靠索引找某个特征对应哪个公式全靠回忆痛苦到怀疑人生。4.2 归一化应该放在特征提取后还是训练前这个问题几乎每个用特征向量做后续分析的人都会问。答案是归一化永远在特征拼接完成之后、交给模型训练之前做不要在特征提取阶段做。原因很直接归一化需要用到数据集的全局均值和标准差你把所有样本都拿到手之前这个均值是不存在的。如果每个样本单独归一化等于把每个样本的幅值信息彻底抹掉了——而幅值本身就是很重要的故障信息。具体做法是等特征矩阵X组装完成后用zscore做标准化% 在所有样本的特征都提取完后执行 X feature_matrix; % Nx26N是样本数 [X_norm, mu, sigma] zscore(X); % 每列减去均值除以标准差zscore返回三个值后两个是每列的均值和标准差需要在测试阶段使用同样的数值做归一化。很多人在实际部署时忘记保存mu和sigma等到新数据来了无法归一化这就是经典的“模型上线翻车”现场。4.3 特征冗余的坑峭度和峰值因子高度相关你可能已经注意到时域特征里的峰值因子、脉冲因子、裕度因子本质上都是“峰值除以某种幅度统计量”。这导致它们在物理含义上高度冗余。实际测试中峰值因子和脉冲因子的皮尔逊相关系数常常超过0.95两个特征等于一个特征的信息。在这种情况下我不建议把14个特征一股脑全输出。更实际的做法是先用corrcoef计算相关性矩阵找到相关系数超过0.9的特征对只保留其中一个。另一个更自动化的方式是PCA但PCA的结果是原有特征的线性组合可解释性变差不利于你向非技术背景的同事解释“为什么这个特征能反映故障”。所以我在生产项目里基本是手动做特征筛选把相关性高的特征删掉后再做分类。5. 特征提取避坑指南采样率陷阱、短样本FFT分辨率、NaN数据与参考标准5.1 采样率标错导致频率轴整体偏移现象提取出来的重心频率、主峰频率和理论计算的特征频率对不上比如轴承外圈故障频率理论值是210Hz程序算出来是205Hz而且转速越低偏差越大。原因MATLAB里读数据时采样率参数传错了。常见来源是数据采集软件导出的文件头部信息标注的采样率与实际采集卡设置不一致或者读CSV时只读了数据列采样率标称值被当成默认值。解决先用一段已知频率的信号验证采样率。比如拿一个50Hz的工频信号跑fft看峰值是否准确落在50Hz上不对就修改采样率参数。还有一种隐蔽情况采集卡实际采样率是51.2kHz也就是所谓的“整数倍采样”但导出软件标注的是50kHz这会导致约2.4%的频率误差。我的习惯是启动程序后第一件事打印fs做一次校验再继续往下算。5.2 短样本FFT的频率分辨率不够特征频率算不准现象对一段0.1秒的信号做fft频率轴的间隔是10Hz轴承故障特征频率是105Hz结果频谱峰值落在100Hz或110Hz两个都差5Hz算出的边带能量偏差更大。原因FFT频率分辨率等于fs/NN是样本点数。样本时长越短N越小频率网格越粗。解决在pwelch中通过加大NFFT参数来做插值或者回到数据采集阶段确保每段样本至少采集1秒以上的信号。对1024Hz采样率、1秒时长频率分辨率是1Hz足够分辨大多数机械故障特征频率。如果信号本身就只有0.1秒那就放弃FFT路线只用时域特征强行算频谱只会得到一团无法解释的噪声。5.3 数据里有NaN或Inf特征计算结果全是NaN现象time_domain_features返回一堆NaN检查原始波形却发现只有一个采样点异常。原因传感器断线或采集卡饱和时某个采样点的数值可能是NaN、Inf或远大于正常量级的数字。这些异常值会污染整个统计计算尤其是均值、RMS对极端值非常敏感。解决在特征提取入口处做一个数据清洗先把Inf转为NaN再用fillmissing做线性插值或者直接丢弃该段数据。我一般用以下写法x(isinf(x)) NaN; if any(isnan(x)) x fillmissing(x, linear); % 线性插值填坑 warning(检测到NaN已用线性插值处理); end注意这个处理必须在time_domain_features之前执行并且处理后要重新观察一下波形确认插值没有引入不合理的假脉冲。特别是原始信号里如果存在很长的NaN连续段线性插值会在中间制造一条直线这会直接削弱故障脉冲特征。5.4 峭度/均值等特征没有公认阈值自设报警限需要小心现象把峭度算出来后想设报警阈值发现不同标准给出的参考值不一样有的说峭度大于3算故障有的说大于5到底听谁的原因峭度值受信号预处理方式影响很大。原始信号经过高通滤波后低频大能量成分被削弱脉冲特征相对突出峭度值会系统性偏大。资源包里的代码直接用原始信号算峭度如果你的后续处理加了滤波峭度参考值当然要重新标定。解决先明确自己的峭度是“原始信号峭度”还是“滤波后峭度”然后采集正常工况下的足够多样本做基线分布用均值加3倍标准差作为报警线比直接用文献里的固定阈值可靠得多。我一般要求不低于50组正常样本做基线统计。5.5 不同文件的数据长度不一致特征矩阵拼不起来现象批量提取特征时feature_matrix行数总是和file_names对不上或者某个文件的特征向量维度比其他文件多一截。原因部分mat文件里的信号长度明显不同而time_domain_features对长度不敏感但spectral_features的分段逻辑可能导致输出的band_energy维度不同——比如某段信号长度不足有效频段数不足10个。解决在批处理脚本开头加一个统一长度的重采样逻辑把所有信号截断或重采样到固定长度。最简单的做法是统一取前N个点N是某个固定值比如fs * 10即10秒长度。如果信号长度不足N直接丢弃该文件并打印警示。这样能保证特征矩阵严格对齐不会出现后续分类器报维度错误。6. 最后一招把特征提取结果和人工频谱计算做交叉验证特征提取代码跑通了不代表结果的物理意义是对的。我的习惯是每换一批新数据强制做一次交叉验证——拿其中一段信号用MATLAB的交互式工具spectrogram或直接画出幅值谱肉眼数出主峰的频率值再和程序输出的主峰频率对比。如果差异超过一个频率分辨率单位就说明程序里有参数错误不查清楚绝不继续往下做。具体的验证方法是这样取一段信号先用findpeaks找出幅值谱前5个峰值的频率和幅值打印出来和spectral_features的输出对比。如果主峰频率对上了但band_energy有异常可以放大某个频段单独画图检查。这个步骤看起来是多余的但在新数据集上做过一次之后你对这套特征提取代码的信任度会完全不一样。从那以后我每次换数据集都强制走一遍这段验证流程宁可多花半小时也不带着一个不确定的参数跑完整条流水线。希望帮到你。本文还有配套的精品资源点击获取