Matlab时域频域特征提取实战：从峭度到频谱质心的核心方法

简介：一套面向MATLAB用户的时域与频域特征提取实战资源，适合信号处理初学者以及需要做数据分析的学生、工程师。它从原始信号出发，系统讲解并实现了时域统计特征（均值、峰值、均方根、峭度、尖峰因子）与频域特征（频谱中心频率、带宽、功率谱密度、谱熵等）的计算，并提供了完整的特征组合思路，可用于设备故障诊断、语音识别、医学信号分析等应用场景。压缩包共10个文件，包含6个mat数据文件和4个m脚本文件，整体约20.68MB；其中数据文件提供可直接运行的样本，脚本文件覆盖时域统计、频域统计与整体特征提取三个环节，结构清晰、便于二次复用。目前已有2024人学习，代码注释与模块划分对初学者较为友好。通过学习本套资源，读者能掌握从信号读入、预处理、特征计算到特征向量输出的完整流程，并可将提取的多维特征直接用于后续分类或回归任务，是一份理论结合实践的实用资料。

1. 认识“matlab时域频域特征提取.zip”：为什么值得拆开而不是直接跑

很多人拿到类似“matlab时域频域特征提取.zip”这样的代码包，第一反应就是解压、run、出图，然后完事。但如果你真的要做故障诊断、设备状态监测或者生物信号分析，就会发现事情远没有这么简单。这个包的核心是把一段原始信号浓缩成一组数字：时域均值、方差、峭度，频域幅值谱峰值、频带能量、频谱质心……这些数字看似平淡，却是后续机器学习模型或阈值判断的真正输入。决定模型好坏的不是分类算法，而是你从信号里提取了什么特征。这篇笔记里，我用自己常用的一套Matlab逻辑，把时域和频域特征提取的选型理由、核心代码、参数边界和踩坑点一次讲透。新手可以照着我给的函数跑起来，熟手也能在参数和边界上少走点弯路。

2. 时域频域特征到底在提取什么：定义与选型依据

2.1 时域特征：从统计量到无量纲指标

时域特征直接从原始信号的时间序列上计算，不经过傅里叶变换，所以计算成本低、物理意义直观。最基本的几个是均值、峰峰值、均方根和标准差。均值反映信号的直流分量，标准差反映波动幅度，均方根则同时包含幅值和能量信息，在振动和电流信号里几乎是必算项。峰峰值对冲击类故障很敏感，但也容易被离群值带偏。

再往上走就是峭度、波形因子、峰值因子、脉冲因子和裕度因子。峭度是四阶中心矩与标准差的四次方之比，信号接近正态分布时峭度接近3；滚动轴承早期故障往往产生稀疏冲击，峭度会显著升高。波形因子是均方根除以整流平均值，峰值因子是峰值除以均方根，这两个对磨损类退化相对稳定。脉冲因子和裕度因子对单次冲击更敏感，但也更容易受噪声影响。

特征不是越多越好。举个例子，标准差和均方根在信号零均值时高度相关，波形因子和峰值因子又可以用均方根和峰值互相推导。我在模拟项目X里试过同时放入12个时域特征，结果相关矩阵里好几对系数超过0.95，模型学到的冗余远比信息多。所以选时域特征前，先确认信号是否做过零均值化，零均值后标准差和均方根可以只留一个。

计算时域特征时，数据的长度对数值稳定性影响很大。峭度需要四阶矩，样本量太少时方差本身就不稳，峭度更容易爆炸。实践里我一般要求每条样本至少包含20个以上的信号周期。比如转频25Hz的信号，采样2秒就是50个周期，足够；如果只有0.2秒，5个周期，峭度基本没法用。另外，去直流后再算均方根才是真正的交流有效值，否则算出来的RMS里混着直流偏置，故障分量会被掩盖。

2.2 频域特征：从幅值谱到频带能量

频域特征通过对信号做FFT得到。对实数信号来说，N点FFT产生N/2+1条有效谱线，频率分辨率是fs/N。幅值谱的纵轴是各频率分量的幅值，功率谱密度则把幅值平方归一化后按频率分布。实际使用中，我不太关心每条谱线的精确值，更常用几个聚合指标：主频幅值、频谱质心、谱峰位置，以及特定频带内的能量占比。

