
简介面向需要使用MATLAB开展振动信号分析与功率谱密度计算的工程技术人员与信号处理学习者压缩包提供了一个可直接运行的算法实例重点解决加速度信号的sin波模拟与PSD求解问题。资源覆盖从信号生成、窗函数处理、FFT变换到功率谱归一化、对数展示的完整流程适用于机械振动监测、地震信号分析等场景。包内共2个m文件整体约1KB代码简洁便于快速阅读和二次修改。已有284人学习下载适合希望用短小实例快速上手MATLAB加速度PSD分析的用户。通过运行这两个脚本能够直观掌握基于fft的加速度信号功率谱估计方法并将频率轴构建、功率谱归一化等关键步骤迁移到自身数据中。资源虽小但完整呈现了从时域建模仿真到频域特征提取的链路对提升信号处理理论与工程实践结合能力具有不错的参考价值。1. 为什么加速度信号的 PSD 要这样算从 GRMS 指标反推需求做振动试验的人经常会在试验大纲里看到 GRMS 这个指标它是随机振动中总均方根加速度单位是 g 或 m/s²。拿到一套 MATLAB 的加速度功率谱密度算例里面只有GRMS1_GRMS2.m和Buffet_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再乘频率分辨率 dfFs/N两个写法数值完全等价区别只是归一化放在前面还是后面。容易翻车的是单边谱合并FFT 输出以 Fs/2 对称绘图时通常只保留 0 到 Fs/2 的 N/21 条谱线除 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/21); % 取单边谱 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 法一行代码实现分段平均 PSDpwelch把周期图法包装成了分段平均流程将数据切成长度相等且有重叠的段每段加窗做 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 结果精确反推 GRMSWelch 法返回的频率轴可能不是严格等间隔的吗实际是严格等间隔的但频率间隔与窗长、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 曲线是怎么统计出来的。本文还有配套的精品资源点击获取