简介一个面向MATLAB信号处理与振动分析场景的算法示例聚焦加速度信号模拟与功率谱密度PSD求解压缩包共2个M文件大小仅1KB包含从正弦波生成、加窗预处理、FFT频谱变换到PSD归一化计算的完整代码流程。适合分析振动、冲击或噪声数据的工程师、科研人员也适合准备深入学习频域方法的MATLAB用户已有284人学习下载。学习该示例可掌握加速度信号仿真与PSD计算核心步骤构造时间轴与正弦信号、选择窗函数抑制频谱泄漏、对FFT结果归一化并绘制对数功率谱同时代码中引入GRMS指标计算有助于理解加速度功率谱在工程振动评估、地震监测与机械故障诊断中的实际应用。两个M文件分别承担信号构造与统计算法演示注释清晰、逻辑连贯便于按行阅读和二次修改通过调整频率、采样率等参数还能观察不同信号对PSD估计的影响深入理解窗函数、幅值归一化和频率分辨率等概念。1. 从加速度信号到功率谱密度为什么振动分析绕不开 PSD拿到一段加速度计采集的时域信号直接做 FFT 看幅值谱是大多数 MATLAB 新手的第一反应。但放到实际工程里无论是汽车 NVH、机床主轴振动监测还是桥梁健康检测你会发现正规报告里几乎只用功率谱密度PSDPower Spectral Density很少有人直接贴幅值谱。原因很直接幅值谱告诉你某个频率上的振动有多强而 PSD 告诉你这个频率附近单位频带内携带了多少功率两者差了窗函数、分辨率带宽和噪声带宽的换算关系。对于随机振动或叠加了噪声的加速度信号PSD 才是统计学上稳定、可对比、可直接换算到物理单位如 g²/Hz 或 (m/s²)²/Hz的指标。本篇文章围绕加速度信号加速度功率谱这条主线讲清楚从采集数据到一份可用 PSD 曲线的最小实现路径重点落在 MATLAB 的 pwelch 函数和自功率谱估计Auto PSD即标题中的 apsd上。标题里提到的 sin 算法常见的解读有两种一是激励源为正弦扫频信号二是想表达单频正弦信号的 PSD 估计。本文按后者展开顺带覆盖前者。适合正在做振动测试数据处理、想搞清楚 PSD 和 FFT 区别、以及需要把 MATLAB 的 PSD 结果换算成工程单位的工程师。全文配套代码以 MATLAB 2023a 为基础编写低版本 R2017b 以上都能直接运行。2. PSD 估计的数学基础与 MATLAB 里的 apsd 实现路径2.1 自功率谱与互功率谱apsd 到底在算什么功率谱密度估计分两类自功率谱密度Auto PSD和互功率谱密度Cross PSD。标题中的 apsd 几乎可以确定是 Auto PSD 的缩写。自功率谱描述的是单个信号自身能量在频域的分布而互功率谱描述两个信号在频域上的相关程度。加速度信号分析中我们关心的是某个测点的振动能量集中在哪些频段自然用的是自功率谱。从定义上看自功率谱是信号自相关函数的傅里叶变换。但工程上没人真去先算自相关再做变换而是用周期图法Periodogram或其改进版本。周期图法的核心计算式是PSD(f) |X(f)|² / (fs * N * 窗能量修正)其中 X(f) 是加窗后信号的 FFT 结果fs 是采样率N 是 FFT 点数。分母里的 fs 把功率从每比特归一化到每赫兹这就是 PSD 和幅值谱最本质的区别——PSD 是密度不是幅值。如果输入的时域信号 x(t) 单位是 m/s²那么计算出的 PSD 单位就是 (m/s²)²/Hz。用重力加速度 g9.80665 m/s²归一化后常见的工程单位是 g²/Hz。很多人拿到 MATLAB 的 pwelch 输出直接画图纵轴数值小到 1e-4 量级然后怀疑代码写错了——其实只是单位是 (m/s²)²/Hz 而已换算成 g²/Hz 需要除以 96.17即 9.80665²。2.2 Welch 法为什么是加速度 PSD 的默认选择MATLAB 中计算 PSD 的函数有好几个periodogram、pwelch、cpsd、mscohere 等。其中 pwelch 用的是 Welch 重叠段平均法原理是把长信号分成若干段、每段加窗、做 FFT、求功率、再对多段结果取平均。这样做有两个直接好处。第一是方差降低。单段周期图的方差很大谱曲线毛刺多不稳定而 Welch 法对 L 段独立数据取平均方差近似降为原来的 1/L。加速度信号往往长达几十秒甚至几分钟完全不缺分段的数据量为什么不利用这个统计优势。第二是控制频谱泄漏。直接对整段数据做 FFT矩形窗的旁瓣衰减只有 -13 dB远端泄漏严重换用 Hann 窗后旁瓣衰减能到 -31 dBHanning 加窗配合 50% 重叠是振动测试的默认配置。% 最小化的 Welch PSD 估计示例 fs 2048; % 采样率 2048 Hz t 0:1/fs:10-1/fs; % 10 秒时长 x 0.5*sin(2*pi*50*t) 0.3*sin(2*pi*120*t); % 模拟加速度信号 x x 0.2*randn(size(x)); % 叠加高斯白噪声模拟真实采集 [pxx, f] pwelch(x, hann(1024), 512, 1024, fs); figure; plot(f, 10*log10(pxx)); xlabel(频率 (Hz)); ylabel(功率谱密度 (dB/Hz)); grid on;这段代码中pwelch 的核心参数是窗函数 hann(1024) 表示每段 1024 个点重叠 512 个点50%NFFT 取 1024fs 为 2048。输出 pxx 是单边功率谱密度向量f 是对应的频率向量。10*log10 是为了转为 dB 显示方便同时看清 50 Hz 和 120 Hz 两个谱峰以及噪声基底。如果直接画线性幅值噪声底和谱峰的对比会非常不明显。注意pwelch 默认输出是单边 PSD即频率范围从 0 到奈奎斯特频率 fs/2。如果输入信号是复数才需要改用双边谱加速度信号永远是实数不用考虑这个分支。2.3 从幅值谱换算到 PSD 的手动对照有些场景下项目组还在用 FFT 幅值谱做初判这时可以手动做一次换算来交叉验证 pwelch 的结果避免代码库里有两套口径不一致的处理函数。对长度为 N 的实信号 x加窗后的 FFT 幅值谱为 |X(k)|对应的单边 PSD 为PSD(k) 2 * |X(k)|² / (fs * sum(win.^2))系数 2 是单边谱的功率折叠因子直流分量和奈奎斯特频率处系数为 1但这两点不影响宽带分析分母中的 sum(win.^2) 是窗函数的能量修正。Hann 窗的 sum(win.^2) ≈ 0.375N矩形窗则为 N。这个换算公式在手工核对时非常有用——如果 pwelch 的结果和这个公式算出的数量级对不上基本可以确定是单位或者归一化因子的问题。% 手动换算与 pwelch 结果对照 win hann(1024); N length(win); X fft(x(1:1024) .* win, N); X_p X(1:N/21); psd_manual 2 * abs(X_p).^2 / (fs * sum(win.^2)); psd_manual(1) psd_manual(1) / 2; % DC 点不做功率折叠 psd_manual(end) psd_manual(end) / 2; % 奈奎斯特点不做功率折叠 f_manual (0:N/2) * fs / N; figure; plot(f, 10*log10(pxx), b); hold on; plot(f_manual, 10*log10(psd_manual), r--); legend(pwelch, 手动 FFT 换算);对比结果中两条曲线应该几乎重叠差异仅在浮点精度范围内。这个对照本身就是一个很好的验证手段——它确认了加窗、归一化、单边折叠三个环节都处理正确了。3. 用 pwelch 计算加速度功率谱的 MATLAB 最小实例3.1 从 .mat 或 CSV 导入加速度数据的第一步实际项目中加速度数据往往以 CSV、TXT 或 .mat 文件形式存放。工程上常见的数据格式是第一列时间戳、第二列加速度值单位可能是 g 也可能是 m/s²先确认单位再处理否则后面所有换算都白做。导入用 readmatrix 或 load 都很方便但要注意 readmatrix 对文件头、分隔符的处理有时会静默出错导入后做一个 sanity check 很有必要。% 从 CSV 导入加速度信号并做基本检查 data readmatrix(accel_data.csv); % 两列时间(s), 加速度(m/s²) t_raw data(:, 1); x_raw data(:, 2); fs 1 / mean(diff(t_raw)); % 由时间列反推采样率 fprintf(采样率 %.2f Hz数据长度 %d 点\n, fs, length(x_raw)); figure; plot(t_raw, x_raw); xlabel(时间(s)); ylabel(加速度(m/s²)); title(原始加速度信号);导入后先看时域波形这是最便宜的异常检测手段饱和削顶、零漂、毛刺突变、丢数全都能在时域图上一眼看出。数据无异常再做去趋势因为加速度计本身的零漂会让信号带一个直流偏置直流分量在 PSD 的 0 Hz 处形成一个很大的谱峰会掩盖低频段的真实信息。去趋势操作用 detrend 函数默认去掉均值也就是把零漂移除。如果数据有明显的线性趋势比如加速度计缓慢温漂用 detrend(x_raw, linear)。这里不要用高通滤波器替代去趋势因为滤波器有相位延迟和暂态效应对后续 PSD 估计的频域形状影响比去趋势大得多。3.2 完整可运行的加速度 PSD 分析脚本下面给出一份完整可运行的加速度功率谱分析脚本覆盖从导入到出图、再到峰值频率自动提取的完整链路。这个脚本的定位是骨架代码实际使用时替换数据源、调整参数即可。% 完整加速度 PSD 分析脚本 % 输入加速度时域信号 x采样率 fs % 输出PSD 曲线、峰值频率列表 fs 2048; x detrend(x_raw); % 去除零漂和线性趋势 % 参数配置 segment 2048; % 每段点数对应 1 秒数据 overlap 0.5; % 重叠率 50% nfft 2048; % FFT 点数与段长相同即可 win hann(segment, periodic); % 计算 PSD [pxx, f] pwelch(x, win, round(segment*overlap), nfft, fs); % 转成工程单位从 (m/s²)²/Hz 到 g²/Hz g 9.80665; pxx_g pxx / g^2; % 找出谱峰忽略 DC 附近 1 Hz 以下 f_min_idx find(f 1, 1); [pks, locs] findpeaks(pxx_g(f_min_idx:end), f(f_min_idx:end), ... MinPeakHeight, max(pxx_g)*0.05, MinPeakDistance, 5); figure; semilogx(f(f_min_idx:end), pxx_g(f_min_idx:end)); hold on; plot(locs, pks, rv, MarkerSize, 6); xlabel(频率 (Hz)); ylabel(PSD (g²/Hz)); title(加速度功率谱密度估计); grid on; % 输出峰值频率 for i 1:length(locs) fprintf(谱峰 #%d: %.2f Hz, PSD %.4e g²/Hz\n, i, locs(i), pks(i)); end脚本里几个参数的选取逻辑如下。segment 取 2048 意味着每段数据时长 1 秒频率分辨率为 fs/nfft 1 Hz这正是大多数机械设备振动分析的经验起点。nfft 取与 segment 相同即可取更大的 nfft比如 4096不会提升有效分辨率只会让谱线更密看起来更平滑实际信息量不变。重叠率 0.5 是 Welch 原文的建议值配合 Hann 窗综合来看谱估计偏差和方差平衡最好。findpeaks 的 MinPeakDistance 设为 5 Hz是为了避免同一个谱峰被误报为多个——如果轴系转频和倍频间距小于 5 Hz这个阈值需要调小但那种场景建议先用转速跟踪再分析。3.3 参数选择对 PSD 结果影响的具体对照PSD 估计里最经典的权衡是频率分辨率与方差之间的此消彼长。频率分辨率由 fs/nfft 决定nfft 越大分辨率越高但 nfft 越大意味着分段越少平均次数减少谱曲线方差变大。这个 trade-off 必须通过具体数字理解否则参数就是乱调的。段长点数分辨率 (Hz)重叠 50% 时的平均段数10s 数据谱线平滑度512439好1024219较好204819一般40960.54较差假设采样率 2048 Hz、数据长度 10 秒段长取 512 点时频率分辨率只有 4 Hz如果两个相邻的谱峰比如 50 Hz 和 52 Hz 的边频带相距小于 4 Hz就无法分辨而段长取 4096 时分辨率 0.5 Hz但只分 4 段平均谱线毛刺明显小峰值可能被噪声底淹没。实操建议从段长 1024 起步跑完看一眼曲线如果谱峰太宽分不清边频增加段长如果曲线毛刺太多、谱峰位置跳动减小段长或提高重叠率到 75%。加窗对 PSD 的影响主要体现在谱泄漏上——矩形窗谱峰最尖锐但旁瓣大Hann 窗谱峰略宽但旁瓣小在加速度信号常常包含随机分量分析中 Hann 窗是默认首选不要在没把握的情况下换成 Kaisar 或 Chebyshev 窗除非你已经能解释窗函数的旁瓣衰减指标和等效噪声带宽的含义。4. 加速度信号 PSD 的单位换算、频段积分的实战处理4.1 从 g²/Hz 到 RMS 值的积分换算PSD 曲线本身是密度函数积分才有物理意义。PSD 曲线下的面积等于信号在对应频带内的方差开根号就是 RMS有效值。对于用于振动评估的加速度信号RMS 值是比幅值谱峰值更稳定、更有工程意义的指标。% 计算指定频带内的 RMS 值 f_low 10; % 下限频率 Hz f_high 500; % 上限频率 Hz idx (f f_low) (f f_high); df f(2) - f(1); % 频率分辨率 band_power sum(pxx_g(idx) * df); % 频带内功率面积 band_rms_g sqrt(band_power); fprintf(%d-%d Hz 频带内 RMS 加速度 %.4f g RMS\n, f_low, f_high, band_rms_g);注意这里直接用 sum(pxx * df) 做积分前提是频率轴等间距。pwelch 输出的频率轴在线性坐标下均匀分布0 到 fs/2所以这种累加方式没问题。如果频率轴是对数分布的——因为画图用 semilogx 而误以为数据也是对数分布的——那就不能用这个算法必须先插值到线性轴再积分。这种错误在项目代码里屡见不鲜。另外要强调RMS 积分结果对频带边界极其敏感。如果 f_low9 Hz但频率分辨率只有 1 Hz那 9 Hz 这个索引可能并不存在取整到 f9 或 f10 会直接影响积分结果误差可达 5%-10%。要精确控制积分频带建议用 trapz 做梯形积分它能接受非对齐的积分区间。% 更精确的梯形积分方法 band_rms_g sqrt(trapz(f(idx), pxx_g(idx)));trapz 和 sum(df) 的差别在于梯形积分计入了频带边缘两个半格的影响精度略高。对于惯性导航、军工产品的振动指标验收这个精度很关键。4.2 加速度 PSD 与速度、位移谱的换算加速度 PSD 可以换算成速度 PSD 和位移 PSD。这在工程上特别有用——同一个振动信号加速度谱在高频段突出速度谱强调中频位移谱突出低频三张谱结合起来才能完整描述振动特性。换算公式很简单速度 PSD(f) 加速度 PSD(f) / (2πf)²位移 PSD(f) 加速度 PSD(f) / (2πf)⁴。换算在频域逐点做除法不需要重新采集数据。% 加速度 PSD 换算为速度 PSD 和位移 PSD omega 2*pi*f; psd_v pxx_g .* (g^2) ./ omega.^2; % (m/s)²/Hz psd_d pxx_g .* (g^2) ./ omega.^4; % m²/Hz % 使用对数坐标对比三条曲线 figure; loglog(f(f_min_idx:end), pxx_g(f_min_idx:end), b); hold on; loglog(f(f_min_idx:end), psd_v(f_min_idx:end), r); loglog(f(f_min_idx:end), psd_d(f_min_idx:end), g); legend(加速度 PSD (g²/Hz), 速度 PSD ((m/s)²/Hz), 位移 PSD (m²/Hz)); xlabel(频率 (Hz)); ylabel(PSD);低频段做位移换算时要特别小心f 趋于 0 时除以 f⁴ 会让数值趋近于无穷大任何微小的加速度零漂都在位移谱中被无限放大。实际处理中通常设置一个最低换算频率比如 1 Hz低于该频率的位移谱直接不输出。这也侧面解释了为什么位移谱很少直接从加速度积分得到——时域双重积分的直流漂移问题更难控制。4.3 多段测量数据 PSD 的平均与置信区间评估单次测量的 PSD 方差较大工程上常对同一工况重复测量多次然后对 PSD 取平均同时计算置信区间来评估估计精度。pwelch 函数本身的分段平均是数据内的平均多次测量之间的平均是数据间的平均两者不冲突。% 多段测量数据的 PSD 平均 n_meas 5; % 重复测量次数 psd_all zeros(n_meas, length(f)); for i 1:n_meas % 假设 x_cell{i} 存第 i 次测量的数据fs 相同 [psd_all(i, :), f] pwelch(x_cell{i}, hann(2048,periodic), 1024, 2048, fs); end psd_mean mean(psd_all, 1); % 平均 PSD psd_std std(psd_all, 0, 1); % 标准差 df f(2) - f(1); % 绘制带 ±1σ 阴影带的 PSD figure; semilogx(f, 10*log10(psd_mean), b, LineWidth, 1.5); hold on; fill([f fliplr(f)], ... [10*log10(psd_mean psd_std) fliplr(10*log10(psd_mean - psd_std))], ... b, FaceAlpha, 0.2, EdgeColor, none); xlabel(频率 (Hz)); ylabel(PSD (dB Hz^{-1}));置信带的宽度直接反映谱估计是否稳定。如果 ±1σ 带宽超过 5 dB说明测量次数不够或数据本身非平稳这时候加大平均次数比调窗函数更有效。对数坐标下 PSD 的标准差近似和频率无关所以画出来是一条均匀带宽——如果看到带宽随频率剧变说明信号在某些频段不满足平稳性假设PSD 本身可能不再适用需要改用短时傅里叶变换时频分析。5. 加速度功率谱实测中的典型误区和验证技巧5.1 频谱泄漏为什么 50 Hz 的谱峰会长胖实测数据很少是整周期截断的对非整周期信号直接做 FFT能量就会泄漏到相邻频点上表现为谱峰变宽、基底抬高。Pwelch 的加窗分段已经大幅缓解了这个问题但很多人犯的错误是——只在 pwelch 里指定了窗却忘了检查 nfft 是否覆盖了窗函数的完整长度。如果 nfft 大于窗长MATLAB 会做零填充等效频率分辨率提高但实际物理分辨率不变如果 nfft 小于窗长信号被截断等于隐式换了矩形窗之前选的窗函数白费了。一个快速验证是否存在严重泄漏的方法看 PSD 曲线的谱峰是否呈现钟形平顶的形状。如果 50 Hz 处的谱峰左右呈明显的对称展宽、基底明显抬高基本可以判断分辨率不足或窗函数旁瓣太大。另一种情况是谱峰不对称、一侧有明显拖尾此时怀疑数据里有衰减振荡分量或频率漂移不完全是泄漏问题。更直接的验证是用合成信号校准生成一个幅度已知、频率为 50 Hz 的正弦波加到白噪声上用同样的参数跑一遍 pwelch看 50 Hz 处 PSD 的面积即该频带内的 RMS是否和合成时的幅值对得上。合成了就心里有底因为实际数据里永远不可能告诉你真实值是多少。% 合成信号校准 PSD 流程 fs 2048; t 0:1/fs:30-1/fs; x_syn 0.5*sin(2*pi*50*t) 0.1*randn(size(t)); [pxx_syn, f_syn] pwelch(x_syn, hann(4096,periodic), 2048, 4096, fs); % 积分 48-52 Hz 频带内的功率 idx_syn (f_syn 48) (f_syn 52); rms_est sqrt(trapz(f_syn(idx_syn), pxx_syn(idx_syn))); fprintf(合成 50Hz 正弦 RMS 估计 %.4f理论值 %.4f\n, rms_est, 0.5/sqrt(2));理论 RMS 是 0.5/√2 ≈ 0.3536。如果 rms_est 和这个值偏差超过 2%优先检查窗函数能量归一化是否正确、重叠率是否为 0.5、以及 pwelch 输出是否被误当作双边谱处理。5.2 PSD 结果一致性验证半谱与全谱的对照Matlab 的 pwelch 默认输出单边谱但有些早期代码或从 Python scipy 迁移过来的工程师习惯用双边谱。两边谱的总功率是一样的只是分布在 0 ~ fs/2 还是 0 ~ fs 的区别。如果怀疑手里的脚本把单双边搞混了验证方法非常粗暴用 cumtrapz 从 0 到奈奎斯特频率积分单边 PSD得到全频带 RMS然后把原始时域信号直接算 RMS两者应该非常接近。% 全频带 RMS 对照验证 rms_from_psd sqrt(trapz(f, pxx_g)); % 从 PSD 积分 rms_from_time rms(x / g); % 从时域直接算除以 g 转成 g 单位 fprintf(PSD 积分 RMS %.4f g时域 RMS %.4f g偏差 %.2f%%\n, ... rms_from_psd, rms_from_time, (rms_from_psd-rms_from_time)/rms_from_time*100);偏差超过 3% 时不要犹豫直接从这几个点排查x 是否 detrend 过直流偏置会拉高时域 RMS 但不贡献到 0 Hz 以外的频段、pxx_g 单位换算是否正确少除了 g² 偏差会大到 7 个数量级一眼就能看出来、pwelch 的窗函数是不是用的 periodic 而不是 symmetric 变体。这两个窗变体长度差一个点能量修正差 0.1% 左右虽然不至于导致 3% 的偏差但在高标准计量场景下不可忽略。5.3 把 PSD 结果导出到报告和后续处理的技巧最后落一个实际工作中几乎每次都需要的操作——把 PSD 结果导出为标准 CSV 带表头文件顺便生成一张适合贴进报告的高质量图。很多人用 saveas 直接存 .fig 或者 .png但分辨率一放到 Word 里就糊。更可靠的方案是导出矢量图SVG 或 PDF或者在 MATLAB 里先调整好 Figure 属性再 print 到 300 dpi 的 PNG。% 高质量导出 % 1. 数据导出 out_table table(f(:), pxx_g(:), VariableNames, {Frequency_Hz, PSD_g2_Hz}); writetable(out_table, acceleration_psd_export.csv); % 2. 矢量图导出推荐 SVG figure; semilogx(f, 10*log10(pxx_g), b, LineWidth, 1.2); xlabel(频率 (Hz), FontSize, 11); ylabel(PSD (g²/Hz), FontSize, 11); grid on; set(gcf, PaperPositionMode, auto); exportgraphics(gcf, acceleration_psd_result.svg, ContentType, vector); % 3. 300 dpi 位图导出用于 PPT 快速插入 exportgraphics(gcf, acceleration_psd_result.png, Resolution, 300);exportgraphics 是 MATLAB R2020a 之后的推荐导出函数能正确保留坐标轴标注、对数坐标、线型细节。低于这个版本就用 print(handle, -dsvg, file.svg)。实际项目报告中还有一个讲究对数坐标下不要把频率轴从 0 开始画从 1 Hz 或 10 Hz 开始会显得谱峰区域饱满得多——这不是作弊只是展示习惯但审阅人看着舒服通过率就高。本文还有配套的精品资源点击获取