频谱质心的数学含义是频率的加权平均，频率权重用幅值或功率谱值。破损齿轮的频谱质心往往会从正常时的低频区往高频区偏移，这个指标比单看峰值更稳健。频带能量则是把FFT结果按你关心的区间切段求和，常用于电机电流信号：比如把0-50Hz、50-100Hz、100-200Hz分别累加，得到3个特征。注意频带划分必须结合设备转频和故障特征频率，不能随便切。

频域特征提取前，必须知道采样率fs和样本长度N。fs决定了能看到的最高频率（奈奎斯特频率fs/2），N决定了频率分辨率df=fs/N。如果你想分辨2Hz以内的边频带，就需要保证N在时域上覆盖至少1/2秒（对应df=2Hz）。这个在后续参数调整里会反复遇到。

2.3 特征选型的三步决策法

第一步，明确你的信号类型和目标。振动信号优先考虑峭度、RMS、频谱质心和特征频带能量；电流信号更关注基波附近边频带能量；声学信号可以先算频谱质心和频带能量。信号类型决定了哪些特征有物理意义，那些纯粹数学上成立的指标如果和故障机理无关，最后只会变成噪声。

第二步，先算一个宽特征集，比如时域10个加频域5个，然后用相关性矩阵剔除相关系数大于0.95的重复特征。这一步不需要复杂算法，Matlab的corrcoef就够了。排序后把冗余特征删掉，保留与故障机理直接相关的那些。这一步做完，特征维度通常能从15个降到8个左右。

第三步，用箱线图或散点图看正常样本和故障样本在每个特征上的分布重叠度。如果两类样本在某个特征上完全重叠，直接删掉；如果明显分开，保留。我一般会写一个自动筛选：对每个特征计算两类均值差除以标准差，得到分离度指标，大于1的才进入模型。这三步做下来，特征维度通常从15个降到6到8个，模型准确率反而上升。别迷信“特征越多越丰富”，落地时训练数据有限，冗余特征只会加剧过拟合。

3. 用Matlab实现时域频域特征提取：核心函数与最小可运行代码

3.1 搭一个可复用的特征提取框架

我习惯封装成两个函数：一个算时域特征，一个算频域特征，外层脚本负责读数据和组装表格。这样换数据集时不用改特征代码，换特征时不用动数据读取逻辑。

下面这个函数是主入口，输入一段信号x和采样率fs，输出一个1行多列的特征向量。

function feat = extract_features(x, fs) % 输入: % x : 列向量信号，单位不限，建议先校准到物理单位 % fs : 采样率，单位 Hz % 输出: % feat : 结构体，包含时域和频域特征 x = x(:); % 强制转成列向量 x = x - mean(x); % 去直流，避免直流分量淹没交流特征 feat.t_rms = rms(x); feat.t_std = std(x); feat.t_peak = max(abs(x)); feat.t_kurt = kurtosis(x); feat.f_main = calc_main_freq(x, fs); feat.f_centroid = calc_freq_centroid(x, fs); end

这段代码里，去直流是第一个关键动作。如果不做这一步，直流分量会出现在FFT的0Hz处，而且时域rms值包含了直流偏置，后续所有无量纲指标都会失真。kurtosis函数用的是Matlab统计工具箱自带实现，如果你用的是基础版，可以改用mean((x-mean(x)).^4)/(std(x)^4)。

参数说明：feat是结构体，比数组直观；后期如果要拼成训练矩阵，用struct2array转出值即可。calc_main_freq 和 calc_freq_centroid 需要我们自己定义，下面两节会给出。

3.2 时域特征提取函数：代码逐段拆解

这是时域特征的核心函数，我一次算全，方便你看清楚每个特征的边界。

function f = calc_time_features(x) % 时域特征提取 % 输入 x 为去直流后的信号列向量 n = length(x); f.mean = mean(x); % 均值（去直流后应接近0） f.rms = sqrt(sum(x.^2)/n); % 均方根，手工计算避免依赖工具箱 f.std = std(x, 0); % 标准差，0表示除以 n-1 f.peak = max(abs(x)); % 峰值 f.peak_to_peak = max(x) - min(x); % 峰峰值 f.crest = f.peak / f.rms; % 峰值因子 f.kurt = mean((x - f.mean).^4) / (f.std^4); % 峭度 f.form = f.rms / (mean(abs(x)) + eps); % 波形因子，eps防除零 end

