企业官网建设流程全解析

简介：面向信号处理学习者和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{end+1} = 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，重叠512，NFFT 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×13，band_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; % Nx26，N是样本数 [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/N，N是样本点数。样本时长越短，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有异常，可以放大某个频段单独画图检查。这个步骤看起来是多余的，但在新数据集上做过一次之后，你对这套特征提取代码的信任度会完全不一样。从那以后我每次换数据集，都强制走一遍这段验证流程，宁可多花半小时，也不带着一个不确定的参数跑完整条流水线。希望帮到你。

本文还有配套的精品资源，点击获取