加速度信号PSD计算与GRMS校核:FFT到Welch法实战解析
2026/9/15 22:31:58 网站建设 项目流程

简介:面向需要使用MATLAB开展振动信号分析与功率谱密度计算的工程技术人员与信号处理学习者,压缩包提供了一个可直接运行的算法实例,重点解决加速度信号的sin波模拟与PSD求解问题。资源覆盖从信号生成、窗函数处理、FFT变换到功率谱归一化、对数展示的完整流程,适用于机械振动监测、地震信号分析等场景。包内共2个m文件,整体约1KB,代码简洁,便于快速阅读和二次修改。已有284人学习下载,适合希望用短小实例快速上手MATLAB加速度PSD分析的用户。通过运行这两个脚本,能够直观掌握基于fft的加速度信号功率谱估计方法,并将频率轴构建、功率谱归一化等关键步骤迁移到自身数据中。资源虽小,但完整呈现了从时域建模仿真到频域特征提取的链路,对提升信号处理理论与工程实践结合能力具有不错的参考价值。

1. 为什么加速度信号的 PSD 要这样算:从 GRMS 指标反推需求

做振动试验的人经常会在试验大纲里看到 GRMS 这个指标,它是随机振动中总均方根加速度,单位是 g 或 m/s²。拿到一套 MATLAB 的加速度功率谱密度算例,里面只有GRMS1_GRMS2.mBuffet_sin.m两个脚本,乍看名字像是在算时域 RMS 和造正弦波,实际上这两件事正是 PSD 分析的两端:信号怎么模拟、能量怎么校核。故障诊断和结构疲劳评估里经常遇到这样的场景,时域波形看起来差不多的两段加速度数据,转到频域后一个能量集中在 20 Hz 窄带,一个铺散在 200 Hz 宽带,判断依据只能是 PSD 曲线和 GRMS 数值。下面就把从 sin 信号构造、FFT 估计、窗函数处理到 PSD 积分校核的完整链路走一遍,适合正在处理振动台数据、想把加速度功率谱密度彻底算明白的工程师。

2. PSD 估计的数学底子:FFT、窗函数与单边谱归一化

2.1 周期图法在算什么:从时域总能量到频域谱密度

PSD 的单位是加速度单位的平方除以频率,常见的是 g²/Hz 或 (m/s²)²/Hz。把 0 到 Fs/2 频率范围内的 PSD 曲线积分再开方,就得到总均方根加速度 GRMS,这一性质让 PSD 成为随机振动试验中最常用的验收曲线。周期图法是最直接的一种估计思路:对 N 点加速度序列加窗,做 FFT,取幅值平方后除以 Fs·N,得到双边功率谱密度。很多人在这一步习惯直接写abs(fft(x)).^2/N,再乘频率分辨率 df=Fs/N,两个写法数值完全等价,区别只是归一化放在前面还是后面。容易翻车的是单边谱合并:FFT 输出以 Fs/2 对称,绘图时通常只保留 0 到 Fs/2 的 N/2+1 条谱线,除 DC 和 Nyquist 频率外每条谱线能量要乘 2,否则 GRMS 校核结果会偏小大约根号 2 倍。

% 周期图法估计加速度信号 PSD,输入 x 为加速度时域序列 x = x(:); % 统一为列向量 N = length(x); % 采样点数 Fs = 1000; % 采样率,单位 Hz w = hann(N, 'periodic'); % 周期汉宁窗,降低频谱泄漏 xw = x .* w; % 加窗 X = fft(xw); % FFT 得到复数频谱 P2 = abs(X).^2 / (Fs * N); % 双边 PSD,单位 (m/s2)^2/Hz P1 = P2(1:N/2+1); % 取单边谱 P1(2:end-1) = 2 * P1(2:end-1); % 合并负频率能量 f = (0:N/2).' * Fs / N; % 单边频率轴

代码里hann(N,'periodic')是周期窗形式,适合谱估计;普通hann(N)是对称窗,边界不完全归零,周期延拓时会产生细微跳变。Fs*N这个分母决定了 PSD 的数值尺度,换采样率或换点数后 GRMS 积分结果应当不变,这是一个很好的自检手段。把 Fs 和 N 成倍改变后重算,如果 GRMS 与原来一致,说明归一化写对了。

2.2 加窗为什么是必须的:频谱泄漏与旁瓣衰减

对有限长信号直接做 FFT,相当于用矩形窗截断原始序列。矩形窗的频率响应旁瓣只衰减约 13 dB,当信号能量较大又不在整数频率点上时,这些旁瓣会污染邻近频段,看起来像多出一堆小峰值,这就是频谱泄漏。汉宁窗是第一选择,旁瓣衰减约 31 dB,主瓣宽度适中;汉明窗旁瓣衰减更陡,但第一旁瓣特性与 hann 不同,常用于语音这类时序分析。平坦顶窗峰值测量最准,代价是主瓣很宽,频率分辨能力差,适合标定正弦幅值而不是观察宽带随机激励。