第一行注释就把“去直流后的输入”作为前提，这个约定要在整个工程里统一。rms我故意不用Matlab自带的rms函数，而是写sqrt(sum(x.^2)/n)，因为在批处理循环里调用自带函数有额外开销，且手工写法能精确控制分母是n还是n-1。std(x,0)的第二个参数0代表除以n-1，如果你希望和RMS保持同一口径，可以改成std(x,1)，这一点很多人会踩坑。

eps的出现是因为整流平均值在某些极限样本里可能接近0，加了它之后波形因子不会变成Inf。实际中如果一段信号几乎为0，波形因子本身没有意义，这时候应该更早地在数据清洗阶段把它过滤掉，而不是靠eps兜底。

补充一下峭度的边界：峭度对异常尖峰非常敏感，一个毛刺能让峭度从3跳到几十。所以在算峭度之前，我建议先对信号做中值滤波或剔除明显的离群点。但注意，轴承故障的冲击本来就是你要提取的信息，不能无脑滤波。工程上通常只剔除超过6倍标准差的样本点，保留4到6倍之间的冲击。

3.3 频域特征提取函数：FFT与频带能量

频域部分的核心是正确使用FFT，并处理直流和窗函数。

function f = calc_freq_features(x, fs) % 频域特征提取 % 输入: x 为去直流后的信号, fs 为采样率 n = length(x); X = fft(x) / n; % 归一化幅值谱，除以n后幅值正确 mag = abs(X(1:floor(n/2)+1)); % 取单边谱 mag(2:end-1) = 2 * mag(2:end-1); % 单边谱幅值修正 freq = (0:floor(n/2)) * fs / n; % 频率轴 [fval, idx] = max(mag(2:end)); % 找主峰，跳过0Hz f.f_main = freq(idx+1); f.f_centroid = sum(freq .* mag) / sum(mag); % 频谱质心 f.band1 = sum(mag(freq >= 0 & freq < 50)); % 低频段能量 f.band2 = sum(mag(freq >= 50 & freq < 100)); f.band3 = sum(mag(freq >= 100 & freq < 200)); f.energy = sum(mag.^2); % 总能量 end

注意这里有几个关键参数。FFT归一化：fft(x)/n是对长度归一化，使得正弦幅值在频谱上能恢复真实幅值。单边谱修正：对实数信号，除了0Hz和奈奎斯特频率外，每个幅值要乘2，因为负频率的能量折叠到了正频率。freq = (0:floor(n/2)) * fs / n 给出了每条谱线对应的物理频率，单位Hz。

频带划分里的50、100、200是我随便写的示例，实际必须根据你的设备特征频率来改。比如电机电流信号，工频50Hz，你关心的是50Hz附近的边频，那频带应该设成45~55Hz、55~65Hz这样，而不是统统从0到200。在代码里 band1 = sum(mag(freq >= 0 & freq < 50)); 这种写法会包含0Hz，所以尽管前面去了直流，但频谱泄漏仍可能让0Hz能量漏进来。更稳妥的是用 freq > 0.5 之类的小阈值跳过0Hz。

频谱质心 sum(freq .* mag)/sum(mag) 用幅值加权，如果你更关心能量分布，可以用 sum(freq .* mag.^2)/sum(mag.^2)，两种各有侧重。幅值加权对主频附近的变化更敏感，能量加权对宽带噪声更敏感。

注意：频带划分必须结合设备特征频率，不能照搬示例的50、100、200。先把设备转频、故障特征频率算清楚，再写进band参数。

3.4 批量处理多个样本：循环与表格输出

在实际项目里，你不会只处理一段信号，而是面对一个文件夹里几百个csv。下面这段脚本循环读取样本、提取特征、拼成训练表格。

