MATLAB加速度信号功率谱密度计算与pwelch应用实践

简介：一个面向MATLAB信号处理与振动分析场景的算法示例，聚焦加速度信号模拟与功率谱密度（PSD）求解；压缩包共2个M文件，大小仅1KB，包含从正弦波生成、加窗预处理、FFT频谱变换到PSD归一化计算的完整代码流程。适合分析振动、冲击或噪声数据的工程师、科研人员，也适合准备深入学习频域方法的MATLAB用户；已有284人学习下载。学习该示例可掌握加速度信号仿真与PSD计算核心步骤：构造时间轴与正弦信号、选择窗函数抑制频谱泄漏、对FFT结果归一化并绘制对数功率谱；同时代码中引入GRMS指标计算，有助于理解加速度功率谱在工程振动评估、地震监测与机械故障诊断中的实际应用。两个M文件分别承担信号构造与统计算法演示，注释清晰、逻辑连贯，便于按行阅读和二次修改；通过调整频率、采样率等参数，还能观察不同信号对PSD估计的影响，深入理解窗函数、幅值归一化和频率分辨率等概念。

1. 从加速度信号到功率谱密度：为什么振动分析绕不开 PSD

拿到一段加速度计采集的时域信号，直接做 FFT 看幅值谱，是大多数 MATLAB 新手的第一反应。但放到实际工程里，无论是汽车 NVH、机床主轴振动监测，还是桥梁健康检测，你会发现正规报告里几乎只用功率谱密度（PSD，Power 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。用重力加速度 g（9.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 dB，Hanning 加窗配合 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 取 1024，fs 为 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/2+1); 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 必须通过具体数字理解，否则参数就是乱调的。

假设采样率 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_low=9 Hz，但频率分辨率只有 1 Hz，那 9 Hz 这个索引可能并不存在，取整到 f=9 或 f=10 会直接影响积分结果，误差可达 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 开始会显得谱峰区域饱满得多——这不是作弊，只是展示习惯，但审阅人看着舒服，通过率就高。

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