简介面向需要开展加速度信号频域分析的研究生、工程师与MATLAB初学者这套实例演示了从正弦波模拟到加速度功率谱密度PSD计算的完整链路。压缩包共2个文件全部为.m源码整体仅1KB代码精简却覆盖了信号处理的核心环节适合逐行研读与二次开发。实例中利用sin函数构造模拟加速度信号并在时间域完成数据预处理随后借助快速傅里叶变换FFT得到频谱再通过模平方归一化获得功率谱估计同时构建对应频率轴最终以对数坐标绘制PSD曲线便于观察谐波与噪声特征。此外代码中还给出了GRMS均方根加速度的计算逻辑可直接用于振动台数据、机械故障诊断或地震监测等场景。目前已有284人学习下载对希望快速掌握PSD算法实现、理解加速度信号频域特性的读者而言是一份简洁实用的参考脚本。1. 加速度PSD分析在MATLAB里到底在算什么振动台做耐久试验时传感器贴在试件上采集到的是一串加速度时间序列单位是 m/s²。直接看时域波形只能判断幅值大小回答不了“能量集中在哪个频率”这个问题。加速度功率谱密度acceleration PSD把时域信号变换到频域描述单位频率宽度内的加速度能量分布单位是 (m/s²)²/Hz。NVH工程师看它判断共振峰可靠性工程师拿它做随机振动疲劳输入控制工程师用它验证激振谱形。MATLAB 里最常用的计算函数是pwelch基于 Welch 平均周期图法加窗、分段、FFT、平均一条龙。适合手里已有加速度数据、想快速算出 PSD 曲线但不想从零手写离散傅里叶变换的人。这篇不绕理论推导直接讲参数怎么设、代码怎么写、结果怎么验证。2. APSD 核心参数窗函数、重叠率与分辨率设置加速度信号跟一般电压信号有个明显区别动态范围大冲击成分多而且工程上关注的多是宽频随机成分叠加若干窄带共振峰。如果直接用整段数据做一次 FFT频谱泄漏和噪声方差会让人没法看。Welch 方法的核心思想是把长信号切成多段每段加窗后做 FFT再把所有段的功率谱平均。这个流程里有四个参数决定最终曲线质量窗函数类型、每段长度、重叠率、FFT 点数。2.1 窗函数为什么是 PSD 的第一步FFT 隐含假设信号是周期延拓的而实际信号的段边界几乎一定不连续截断造成频谱泄漏。窗函数的作用是让段两端平滑衰减到零减小边界突变。加速度信号分析中Hanning 窗是默认选择主瓣宽度适中旁瓣衰减 31 dB对连续谱的泄漏抑制效果远好于矩形窗。只有当信号本身是瞬态冲击比如锤击试验的脉冲响应才优先考虑矩形窗或力窗避免把有效能量也衰减掉。fs 2048; % 采样率单位 Hz t 0:1/fs:10-1/fs; % 10 秒加速度数据 accel 0.5*sin(2*pi*50*t) 1.2*sin(2*pi*120*t) 0.1*randn(1, length(t)); win hanning(2048); % 2048 点 Hanning 窗对应 1 秒数据 noverlap 1024; % 50% 重叠 nfft 2048; % FFT 点数与窗长一致 [pxx, f] pwelch(accel, win, noverlap, nfft, fs);win向量直接传给pwelch函数内部会自动把数据按窗长分段。Hanning 窗在这里的作用是把每段首尾压到接近零代价是主瓣变宽也就是频率分辨率会下降。如果信号里有两个频率很近的共振峰比如 49 Hz 和 51 Hz窗长 2048 点、采样率 2048 Hz 时的频率分辨率只有 1 Hz两个峰基本分不开这时候需要加长窗段。2.2 重叠率如何影响方差与自由度加窗之后段首尾的数据被压小信息利用率降低所以相邻段之间通常有一段重叠。50% 重叠是工程上性价比最高的选择除了第一段和最后一段每个数据点几乎被用了两次而段之间的相关性还不至于过高。75% 重叠能进一步降低谱估计方差但计算量翻倍而且相邻段高度相关平均带来的方差改善远低于理论值。随机振动测试标准里经常要求 50% 或 75% 重叠测试目的不同选择不同。重叠率的物理意义体现在平均次数上。10 秒数据、窗长 1 秒、50% 重叠实际分段数是 19 段功率谱在每根频率线上做了 19 次平均。谱估计的相对标准差大约是 1/sqrt(平均次数)19 次平均对应约 23% 的标准差。如果换成 75% 重叠分段数变成 37标准差降到 16%曲线更平滑但谱线之间的独立性变差。对加速度 PSD 来说随机振动背景上叠加的共振峰是确定性的不会因为平均次数少而消失所以追求的是背景噪声平滑不是把峰抹平。2.3 分辨率与频率轴的换算关系频率分辨率等于 fs 除以窗长点数。窗长 2048 点、采样率 2048 Hz分辨率是 1 Hz。工程上常见的做法是先确定目标分辨率再反推窗长。比如随机振动控制标准里常用 2 Hz 分辨率采样率 2048 Hz 时窗长取 1024 点就够。nfft可以比窗长大多出来的部分会在 FFT 前自动补零补零只让谱线变密不会提高真实分辨率因为分辨率是窗长决定的不是 FFT 点数决定的。参数典型值对结果的影响组合建议窗函数hanning泄漏抑制与主瓣宽度折中连续随机信号默认窗长1024~4096 点决定频率分辨率分辨率 fs / 窗长重叠率50%~75%影响平滑度和计算量默认 50%数据短取 75%nfft≥ 窗长只加密谱线不提高分辨率设为与窗长相同的 2 的幂采样率由硬件决定决定奈奎斯特频率分析上限 fs/2pwelch返回的f轴范围是 0 到 fs/2对应单边谱。加速度 PSD 的纵轴习惯上用 dB 表示10*log10(pxx)转成 dB 后画图动态范围大的信号更容易看出趋势。写代码时要注意pwelch默认输出是功率谱密度不是功率谱两者的区别在于前者除以了频率分辨率单位是每赫兹的功率这一点跟直接对 FFT 幅值取平方是完全不同的。用pwelch而不是abs(fft(x)).^2的原因就在这里后者既没做平均也没考虑窗函数能量增益算出来的谱噪声极重。3. pwelch 实现加速度 PSD 的完整 MATLAB 过程从一段加速度时间序列到一张可用的 PSD 曲线中间有个固定套路数据清洗、去趋势、定参数、算谱、验证。很多初学者直接拿原始信号算 PSD结果在低频段看到巨大的能量斜坡其实是传感器零漂或者温度漂移造成的趋势项。3.1 从加速度时间序列到 PSD 的最小代码数据清洗的第一步是看时域波形有没有异常尖峰比如线缆松动瞬间的脉冲。这类脉冲在频域会表现为宽频能量抬高污染整个 PSD。常见做法是先把明显超出量程或超出均值 5 倍标准差的点找出来做插值替换或者直接删除该段数据。% 数据导入与预处理 data load(accel_run01.txt); t data(:,1); % 时间列 accel_raw data(:,2); % 加速度列单位 m/s^2 % 去除趋势项消除传感器零漂 accel detrend(accel_raw, linear); % 剔除异常尖峰超过均值±5倍标准差则替换为局部中值 mu mean(accel); sd std(accel); idx abs(accel - mu) 5*sd; accel(idx) movmedian(accel, 101, omitnan);detrend函数默认拟合一次多项式并减去对加速度信号来说线性趋势是最常见的干扰成分。movmedian的窗口长度设成 101 个采样点约 50 ms对于 2048 Hz 采样率来说既不会抹掉真实的冲击特征又能把孤立尖峰压掉。在随机振动试验里5 倍标准差这个阈值基本不会误伤正常数据因为随机信号本身就是高斯分布的超过 5 倍标准差的事件概率极低。清洗完成后进入核心计算参数沿用上一节的原则。这里把窗长、重叠、nfft 三个参数用变量提出来方便后续批量调整。% PSD 计算参数 seg_len 2048; % 每段 1 秒 noverlap round(seg_len * 0.5); % 50% 重叠 nfft seg_len; % FFT 点数等于窗长 window hanning(seg_len); % 计算单边加速度功率谱密度 [pxx, f] pwelch(accel, window, noverlap, nfft, fs); % 转 dB 并绘图 figure; semilogx(f, 10*log10(pxx), LineWidth, 1.2); xlabel(频率 (Hz)); ylabel(加速度 PSD (dB re 1 (m/s^2)^2/Hz)); grid on; xlim([1 fs/2]);这段代码里pwelch的第五个参数fs必须传否则输出的f轴是归一化频率范围 0 到 1对工程分析没有直接意义。semilogx画对数频率轴是因为机械振动分析关心的是低频段共振峰线性频率轴会把 1~100 Hz 的细节全挤在左边。纵轴单位默认是线性功率密度转成 dB 后更容易读峰谷差异。3.2 加窗、重叠与 FFT 点数的组合策略三个参数不是孤立的。窗长决定了重叠的数量重叠率又反过来说明了有效平均次数。MATLAB 里有一个被忽略的函数spectrogram本质就是逐段做 STFT取幅值平方再平均跟pwelch是同一套底层逻辑区别在于pwelch默认对每个频点做了窗函数的能量归一化。窗函数能量归一化这一点很容易踩坑。pwelch内部会用窗函数的均方值做归一化确保白噪声输入下测得的 PSD 幅值是准的。如果你手动实现 Welch 算法忘了除以sum(window.^2)/seg_len结果会整体偏小偏小的倍率正好是窗函数的相干增益。这也是为什么我建议直接用pwelch而不是自己写循环。对于短数据比如只有 2 秒的加速度信号分段的矛盾就出来了。想要 1 Hz 分辨率就得用 2048 点窗长2048 Hz 采样率下正好 1 秒2 秒数据最多分 2 段50% 重叠也只能出来 3 段平均。这种场景下再纠结平滑度没有意义只能接受曲线粗糙的现实。更实际的做法是降低分辨率要求窗长缩到 1024 点分辨率放宽到 2 Hz段数翻倍谱线方差改善明显。3.3 单边谱与物理单位换算pwelch默认输出单边谱直流分量和奈奎斯特频率处的能量只算一次其余频点的能量乘了 2把负频率部分折叠过来。这是 PSD 分析的标准做法。加速度 PSD 单位是 (m/s²)²/Hz如果想换算成 g²/Hz1 g 等于 9.80665 m/s²直接除以 9.80665 的平方。实际项目里还有一类需求是把加速度 PSD 积分成速度或位移 PSD。加速度除以 (2πf) 的平方得到位移但这里要注意频率 f 是从 0 开始的低频段会出现除法爆炸所以积分前必须设置下限频率。% 加速度 PSD 转速度 PSD按频点逐点计算 f_min 1; % 积分下限频率单位 Hz valid f f_min; psd_v zeros(size(pxx)); psd_v(valid) pxx(valid) ./ (2*pi*f(valid)).^2;速度 PSD 单位变成 (m/s)²/Hz。这个换算常用于计算振动量级的 RMS 值对 PSD 曲线做数值积分后开平方得到的就是该频带内的加速度 RMS。RMS 值是试验控制里最基础的指标写着 3.3 g RMS 的试验条件本质就是对 PSD 曲线从 20 Hz 积分到 2000 Hz 再开根号。4. 正弦激励信号 PSD 分析的不同场景标题里的 sin 指的是正弦信号。实验室里正弦扫频和正弦定频试验是除随机振动外最常见的两种激励方式它们的 PSD 分析目标与随机振动完全不同。随机振动关心频带内的能量分布正弦试验关心的是某个特定频率上的幅值大小和谐波失真情况。4.1 正弦信号的 PSD 理论值与实测差异理想单频正弦信号 x(t) A·sin(2πf₀t) 的功率集中在 f₀ 一个频点上其功率为 A²/2。PSD 是功率对频率的密度函数理想情况下 f₀ 处是一个无穷窄的谱线密度值趋向无穷大。但实际 FFT 是离散的能量会泄漏到相邻频点所以实测 PSD 在 f₀ 附近形成一个尖峰峰的高度跟窗函数和频率分辨率有关。窗长 1 秒、矩形窗时PSD 峰值的宽度就是 1 Hz 的谱线间隔峰高约等于 A²/2 除以频率分辨率。换成 Hanning 窗主瓣变宽到 2 Hz峰高减半但旁瓣衰减快得多。这给了一个重要结论正弦信号的 PSD 峰值不是绝对的物理量它受参数影响很大。报告里写正弦试验结果时必须同时列出窗函数和分析带宽否则数值没有可比性。% 正弦定频信号 PSD 分析 fs 8192; t 0:1/fs:5-1/fs; f0 100; % 激励频率 100 Hz accel_sin 5.0 * sin(2*pi*f0*t); % 幅值 5 m/s^2 % 参数对比 seg1 4096; [pxx1, f1] pwelch(accel_sin, hanning(seg1), round(seg1*0.5), seg1, fs); % 提取峰值及峰值频率 [pk, loc] max(10*log10(pxx1)); f_peak f1(loc); fprintf(峰值频率 %.2f Hz峰值 %.2f dB\n, f_peak, pk);这个例子输出峰值频率应为 100 Hz这里max函数找的是整条曲线上最高点。如果信号里叠加了噪声峰值位置仍然稳定在 100 Hz不会漂移这是 PSD 相比时域过零检测的一个优势。计算正弦信号 RMS 应该用幅值除以 sqrt(2)不要从 PSD 积分反推因为峰值受窗函数影响反推回去误差很大。4.2 泄漏与窗函数在正弦场景下的表现正弦信号是典型的窄带信号对泄漏比随机信号敏感得多。频率 f₀ 正好落在 FFT 谱线上的整数倍位置时矩形窗也能给出准的峰值。f₀ 偏离谱线半个格时矩形窗的泄漏可以把主瓣能量分散到十几个频点峰值掉几 dB。Hanning 窗对频率偏移的容忍度好很多即使 f₀ 落在两跟谱线正中间峰值损失也在 1.5 dB 以内。所以正弦 PSD 分析的默认窗是 Hanning除非你明确知道 f₀ 与频谱网格对齐且需要最高的频率精度。某些专业模态分析软件用 flattop 窗它的主瓣更宽平幅值精度极高代价是频率分辨率更差。正弦试验的幅值精度比频率精度更敏感时可以用 flattop 窗验证一次。多正弦叠加的工况下两个频率间隔小于分辨率的一半就无法分辨了。比如 100 Hz 和 101 Hz 同时激励需要 0.5 Hz 分辨率也就是 2 秒的窗长。数据长度够的情况下加大窗长比所有窗函数都管用。5. 最小可复现代码与验证方法把前面的内容压缩成一个可以直接运行的脚本附带验证逻辑方便换自己的数据时快速上手。5.1 用合成信号验证算法正确性验证 PSD 算法是否写对不依赖真实数据也能做。构造一段已知功率的白噪声和已知幅值的正弦信号算完 PSD 以后和理论值对比。白噪声用randn生成方差设成 1理论 PSD 是常数 1单位 (m/s²)²/Hz。实际算出来的 PSD 会在 1 附近波动波动范围 1±0.3 以内说明参数设置基本正确。fs 2048; N fs * 20; % 20 秒数据 rng(42); noise randn(1, N); % 单位白噪声 [pxx_n, f_n] pwelch(noise, hanning(1024), 512, 1024, fs); % 检查频带内数值是否接近 1 mean_val mean(pxx_n(10:end)); fprintf(白噪声 PSD 均值: %.3f (理论值 1.0)\n, mean_val);这里取 10 到末尾频段做平均是为了避开低频段可能存在的窗函数泄漏影响。rng(42)固定随机数种子确保每次运行结果一致方便比对。如果平均值偏离 1 超过 50%优先检查采样率是不是填错了其次是窗长和 nfft 是否一致。5.2 能级检查与频段积分真实数据的 PSD 是否正确可以用能级积分验证。计算 PSD 曲线在 10 Hz 到 100 Hz 频段内的积分再开方得到该频段 RMS 值。用滤波后的时域信号直接算 RMS 做对比两者应该在 3% 以内吻合。% 频段 RMS 对比验证 band f 10 f 100; rms_psd sqrt(trapz(f(band), pxx(band))); % 时域直接法 [b, a] butter(4, [10 100]/(fs/2), bandpass); filtered filtfilt(b, a, accel); rms_time rms(filtered); fprintf(PSD积分RMS: %.4f时域RMS: %.4f\n, rms_psd, rms_time);两组数值相对误差在 3% 以内说明 PSD 没有量纲错误。误差偏大时检查pwelch是否用了单边谱如果误用了双边谱且没有乘以 2积分结果会瘦一半。trapz是梯形积分处理不等间距谱线时也用线性插值这里频率轴是均匀的直接传整个频段索引即可。最后提醒一个常见坑pwelch传多个数据列时会把每列独立计算返回的pxx尺寸是频率点数乘列数。列数据如果是一列时间戳一列加速度直接传给pwelch会把时间戳当加速度也算一遍谱结果曲线完全不对。先做pwelch(accel_signal, ...)而不是整个数据矩阵。本文还有配套的精品资源点击获取
简介面向需要开展加速度信号频域分析的研究生、工程师与MATLAB初学者这套实例演示了从正弦波模拟到加速度功率谱密度PSD计算的完整链路。压缩包共2个文件全部为.m源码整体仅1KB代码精简却覆盖了信号处理的核心环节适合逐行研读与二次开发。实例中利用sin函数构造模拟加速度信号并在时间域完成数据预处理随后借助快速傅里叶变换FFT得到频谱再通过模平方归一化获得功率谱估计同时构建对应频率轴最终以对数坐标绘制PSD曲线便于观察谐波与噪声特征。此外代码中还给出了GRMS均方根加速度的计算逻辑可直接用于振动台数据、机械故障诊断或地震监测等场景。目前已有284人学习下载对希望快速掌握PSD算法实现、理解加速度信号频域特性的读者而言是一份简洁实用的参考脚本。1. 加速度PSD分析在MATLAB里到底在算什么振动台做耐久试验时传感器贴在试件上采集到的是一串加速度时间序列单位是 m/s²。直接看时域波形只能判断幅值大小回答不了“能量集中在哪个频率”这个问题。加速度功率谱密度acceleration PSD把时域信号变换到频域描述单位频率宽度内的加速度能量分布单位是 (m/s²)²/Hz。NVH工程师看它判断共振峰可靠性工程师拿它做随机振动疲劳输入控制工程师用它验证激振谱形。MATLAB 里最常用的计算函数是pwelch基于 Welch 平均周期图法加窗、分段、FFT、平均一条龙。适合手里已有加速度数据、想快速算出 PSD 曲线但不想从零手写离散傅里叶变换的人。这篇不绕理论推导直接讲参数怎么设、代码怎么写、结果怎么验证。2. APSD 核心参数窗函数、重叠率与分辨率设置加速度信号跟一般电压信号有个明显区别动态范围大冲击成分多而且工程上关注的多是宽频随机成分叠加若干窄带共振峰。如果直接用整段数据做一次 FFT频谱泄漏和噪声方差会让人没法看。Welch 方法的核心思想是把长信号切成多段每段加窗后做 FFT再把所有段的功率谱平均。这个流程里有四个参数决定最终曲线质量窗函数类型、每段长度、重叠率、FFT 点数。2.1 窗函数为什么是 PSD 的第一步FFT 隐含假设信号是周期延拓的而实际信号的段边界几乎一定不连续截断造成频谱泄漏。窗函数的作用是让段两端平滑衰减到零减小边界突变。加速度信号分析中Hanning 窗是默认选择主瓣宽度适中旁瓣衰减 31 dB对连续谱的泄漏抑制效果远好于矩形窗。只有当信号本身是瞬态冲击比如锤击试验的脉冲响应才优先考虑矩形窗或力窗避免把有效能量也衰减掉。fs 2048; % 采样率单位 Hz t 0:1/fs:10-1/fs; % 10 秒加速度数据 accel 0.5*sin(2*pi*50*t) 1.2*sin(2*pi*120*t) 0.1*randn(1, length(t)); win hanning(2048); % 2048 点 Hanning 窗对应 1 秒数据 noverlap 1024; % 50% 重叠 nfft 2048; % FFT 点数与窗长一致 [pxx, f] pwelch(accel, win, noverlap, nfft, fs);win向量直接传给pwelch函数内部会自动把数据按窗长分段。Hanning 窗在这里的作用是把每段首尾压到接近零代价是主瓣变宽也就是频率分辨率会下降。如果信号里有两个频率很近的共振峰比如 49 Hz 和 51 Hz窗长 2048 点、采样率 2048 Hz 时的频率分辨率只有 1 Hz两个峰基本分不开这时候需要加长窗段。2.2 重叠率如何影响方差与自由度加窗之后段首尾的数据被压小信息利用率降低所以相邻段之间通常有一段重叠。50% 重叠是工程上性价比最高的选择除了第一段和最后一段每个数据点几乎被用了两次而段之间的相关性还不至于过高。75% 重叠能进一步降低谱估计方差但计算量翻倍而且相邻段高度相关平均带来的方差改善远低于理论值。随机振动测试标准里经常要求 50% 或 75% 重叠测试目的不同选择不同。重叠率的物理意义体现在平均次数上。10 秒数据、窗长 1 秒、50% 重叠实际分段数是 19 段功率谱在每根频率线上做了 19 次平均。谱估计的相对标准差大约是 1/sqrt(平均次数)19 次平均对应约 23% 的标准差。如果换成 75% 重叠分段数变成 37标准差降到 16%曲线更平滑但谱线之间的独立性变差。对加速度 PSD 来说随机振动背景上叠加的共振峰是确定性的不会因为平均次数少而消失所以追求的是背景噪声平滑不是把峰抹平。2.3 分辨率与频率轴的换算关系频率分辨率等于 fs 除以窗长点数。窗长 2048 点、采样率 2048 Hz分辨率是 1 Hz。工程上常见的做法是先确定目标分辨率再反推窗长。比如随机振动控制标准里常用 2 Hz 分辨率采样率 2048 Hz 时窗长取 1024 点就够。nfft可以比窗长大多出来的部分会在 FFT 前自动补零补零只让谱线变密不会提高真实分辨率因为分辨率是窗长决定的不是 FFT 点数决定的。参数典型值对结果的影响组合建议窗函数hanning泄漏抑制与主瓣宽度折中连续随机信号默认窗长1024~4096 点决定频率分辨率分辨率 fs / 窗长重叠率50%~75%影响平滑度和计算量默认 50%数据短取 75%nfft≥ 窗长只加密谱线不提高分辨率设为与窗长相同的 2 的幂采样率由硬件决定决定奈奎斯特频率分析上限 fs/2pwelch返回的f轴范围是 0 到 fs/2对应单边谱。加速度 PSD 的纵轴习惯上用 dB 表示10*log10(pxx)转成 dB 后画图动态范围大的信号更容易看出趋势。写代码时要注意pwelch默认输出是功率谱密度不是功率谱两者的区别在于前者除以了频率分辨率单位是每赫兹的功率这一点跟直接对 FFT 幅值取平方是完全不同的。用pwelch而不是abs(fft(x)).^2的原因就在这里后者既没做平均也没考虑窗函数能量增益算出来的谱噪声极重。3. pwelch 实现加速度 PSD 的完整 MATLAB 过程从一段加速度时间序列到一张可用的 PSD 曲线中间有个固定套路数据清洗、去趋势、定参数、算谱、验证。很多初学者直接拿原始信号算 PSD结果在低频段看到巨大的能量斜坡其实是传感器零漂或者温度漂移造成的趋势项。3.1 从加速度时间序列到 PSD 的最小代码数据清洗的第一步是看时域波形有没有异常尖峰比如线缆松动瞬间的脉冲。这类脉冲在频域会表现为宽频能量抬高污染整个 PSD。常见做法是先把明显超出量程或超出均值 5 倍标准差的点找出来做插值替换或者直接删除该段数据。% 数据导入与预处理 data load(accel_run01.txt); t data(:,1); % 时间列 accel_raw data(:,2); % 加速度列单位 m/s^2 % 去除趋势项消除传感器零漂 accel detrend(accel_raw, linear); % 剔除异常尖峰超过均值±5倍标准差则替换为局部中值 mu mean(accel); sd std(accel); idx abs(accel - mu) 5*sd; accel(idx) movmedian(accel, 101, omitnan);detrend函数默认拟合一次多项式并减去对加速度信号来说线性趋势是最常见的干扰成分。movmedian的窗口长度设成 101 个采样点约 50 ms对于 2048 Hz 采样率来说既不会抹掉真实的冲击特征又能把孤立尖峰压掉。在随机振动试验里5 倍标准差这个阈值基本不会误伤正常数据因为随机信号本身就是高斯分布的超过 5 倍标准差的事件概率极低。清洗完成后进入核心计算参数沿用上一节的原则。这里把窗长、重叠、nfft 三个参数用变量提出来方便后续批量调整。% PSD 计算参数 seg_len 2048; % 每段 1 秒 noverlap round(seg_len * 0.5); % 50% 重叠 nfft seg_len; % FFT 点数等于窗长 window hanning(seg_len); % 计算单边加速度功率谱密度 [pxx, f] pwelch(accel, window, noverlap, nfft, fs); % 转 dB 并绘图 figure; semilogx(f, 10*log10(pxx), LineWidth, 1.2); xlabel(频率 (Hz)); ylabel(加速度 PSD (dB re 1 (m/s^2)^2/Hz)); grid on; xlim([1 fs/2]);这段代码里pwelch的第五个参数fs必须传否则输出的f轴是归一化频率范围 0 到 1对工程分析没有直接意义。semilogx画对数频率轴是因为机械振动分析关心的是低频段共振峰线性频率轴会把 1~100 Hz 的细节全挤在左边。纵轴单位默认是线性功率密度转成 dB 后更容易读峰谷差异。3.2 加窗、重叠与 FFT 点数的组合策略三个参数不是孤立的。窗长决定了重叠的数量重叠率又反过来说明了有效平均次数。MATLAB 里有一个被忽略的函数spectrogram本质就是逐段做 STFT取幅值平方再平均跟pwelch是同一套底层逻辑区别在于pwelch默认对每个频点做了窗函数的能量归一化。窗函数能量归一化这一点很容易踩坑。pwelch内部会用窗函数的均方值做归一化确保白噪声输入下测得的 PSD 幅值是准的。如果你手动实现 Welch 算法忘了除以sum(window.^2)/seg_len结果会整体偏小偏小的倍率正好是窗函数的相干增益。这也是为什么我建议直接用pwelch而不是自己写循环。对于短数据比如只有 2 秒的加速度信号分段的矛盾就出来了。想要 1 Hz 分辨率就得用 2048 点窗长2048 Hz 采样率下正好 1 秒2 秒数据最多分 2 段50% 重叠也只能出来 3 段平均。这种场景下再纠结平滑度没有意义只能接受曲线粗糙的现实。更实际的做法是降低分辨率要求窗长缩到 1024 点分辨率放宽到 2 Hz段数翻倍谱线方差改善明显。3.3 单边谱与物理单位换算pwelch默认输出单边谱直流分量和奈奎斯特频率处的能量只算一次其余频点的能量乘了 2把负频率部分折叠过来。这是 PSD 分析的标准做法。加速度 PSD 单位是 (m/s²)²/Hz如果想换算成 g²/Hz1 g 等于 9.80665 m/s²直接除以 9.80665 的平方。实际项目里还有一类需求是把加速度 PSD 积分成速度或位移 PSD。加速度除以 (2πf) 的平方得到位移但这里要注意频率 f 是从 0 开始的低频段会出现除法爆炸所以积分前必须设置下限频率。% 加速度 PSD 转速度 PSD按频点逐点计算 f_min 1; % 积分下限频率单位 Hz valid f f_min; psd_v zeros(size(pxx)); psd_v(valid) pxx(valid) ./ (2*pi*f(valid)).^2;速度 PSD 单位变成 (m/s)²/Hz。这个换算常用于计算振动量级的 RMS 值对 PSD 曲线做数值积分后开平方得到的就是该频带内的加速度 RMS。RMS 值是试验控制里最基础的指标写着 3.3 g RMS 的试验条件本质就是对 PSD 曲线从 20 Hz 积分到 2000 Hz 再开根号。4. 正弦激励信号 PSD 分析的不同场景标题里的 sin 指的是正弦信号。实验室里正弦扫频和正弦定频试验是除随机振动外最常见的两种激励方式它们的 PSD 分析目标与随机振动完全不同。随机振动关心频带内的能量分布正弦试验关心的是某个特定频率上的幅值大小和谐波失真情况。4.1 正弦信号的 PSD 理论值与实测差异理想单频正弦信号 x(t) A·sin(2πf₀t) 的功率集中在 f₀ 一个频点上其功率为 A²/2。PSD 是功率对频率的密度函数理想情况下 f₀ 处是一个无穷窄的谱线密度值趋向无穷大。但实际 FFT 是离散的能量会泄漏到相邻频点所以实测 PSD 在 f₀ 附近形成一个尖峰峰的高度跟窗函数和频率分辨率有关。窗长 1 秒、矩形窗时PSD 峰值的宽度就是 1 Hz 的谱线间隔峰高约等于 A²/2 除以频率分辨率。换成 Hanning 窗主瓣变宽到 2 Hz峰高减半但旁瓣衰减快得多。这给了一个重要结论正弦信号的 PSD 峰值不是绝对的物理量它受参数影响很大。报告里写正弦试验结果时必须同时列出窗函数和分析带宽否则数值没有可比性。% 正弦定频信号 PSD 分析 fs 8192; t 0:1/fs:5-1/fs; f0 100; % 激励频率 100 Hz accel_sin 5.0 * sin(2*pi*f0*t); % 幅值 5 m/s^2 % 参数对比 seg1 4096; [pxx1, f1] pwelch(accel_sin, hanning(seg1), round(seg1*0.5), seg1, fs); % 提取峰值及峰值频率 [pk, loc] max(10*log10(pxx1)); f_peak f1(loc); fprintf(峰值频率 %.2f Hz峰值 %.2f dB\n, f_peak, pk);这个例子输出峰值频率应为 100 Hz这里max函数找的是整条曲线上最高点。如果信号里叠加了噪声峰值位置仍然稳定在 100 Hz不会漂移这是 PSD 相比时域过零检测的一个优势。计算正弦信号 RMS 应该用幅值除以 sqrt(2)不要从 PSD 积分反推因为峰值受窗函数影响反推回去误差很大。4.2 泄漏与窗函数在正弦场景下的表现正弦信号是典型的窄带信号对泄漏比随机信号敏感得多。频率 f₀ 正好落在 FFT 谱线上的整数倍位置时矩形窗也能给出准的峰值。f₀ 偏离谱线半个格时矩形窗的泄漏可以把主瓣能量分散到十几个频点峰值掉几 dB。Hanning 窗对频率偏移的容忍度好很多即使 f₀ 落在两跟谱线正中间峰值损失也在 1.5 dB 以内。所以正弦 PSD 分析的默认窗是 Hanning除非你明确知道 f₀ 与频谱网格对齐且需要最高的频率精度。某些专业模态分析软件用 flattop 窗它的主瓣更宽平幅值精度极高代价是频率分辨率更差。正弦试验的幅值精度比频率精度更敏感时可以用 flattop 窗验证一次。多正弦叠加的工况下两个频率间隔小于分辨率的一半就无法分辨了。比如 100 Hz 和 101 Hz 同时激励需要 0.5 Hz 分辨率也就是 2 秒的窗长。数据长度够的情况下加大窗长比所有窗函数都管用。5. 最小可复现代码与验证方法把前面的内容压缩成一个可以直接运行的脚本附带验证逻辑方便换自己的数据时快速上手。5.1 用合成信号验证算法正确性验证 PSD 算法是否写对不依赖真实数据也能做。构造一段已知功率的白噪声和已知幅值的正弦信号算完 PSD 以后和理论值对比。白噪声用randn生成方差设成 1理论 PSD 是常数 1单位 (m/s²)²/Hz。实际算出来的 PSD 会在 1 附近波动波动范围 1±0.3 以内说明参数设置基本正确。fs 2048; N fs * 20; % 20 秒数据 rng(42); noise randn(1, N); % 单位白噪声 [pxx_n, f_n] pwelch(noise, hanning(1024), 512, 1024, fs); % 检查频带内数值是否接近 1 mean_val mean(pxx_n(10:end)); fprintf(白噪声 PSD 均值: %.3f (理论值 1.0)\n, mean_val);这里取 10 到末尾频段做平均是为了避开低频段可能存在的窗函数泄漏影响。rng(42)固定随机数种子确保每次运行结果一致方便比对。如果平均值偏离 1 超过 50%优先检查采样率是不是填错了其次是窗长和 nfft 是否一致。5.2 能级检查与频段积分真实数据的 PSD 是否正确可以用能级积分验证。计算 PSD 曲线在 10 Hz 到 100 Hz 频段内的积分再开方得到该频段 RMS 值。用滤波后的时域信号直接算 RMS 做对比两者应该在 3% 以内吻合。% 频段 RMS 对比验证 band f 10 f 100; rms_psd sqrt(trapz(f(band), pxx(band))); % 时域直接法 [b, a] butter(4, [10 100]/(fs/2), bandpass); filtered filtfilt(b, a, accel); rms_time rms(filtered); fprintf(PSD积分RMS: %.4f时域RMS: %.4f\n, rms_psd, rms_time);两组数值相对误差在 3% 以内说明 PSD 没有量纲错误。误差偏大时检查pwelch是否用了单边谱如果误用了双边谱且没有乘以 2积分结果会瘦一半。trapz是梯形积分处理不等间距谱线时也用线性插值这里频率轴是均匀的直接传整个频段索引即可。最后提醒一个常见坑pwelch传多个数据列时会把每列独立计算返回的pxx尺寸是频率点数乘列数。列数据如果是一列时间戳一列加速度直接传给pwelch会把时间戳当加速度也算一遍谱结果曲线完全不对。先做pwelch(accel_signal, ...)而不是整个数据矩阵。本文还有配套的精品资源点击获取