files = dir('data/*.csv'); n = length(files); feature_names = {'mean', 'rms', 'std', 'peak', 'peak_to_peak', ... 'crest', 'kurt', 'form', 'f_main', 'f_centroid', ... 'band1', 'band2', 'band3', 'energy'}; features = zeros(n, numel(feature_names)); for i = 1:n data = readmatrix(fullfile('data', files(i).name)); x = data(:, 2); % 假设第二列是信号 fs = 2560; % 采样率，根据实际改 x = x - mean(x); t = calc_time_features(x); f = calc_freq_features(x, fs); features(i, :) = [t.mean, t.rms, t.std, t.peak, t.peak_to_peak, ... t.crest, t.kurt, t.form, f.f_main, ... f.f_centroid, f.band1, f.band2, f.band3, f.energy]; end T = array2table(features, 'VariableNames', feature_names); writetable(T, 'features.csv');

这个脚本里，dir('data/*.csv')依赖当前工作目录，最好用fullfile('data', '*.csv')拼接，避免路径分隔符在不同系统上出错。data(:,2) 是读取信号列，如果你的csv有两列时间戳和信号，就用第2列；只有一列则改成 data(:,1)。fs这里写死2560，只是示例，实际要从数据采集参数里读，不能凭感觉。

批量处理最怕的是某个样本长度不一样。下一章会专门讲这个坑。如果你提前知道了长度不一致，可以在循环里加一句x = x(1:min(len, target_len));做截断或补零，但补零会带来频谱泄漏，截断则可能丢失信息。我一般优先截断，因为样本边界处的信息往往不是关键。

4. 避坑指南：时域频域特征提取的5个典型问题

4.1 采样率不归一：特征值忽大忽小

现象：同一设备同一工况，今天提取的RMS是0.5，明天变成0.8，频谱主频也从100Hz变成101Hz，但设备本身没坏。

原因：数据采集系统在重采样后实际采样率发生了微小偏移，或者你的代码里fs写错了。特征提取函数里所有频域坐标都依赖fs，RMS这类时域特征虽然不直接依赖fs，但当样本长度固定时，fs变了等于观测时长变了，统计量自然变。频域的f_main和f_centroid对fs的错误尤其敏感。

解决：把fs作为参数从采集配置里读，不要写死。如果信号是从不同设备来的，统一重采样到同一个fs。Matlab可以用resample(x, p, q)，按p/q比例重采样。注意重采样后要重新检查零均值，因为重采样滤波器会引入瞬态。另外，处理前先打印一行fprintf('fs = %d Hz, n = %d\n', fs, length(x));，确认输入确实是你期望的。这个习惯能帮你抓到一大半“特征值漂移”的问题。

4.2 窗函数漏选：频谱泄漏把峰值抹平

现象：正弦波的理论FFT应该是单一尖峰，实际得到的谱峰扩展到附近几个频点，主频幅值也偏低。

原因：对有限长度信号直接FFT，等价于把无限长信号乘以矩形窗，矩形窗的旁瓣很高，导致能量泄漏。尤其当信号频率不是频率分辨率的整数倍时，泄漏更严重。

解决：先用汉宁窗或汉明窗对信号加窗，再去FFT。加窗后主瓣变宽，但旁瓣大幅降低。代价是频率分辨率下降，所以加窗后做幅值归一化时要除以窗函数的平均值，否则幅值会偏低。我一般会这么写：

w = hanning(n); xw = x .* w; xw = xw / mean(w); % 保持能量一致

注意，加窗后时域特征不能再从xw上算，时域特征必须用原始信号x。所以流程是先算时域特征，再去趋势、扣窗、做FFT。如果你发现加窗后主频幅值还是对不上，检查一下是否忘了除以窗均值。

提示：加窗后时域特征必须用原始信号计算，不能复用加窗后的序列，否则RMS和峭度会整体偏小。

4.3 没去趋势项：低频分量把频域特征带偏

现象：频段能量里band1（低频段）异常大，频谱质心被拽向0Hz，即使你心里清楚信号不该有这么大的直流或超低频。

原因：传感器漂移、温漂或者信号里有缓慢变化的趋势项。FFT会把趋势项视为极低频分量，能量全部堆在头几根谱线上，还会通过泄漏污染周围频段。