窗函数主瓣宽度(×Fs/N)第一旁瓣衰减典型用途
矩形窗213 dB瞬态冲击或整周期截断
hann431 dB随机振动 PSD 默认选择
hamming443 dB窄带分析
flattop895 dB正弦峰值标定

信号是正弦时,如果采样时长恰好是周期的整数倍,加不加窗对幅值影响不大;随机振动没有周期性可言,任何截断都会引入泄漏,所以 Welch 分段平均里每个分段都要先加窗。我一般先用不重叠的周期图扫一遍,确认没有异常尖峰,再用加窗结果做正式报告曲线,这样能直观看出窗函数对曲线的平滑作用。

2.3 频率分辨率 df 与有效采样时长怎么匹配

频谱里相邻两条谱线的间距 df = Fs / N,想区分两个相隔 Δf 的正弦分量,条件大致是 df 小于 Δf。比如要分辨 49.5 Hz 和 50.5 Hz,采样率 1000 Hz 时 N 至少要大于 1000 点,对应 1 秒以上的数据。df 越小谱线越密,但单根谱线的统计方差会变大;随机振动 PSD 估计不追求单根谱线的绝对精度,更看重整段曲线的统计稳定性。处理实测数据时,先根据目标频率间隔确定最小数据长度,再考虑重叠率与平滑度,不要一上来就用最大点数做单次 FFT,那样曲线毛刺会非常严重。

3. 复现 apsd 实例:Buffet_sin.m 信号生成与 GRMS1_GRMS2.m 积分校核

3.1 用正弦叠加模拟抖振加速度信号

Buffet_sin.m里的 Buffet 指抖振,常见于飞机尾翼、汽车外后视镜这类结构,流体分离激励下的加速度响应往往集中在几个窄带频率附近。用 sin 函数把若干频率分量叠加起来是模拟这类信号最直观的做法:主频给一个较大的幅值和初始相位,次级频率给较小幅值和不同相位,还可以再加一点白噪声模拟传感器噪声与随机气流激励。

% 构造 5 秒抖振模拟加速度信号 Fs = 1000; % 采样率 1000 Hz t = (0:Fs*5-1).' / Fs; % 0 到 4.999 秒,共 5000 点 x = 1.2 * sin(2*pi*20*t + 0.3) + ... % 20 Hz 主分量 0.6 * sin(2*pi*25*t) + ... % 25 Hz 次级分量 0.3 * sin(2*pi*80*t); % 高频分量

参数设置上,20 Hz 主分量代表结构一阶模态附近的抖振,25 Hz 代表稍高的激励带,80 Hz 用来观察高频段的谱线形态。幅值如果按 g 读取,GRMS 算出来也直接是 g,与振动台试验大纲对齐很方便。注意时间轴用(0:Fs*5-1).'生成列向量,避免出现整 5 秒点多采一个样本;初始相位 0.3 rad 是为了让信号不是标准余弦起点,模拟真实测量的非对齐状态。

3.2 时域 RMS 与频域 PSD 积分的 GRMS 校核

GRMS1_GRMS2.m这个名字的含义很直白,脚本里给出了两种 GRMS 计算路径。一条路径是直接在时域对加速度样本求均方根sqrt(mean(x.^2));另一条路径是先把加速度信号变换到频域得到单边 PSD,再数值积分开方sqrt(sum(psd(2:end)) * df)。第二条路径里从第二条谱线开始累加,目的是排除直流分量,传感器偏置、零漂这类直流成分不属于振动能量。

% 两种 GRMS 计算路径对比 grms_time = sqrt(mean(x.^2)); % 时域直接 RMS df = Fs / length(x); % 频率分辨率 grms_freq = sqrt(sum(P1(2:end)) * df); % PSD 积分开方 fprintf('时域 GRMS = %.4f\n频域 GRMS = %.4f\n', grms_time, grms_freq);

这里要特别说明加窗带来的能量差异。第 2 章周期图法里数据乘了 hann 窗,时域 RMS 用的是未加窗的原始信号,两者直接对比会差一个窗能量修正系数。hann 窗的功率修正系数约为 0.375,对应 RMS 要除 sqrt(0.375),大约放大 1.63 倍才能和加窗后的频域积分结果对齐。实际工程里我不建议把修正系数叠进报告曲线,更好的做法是保留未加窗的时域 RMS 做基线,再用修正后的 PSD 做频域积分,两张图在数值上差一个恒定比例,用semilogy画出来检查一致性即可。

3.3 正弦扫频与随机振动:两种激励的谱线形态差异

正弦扫频激励是确定性信号,PSD 上表现为尖锐谱峰,峰值高度与采样点数、窗函数直接相关,扫过某一瞬时频率时能量集中在少数几条谱线上。随机振动激励是宽带过程,PSD 曲线平滑,没有明显孤立尖峰。把正弦扫频信号当随机信号做 PSD 估计,会得到峰值高得离谱的曲线;把随机信号按正弦处理去读单根谱线幅值,方差又非常大。Buffet_sin.m里叠加的多个正弦会让 PSD 图出现多个尖锐峰,与宽带随机曲线一眼就能区分,这也提醒使用者:报告 PSD 前先确认激励类型,不同激励对应不同的处理方式和验收标准。