解决：在频域提取前，先对信号做去趋势。简单做法是x = detrend(x, 'linear')，它对一条直线趋势有效；如果趋势是二次曲线，用polyfit拟合后减去。去趋势后再做零均值，两步不能同时省。比如传感器信号有线性漂移，直接减mean只去掉一个常数，剩余斜坡在FFT看起来还是一个大斜坡成分。我通常的流程是：先detrend再去直流，最后才送进特征函数。

4.4 峭度遭遇离群值：一个毛刺毁掉整条特征

现象：某条样本的峭度从3跳到15，所有后续模型都把这个样本当成异常，但看时域波形只是多了一个采集时的电压毛刺。

原因：峭度对四阶矩敏感，一个离群点产生的权重远大于普通样本点。

解决：在特征提取前做离群点处理。我常用的策略是先计算信号中位数和绝对偏差，把超过中位数 + 6*MAD的样本点视为毛刺，将其替换为前后两个点的均值。对轴承故障这种本身带冲击的信号，阈值不能太低，否则会把真实故障冲击滤掉。如果不确定，可以先画波形把毛刺标记出来，再决定阈值。离群点剔除只作用于时域特征计算之前，频域特征建议用剔除后但未加窗的信号来计算，或者用鲁棒性更好的频带能量而不是单峰幅值。

4.5 样本长度不一致：批量特征矩阵维度错乱

现象：批量循环跑到第37个文件时报错，或者features矩阵出现NaN，而单独跑每个文件都是好的。

原因：不同文件长度不同，fft点数不同，频域特征的频率轴长度自然不同。某些文件中某个频带可能没有足够多的谱线，导致sum空集产生NaN。更隐蔽的是，时域特征的计算结果本身不受长度影响，但后续机器学习要求所有样本特征向量长度一致，长度不一致连接起来就会有问题。

解决：统一样本长度。常见做法是截断到固定长度N，比如所有样本取前1024点。如果样本太短，则扩展是x = [x; zeros(N-n,1)]补零，但补零后要标注该样本是截断还是补零的，后续训练时可作为权重项。我一般更倾向截断，因为补零会引入额外的低频能量，导致频域特征前后不一致。截断时注意保留信号的完整周期，最好让N等于采样率除以转频的整数倍，或者直接用整周期截断。如果实在没法保证，就在特征表里加一列“样本长度”，方便后期排查。

5. 快速自检：用解析信号验证你的特征提取结果

5.1 用正弦波做基准验证

构造一个已知幅值A、频率f0的正弦波，加一点噪声。它的理论RMS是A/sqrt(2)，FFT主频就是f0。拿这个信号过一遍你的特征提取代码，对比理论值和输出值。如果RMS偏差超过0.1%，说明你的归一化或去直流有问题；如果f_main不是f0，说明fs或频率轴写错了。

我常用这个自检脚本：

fs = 1000; t = (0:999)' / fs; x = 2 * sin(2*pi*50*t); % 频率50Hz，幅值2 feat = extract_features(x, fs); assert(abs(feat.t_rms - 2/sqrt(2)) < 0.01); assert(abs(feat.f_main - 50) < 0.1); disp('自检通过');

这个自检脚本我会放在所有特征提取函数修改之后。以前我拿到一个现成的时域频域特征提取包，习惯直接信它的计算，后来发现频域主频差了1.5Hz，排查半天才发现是采样率写反。从那以后，我要求所有特征提取代码必须附带一个这样的解析信号自检。这个习惯帮我省下的调试时间，远超写脚本本身花掉的那十分钟。

5.2 用散点图快速看特征分离度

特征算出来不等于可用，先画一张正常样本和故障样本的二维散点图，比如横轴RMS，纵轴峭度。如果两类点明显分区，这个特征组合值得进模型；如果完全混在一起，趁早换特征。这一步不花多少时间，却能在模型训练前就把不可用的特征淘汰掉。

gscatter(features(:,2), features(:,7), labels); xlabel('RMS'); ylabel('Kurtosis');

labels 是样本类别标签，0代表正常，1代表故障。散点图能直观暴露特征是否有区分度，也能发现离群样本。我通常会在建模前对每一对特征都画一遍这种图，花十五分钟省下后面调参几小时。记住，特征提取不是跑完代码就结束，验证和可视化才是它真正生效的环节。希望帮到你。

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