4. 工程实战:pwelch 与周期图法对比,采样率、窗长与重叠率怎么定

4.1 Welch 法一行代码实现分段平均 PSD

pwelch把周期图法包装成了分段平均流程:将数据切成长度相等且有重叠的段,每段加窗做 FFT,再把所有段的功率谱做平均,方差显著降低,段数越多曲线越平滑。MATLAB 从早期版本到目前主流版本,这个接口的参数形式基本稳定,实践里可以直接照下面格式调用。

% Welch 法加速度 PSD 估计 [pxx, f] = pwelch(x, hann(1024, 'periodic'), 512, 1024, Fs); % 参数依次为:信号、窗函数、重叠点数、FFT 点数、采样率

参数说明:窗长 1024 对应频率分辨率约 Fs/1024 = 0.98 Hz;重叠 512 即 50% 重叠,数据利用率高,相邻分段相关性适中;FFT 点数取 1024 与窗长一致,补零只能做频域插值,不会提高真实频率分辨能力,不必为了曲线更细而盲目加大 nfft。该接口返回的单边 PSD 已经完成了负频率能量合并,单位是工程单位平方/Hz,与手写周期图法的结果可直接对比。

4.2 重叠率、窗长与平滑度的权衡参数表

参数常用值对结果的影响
窗长512 / 1024 / 2048窗越长 df 越小,但分段内非平稳成分会被平均掉
重叠率50% / 75%重叠越高方差越低,对非平稳信号反而有害
FFT 点数≥ 窗长大于窗长只是频域插值,不改变频率分辨能力
分段数越多越好方差下降,但信号局部特征被抹平

随机振动数据处理里我默认 50% 重叠加 hann 窗,先跑一版看曲线毛刺。毛刺太多就把重叠提到 75%,分段数几乎翻倍,方差大约再降一半;如果原本期望看到的窄带峰被抹平了,说明窗长太长或重叠太高,把窗长减半重算。非平稳数据比如扫频试验或转速爬升过程,尽量不要用高重叠,否则时间分辨率丢失,频带变化会被平均成模糊一片。

4.3 从 PSD 曲线读问题:尖峰、平带与噪声本底

实测 PSD 曲线有三种典型形态:孤立尖峰表示存在周期性分量,比如电机转频、齿轮啮合频率;平坦宽带区表示随机激励主导;高频段持续衰减后趋于本底,说明传感器噪声或抗混叠滤波器在起作用。做故障诊断时先定位尖峰频率,再换算成转速、叶片通过频率等物理量;做环境试验时更关心整段曲线是否落在试验大纲容差带内。不要一看到尖峰就判定为故障,先检查是否来自供电工频或结构共振;对较窄的峰可以放大局部频率轴确认谱线宽度,单一频点尖峰与展宽峰对应完全不同的激励源。

5. 进阶校核:用 PSD 面积反推 GRMS 并排查估算误差

5.1 从 pwelch 结果精确反推 GRMS

Welch 法返回的频率轴可能不是严格等间隔的吗?实际是严格等间隔的,但频率间隔与窗长、nfft 的关系容易被写错,所以积分时从返回的频率轴直接取 df 最稳妥。

% 从 pwelch 结果反推 GRMS 校核 df = f(2) - f(1); % 频率间隔从实际频率轴取 grms_est = sqrt(sum(pxx(2:end)) * df); % 去掉直流分量再积分

这里用f(2)-f(1)而不是直接用 Fs/N,是因为 Welch 法里频率轴由窗长和 nfft 共同决定,手写 Fs/N 容易在窗长不等于 nfft 时算错。pxx(2:end)排除直流谱线,机械振动分析中直流分量通常视为测量偏置,不应计入振动总能量。

5.2 常见误用自查清单

  • 单边谱未乘 2:低频段能量偏小,GRMS 偏小约根号 2 倍。
  • 加窗 PSD 与未加窗时域 RMS 直接对比:hann 窗能量修正约 1.63 倍,比较前需要换算。
  • 把 PSD 的纵轴当幅值谱读:PSD 是功率密度,单位是工程单位平方/Hz,与 FFT 幅值谱不同,不能混读。
  • 用 Fs/N 计算 Welch 的 df 而不是取返回频率轴:窗长与 nfft 不一致时结果错位。
  • 未去趋势就做积分:传感器零漂会让低频段能量异常抬高,GRMS 虚大。

5.3 csv 数据导入与 fft 仿真衔接

实测加速度数据经常以 csv 格式保存,导入后用同样流程做 fft 仿真:

data = readmatrix('acc_data.csv'); % 读取 csv x = data(:, 2); % 取加速度通道 x = x - mean(x); % 去直流,消除零漂

若 csv 带表头,用readmatrix('acc_data.csv', 'NumHeaderLines', 1);单位若是 mg,按 1 g = 1000 mg 换算后再计算 PSD,否则 GRMS 数值与振动台大纲对不上。做完这些校核后,我习惯在图注里多标一行窗类型、窗长、重叠率和 df,方便隔几个月回读数据时直接知道这张 PSD 曲线是怎么统计出来的。

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

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询