简介：面向音频信号处理初学者与Matlab使用者的实操型文档，围绕音频文件读取、快速傅立叶变换及分贝图绘制展开，帮助读者理解时域波形与频域频谱的对应关系，并掌握以对数尺度观察声音能量分布的方法。压缩包内共1个doc文件，约154KB，以图文结合的方式整理了音乐片段与安静环境敲击桌面两类音频的分析示例，便于对照不同信号在时域复杂度和频域成分上的差异。相关资源已有2195人学习下载，适合作为课程实验、信号分析入门或音频处理自学的参考材料。文档中的示例覆盖采样率、频率轴构建、FFT对称性处理以及20*log10归一化等关键环节，读者可据此理解分贝图在噪声识别、频率成分检测和响度比较中的实际用途，并迁移到音频编码、降噪与音乐制作等场景。
1. 分贝图不是一条 plot 曲线：先搞清楚你在画哪种 dB
Matlab 实现音频的分贝图，看上去一行plot(20*log10(abs(x)))就完事，真正让人返工的从来不是代码，而是纵轴那三个字母到底以什么为参考。同样一段录音，做异响排查的同事在 Matlab 里读到的峰值是 -12，旁边声级计屏幕上写着 78，两个数字都对，一个说的是 dBFS，一个说的是 dB SPL，中间隔着一整条校准链路。分贝本身是比值取对数，参考量一换，图就不是同一张图了。
这篇写给需要把音频变成可交付分贝图的人：设备噪声记录、产品异响定位、语音前端电平调试、声学课程实验都算。读完之后手上应该有三类图的能力——时间轴上的短时声级曲线、频率轴上的 dB 谱、1/3 倍频程柱状图，并且清楚每条曲线背后的参考量、计权方式和时间常数。不铺开讲声学计量，只保证图上每个数字都能追到代码里的某一行。
2. 从 PCM 采样到 dBFS 曲线：读文件、定参考、切帧
2.1 audioread 读进来的是什么，位深决定了什么
拿到一段 wav，第一件事不是画图，是确认手上这份数据是什么量纲。audioread对整数格式（16/24/32 bit PCM）会把码值归一化到 [-1, 1] 的 double，1.0 对应满量程；但对 32 bit 浮点格式的 wav，它不保证归一化，读进来可能就是原始浮点值。这个差异不分清楚，后面所有 dB 计算全部偏移一个常数。
我一般先跑一段探查代码，把采样率、位深、通道数、峰值、时长一次性打出来：
% 读文件前先看元信息，避免拿错量纲 info = audioinfo('sample.wav'); fprintf('采样率 %d Hz, 位深 %d bit, 通道 %d, 时长 %.2f s\n', ... info.SampleRate, info.BitsPerSample, info.NumChannels, info.Duration); [x, fs] = audioread('sample.wav'); fprintf('数据类别 %s, 峰值 %.4f, 最小值 %.4f\n', class(x), max(abs(x(:))), min(x(:))); if size(x, 2) > 1 x = x(:, 1); % 先只用第一通道，多通道合成放到后面单独讲 end位深直接决定理论动态范围：16 bit 是 96 dB，24 bit 是 144 dB。这不是学术数字，它意味着一段 16 bit 录音如果本底噪声录在 -80 dBFS，你在图上看不到任何低于它的细节，因为量化噪声已经把地板抬起来了。做低电平分析（比如设备待机噪声）时，24 bit 采集是硬门槛，用 16 bit 素材画出来的本底曲线只能算装饰。
还有一点容易被忽略：很多录音设备在模拟前端就做了自动增益或者限幅，文件本身没有削顶，但信号已经被压过。这时候 dBFS 曲线依然是"对"的，只是它反映的是处理后的电平，不是原始声压。要判断有没有被动过，看整段信号的峰值分布——正常录音峰值参差不齐，被压过的信号峰值会贴着某个固定值不动。
2.2 20log10 还是 10log10：峰值、RMS 与 dBFS 的边界
dBFS 的 0 点定义在满量程正弦波的峰值上，所以一个幅度为 1.0 的正弦波，峰值是 0 dBFS，而它的 RMS 是 0.7071，换算过来 -3.01 dBFS。这就是第一个高频误解：有人拿满幅正弦测出来 -3 dB，以为代码错了，其实是把峰值参考和 RMS 参考混用了。
换算关系只有两条：
| 量 | 公式 | 说明 |
|---|---|---|
| 峰值电平 | 20*log10(max(abs(x))) | 与 0 dBFS 峰值参考对应 |
| RMS 电平 | 20*log10(sqrt(mean(x.^2))) | 等价于10*log10(mean(x.^2)) |
两条公式在数学上是一回事（20log10(r) = 10log10(r^2)），区别只在你是先开方再取对数，还是先平方再取对数。写代码时我更倾向后一种，因为均方值可以跨帧累加、可以做时间计权、可以求 Leq，一路算下来不用反复开方，误差也小。
对整段信号算一个 RMS 电平，只在信号近似平稳时有意义。真实音频的动态范围动辄 40 dB 以上，整段 RMS 会把所有细节抹平。所以要画的是短时电平曲线，这就引出分帧。
2.3 分帧、加窗与时间轴对齐：帧长和帧移怎么定
短时电平的分帧参数不是随便拍的，它要和你后面要对比的东西对齐。如果最终图要和声级计的 Fast 档比对，帧长就应该取 125 ms；要和 Slow 档比对就取 1 s。帧移决定曲线的时间分辨率，10 ms 是个常用值，一秒钟 100 个点，图上既不会锯齿也不会糊成一团。
加窗这件事在纯电平计算里其实可以不写，但一旦你希望同一条代码路径既能出时域曲线又能出频谱，加上窗反而更统一。加窗后必须做功率补偿，否则 Hann 窗会把电平压低 4.26 dB，这个数字后面避坑章还会再提一次。
frameLen = round(0.125 * fs); % 125 ms，对齐 Fast 时间计权 hopLen = round(0.010 * fs); % 10 ms 帧移 win = hann(frameLen, 'periodic'); winPow = sqrt(mean(win.^2)); % 功率意义下的窗增益，Hann 约 0.6124 nFrames = floor((numel(x) - frameLen) / hopLen) + 1; lvlDbfs = zeros(nFrames, 1); for k = 1:nFrames idx = (k-1)*hopLen + (1:frameLen); seg = x(idx) .* win / winPow; % 补偿窗函数造成的电平损失 lvlDbfs(k) = 10*log10(mean(seg.^2) + eps); % eps 防止静音段出 -Inf end t = ((0:nFrames-1)*hopLen + frameLen/2) / fs; % 帧中心时刻，不是帧起点三个参数值得单独说。frameLen取 125 ms 是为了后续和 Fast 计权做交叉验证，如果只是想看包络，50 ms 也够用。hopLen取 10 ms 是折中，取 1 ms 会让曲线毛刺变多、文件变大，取 100 ms 会漏掉瞬态。win用'periodic'而不是默认的'symmetric'，因为周期窗在 FFT 场景下能量更准，同一份代码后面要复用到频谱章节。
t用的是帧中心而不是帧起点，这个细节在图上不明显，但当你把 Matlab 曲线和声级计导出的 CSV 按时间轴叠在一起时，半个帧的偏移（62.5 ms）会让两条曲线错开一个明显的台阶，看起来像是设备延迟，其实是画图时对齐方式不一致。
3. 把 dBFS 变成 dB SPL：校准、A 计权与时间计权
3.1 用 94 dB / 1 kHz 校准器反推系统偏移量
dBFS 是相对于数字满量程的，dB SPL 是相对于 20 μPa 的，两者之间差一个由传声器灵敏度、前置放大器增益、ADC 参考电压共同决定的常数。这个常数没法从文件里推出来，只能实测。常见做法是用声校准器在 94 dB SPL、1 kHz 下录一段，反算偏移。
% 校准器录音：94 dB SPL @ 1 kHz [cal, fsCal] = audioread('cal_94dB_1kHz.wav'); calRms = sqrt(mean(cal.^2)); % 归一化单位下的 RMS dbfsCal = 20*log10(calRms); offset = 94 - dbfsCal; % 之后 dB SPL = dBFS + offset fprintf('校准读数 %.2f dBFS，系统偏移 %.2f dB\n', dbfsCal, offset);offset就是整条链路的"黑匣子"钥匙，一次校准、全程复用，只要增益旋钮没动过。实操上有几个细节：校准器套上麦克风后要等两三秒读数稳定；录音时长取 5 到 10 秒，去掉头尾各 1 秒再算 RMS；校准文件本身的采样率要和被测文件一致，不一致就得先重采样，否则后续滤波器系数会错。
还有一种情况是没有校准器。这时候可以退而求其次，用一个已知声压级的稳态声源做一次单点标定，比如某台设备的标称噪声值。这种做法的误差可能有 2 到 3 dB，但至少量纲是对的，图能看。最怕的是既没有校准器，又直接在图上写 dB SPL——那个数字除了自己没人敢信。
3.2 A 计权：频域一次算清，别急着设计滤波器
人耳对低频不敏感、对 2 到 4 kHz 略敏感，A 计权就是把这套等响特性固化成一条频率响应曲线。IEC 61672 给出了解析表达式：
function wdB = aWeightDb(f) % f 为频率向量，单位 Hz，返回对应 A 计权修正量（dB） f2 = f.^2; ra = (12194^2 .* f2.^2) ./ ... ((f2 + 20.6^2) .* sqrt((f2 + 107.7^2) .* (f2 + 737.9^2)) .* ... (f2 + 12194^2)); wdB = 20*log10(ra) + 2.00; % 归一化，使 1 kHz 处为 0 dB wdB(f == 0) = -Inf; % 直流分量不计权 end有了这个函数，A 计权可以用两种方式落地。一种是频域：对每帧做 FFT，把幅度谱乘以10.^(aWeightDb(f)/20)，再 IFFT 回时域。另一种是时域：用双线性变换把模拟传函离散成 IIR 滤波器，再filter一遍。
我一般选频域，原因是它只依赖频率轴，不受采样率变化影响，也不会因为滤波器阶数不够在低频端引入相位失真。代价是需要分帧且要求帧内有足够频率分辨率——帧长 125 ms 在 1 kHz 处只有 8 Hz 分辨率，对低频段（20 到 100 Hz）不够，所以对低频做计权时帧长要拉到 1 s 以上。
只用整段 Leq 时更省事，一次 FFT 搞定：
N = 2^nextpow2(numel(x)); X = fft(x .* hann(numel(x)), N); f = (0:N-1)' * fs / N; mag = abs(X) / sum(hann(numel(x))); % 幅度归一化 magDb = 20*log10(mag + eps); magDbA = magDb + aWeightDb(min(f, fs/2)); % 只对 0~fs/2 有意义 leqA = 10*log10(sum(10.^(magDbA/10))) + offset;这里sum(10.^(magDbA/10))是把各频点的功率相加，等价于 Parseval 定理下的总能量，再取一次对数得到 A 计权总级。注意offset加在最后，因为它是一个与频率无关的常数，先加后加结果一样，但放最后代码更好读。
3.3 Fast / Slow / Impulse：时间计权与 Leq 的滚动计算
时间计权是对"能量"做一阶低通，不是对"分贝值"做低通。这个顺序错了，曲线形状会完全不同。标准时间常数是 Fast 125 ms、Slow 1000 ms、Impulse 35 ms（上升）加 1500 ms（下降），实现上前两个用同一套指数平滑即可：
p2 = x.^2; % 瞬时均方 taus = struct('Fast', 0.125, 'Slow', 1.0); % 单位：秒 for name = fieldnames(taus)' T = taus.(name{1}); alpha = exp(-1/(fs*T)); % 一阶 IIR 系数 y = filter(1-alpha, [1 -alpha], p2); % y[n] = (1-a)x[n] + a*y[n-1] lvl = 10*log10(y + eps) + offset; % 转 dB SPL plot((0:numel(lvl)-1)/fs, lvl); hold on; end Leq = 10*log10(mean(p2)) + offset; % 整段等效声级alpha的物理含义是"上一个采样点还剩多少权重"，采样率越高，alpha越接近 1。用filter而不是手写循环，是因为filter内部对初值处理更干净，而且快一个数量级。eps那一下是保命的，静音段p2全零，log10(0)会返回-Inf，画到图上坐标轴直接崩掉。
Leq（等效连续声级）是不带时间计权的，它就是对整段均方取平均再转 dB。这一点经常被误解：有人拿 Fast 曲线做平均得到 Leq，结果和标准算法差零点几 dB，因为mean(10.^(lvl/10))是均方的平均，理论上等价，但曲线本身已经过指数平滑，等于做了两次滤波。正确做法是回到原始p2求平均。
如果要的是滚动 Leq（比如每秒一个值），用累积和实现，避免每帧重新求和：
winSec = 1; % 1 秒滚动窗 nWin = round(winSec * fs); cs = [0; cumsum(p2)]; % 累积和，O(N) 复杂度 idxEnd = nWin:round(0.1*fs):numel(p2); leqRoll = 10*log10((cs(idxEnd+1) - cs(idxEnd-nWin+1)) / nWin) + offset;4. 频谱类分贝图：spectrogram 的 dB 色标与倍频程谱
4.1 spectrogram 从线性幅度到 dB：参考值定错整张图就废
spectrogram默认输出的是线性幅度或功率，直接imagesc画出来是一片黑，因为绝大多数能量集中在很小的动态范围内。转 dB 之前必须先确定参考量，否则色标上的数字没有物理含义。
参考量有三种常见选择，对应三种图：
| 参考量 | 计算方式 | 图上的含义 |
|---|---|---|
| 满量程 | 20*log10(abs(S)) | dBFS，用于数字域相对比较 |
| 声压 | 上式 +offset | dB SPL，需要校准 |
| 功率谱密度 | 10*log10(Pxx) | dB/Hz，跨窗长可比 |
第三种最容易被忽略，但它才是能横向对比的画法。因为 FFT 的 bin 宽度随帧长变化，同样一段噪声，帧长加倍后每个 bin 分到的能量减半，峰值看起来掉了 3 dB——这不是信号变了，是分辨率变了。功率谱密度除以等效噪声带宽后，这个依赖就消失了。
nfft = 2048; win = hann(nfft, 'periodic'); overlap = round(nfft * 0.75); [S, F, T] = spectrogram(x, win, overlap, nfft, fs, 'psd'); Sdb = 10*log10(S + eps) + offset; % dB SPL^2/Hz figure; imagesc(T, F/1000, Sdb); axis xy; % 低频在下，符合直觉 caxis([max(Sdb(:))-80, max(Sdb(:))]); % 动态范围截到 80 dB colormap(jet); colorbar; xlabel('时间 (s)'); ylabel('频率 (kHz)');'psd'选项让 Matlab 内部完成窗增益和 bin 宽度的归一化，省掉手算。caxis那行决定了图好不好看：不设上限，最大值会吃掉整个色表；设得太窄，本底噪声全被截平。常用范围是峰值往下 60 到 90 dB，具体看你要看的是主分量还是背景纹理。
overlap取 75% 是权衡，重叠越高图越平滑但计算量越大。做瞬态分析（敲击、脉冲）时重叠要拉到 90%，否则短事件会落在两帧之间被漏掉。
4.2 pwelch 算功率谱密度并给出本底噪声
要看一段信号"整体有多响、能量分布在哪"，用pwelch做 Welch 平均比单次 FFT 稳得多。它把信号切成多段、每段加窗做 FFT、再对功率谱求平均，随机噪声的方差能降一个数量级，代价是频率分辨率受段长限制。
nfft = 8192; [Pxx, F] = pwelch(x, hann(nfft,'periodic'), nfft/2, nfft, fs, 'psd'); PxxDb = 10*log10(Pxx + eps) + offset; % 本底噪声估计：取谱的低分位数，避开离散谱线 noiseFloor = prctile(PxxDb, 10); fprintf('本底噪声约 %.1f dB SPL/Hz\n', noiseFloor); semilogx(F, PxxDb); xlim([20 fs/2]); yline(noiseFloor, '--', '本底');用 10% 分位数而不是最小值，是因为最小值受单点异常影响太大，而分位数对离散谱线不敏感。这条本底线在图上很有用：它告诉你这套采集系统自身的底噪在哪，任何低于它的分析结果都是不可信的。
要注意pwelch的窗长和你要看的频率下限匹配。想看 20 Hz，窗长至少要 1 s（分辨率 1 Hz），8192 点在 48 kHz 下只有 0.17 s，20 Hz 附近的分辨率约 6 Hz，勉强够看趋势，看细节就不行了。这个取舍没有标准答案，我的习惯是低频分析用 1 s 窗，宽带噪声用 0.17 s 窗。
4.3 1/3 倍频程带能量合成与柱状分贝图
倍频程图是声学报告里最常见的形式，因为它的频率轴按对数分布，和人耳感知更接近。1/3 倍频程就是把每个八度切成三段，中心频率按 10^(0.1) 递增，带宽比是 10^(0.1) ≈ 1.26。
% 中心频率：25 Hz 起，按 1/3 倍频程递增 fc = 1000 * 10.^((-13:13)/10); fc = fc(fc >= 20 & fc < fs/2); % 裁掉超范围的部分，含 20 Hz 下限 fLow = fc / 10^0.05; % 下限截止 fHigh = fc * 10^0.05; % 上限截止 [Pxx, F] = pwelch(x, hann(8192,'periodic'), 4096, 8192, fs, 'psd'); df = F(2) - F(1); Lb = zeros(size(fc)); for k = 1:numel(fc) m = (F >= fLow(k)) & (F < fHigh(k)); Lb(k) = 10*log10(sum(Pxx(m)) * df + eps); % 带内均方求和 end Lb = Lb + offset; bar(log10(fc), Lb, 0.6); % 对数横轴用 log10 后手动画 set(gca, 'XTick', log10([31.5 63 125 250 500 1000 2000 4000 8000]), ... 'XTickLabel', {'31.5','63','125','250','500','1k','2k','4k','8k'}); xlabel('中心频率 (Hz)'); ylabel('频带声压级 (dB SPL)');三个参数需要留意。fc的起点决定低端覆盖，20 Hz 以下人耳几乎听不到，但设备振动分析有时需要，按需下调。df是 bin 宽度，乘以它才把功率谱密度还原成带内功率，漏乘这一步所有柱子的高度会低 30 到 40 dB，是个非常隐蔽的错误。bar的横轴用log10(fc)再手动设刻度，是因为 Matlab 的bar对semilogx支持不好，这个做法比用semilogx加stairs更省事。
如果要做 A 计权的倍频程图，不必对信号先滤波，直接把aWeightDb(fc)加到每根柱子上就行，因为每个倍频带内的计权量近似为常数。这个近似在低频端误差稍大，但工程上够用。
5. 避坑排查：分贝图上蹿下跳的 5 个真实原因
5.1 曲线整体低约 96 dB：dBFS 直接当成了 dB SPL
现象：画出来的曲线全部集中在 -90 到 0 之间，而设备标称噪声是 60 多分贝。
原因：没做校准，纵轴还是 dBFS。16 bit 系统的理论动态范围正好 96 dB，所以图上数值和实际声压级差着一个约 90 到 100 的常数，看起来"很有规律"，容易被误认为是设备本底。
解决：跑一遍校准器流程拿到offset，所有输出统一加这个常数，并且在纵轴标签里明确写dB SPL。没有校准条件时至少要标注dBFS，绝不能两个词混着写。
5.2 A 计权后高频翘起、低频下不去：滤波器系数和采样率错配
现象：加了 A 计权，1 kHz 附近确实在 0，但 8 kHz 处反而比未计权时高，低频衰减也不够。
原因：走的是 IIR 设计路线，系数是从某个固定采样率的模拟传函离散来的，换了采样率没重算，或者双线性变换没做频率预畸变，导致高频段传函畸变。
解决：换成频域计权，只依赖f向量，采样率怎么变都对。如果一定要用 IIR，设计完必须用freqz把实际响应画出来，和标准表逐点比对，差 0.5 dB 以上就重来。
5.3 静音段出现 -Inf 或 NaN，曲线断成几截
现象：波形中间有停顿，分贝曲线在那一段画不出来，坐标轴被拉爆。
原因：log10(0)返回-Inf，后续做平均时Inf - Inf又变成NaN，绘图函数遇到这些值直接跳过。
解决：所有取对数的地方统一加eps，或者在加之前用max(p2, realmin)兜底。更稳的做法是在画图前用ylim显式限制显示范围，比如[0 120]，把异常值挡在视图外面。同时留意不要把eps加在 dB 值上——那几乎没作用，必须加在功率或幅度上。
5.4 加窗后整条曲线平移 4 dB 多：补偿用了相干增益
现象：同一段信号，加了 Hann 窗之后电平比不加窗低了将近 4.3 dB，换个窗函数平移量还不一样。
原因：窗函数对幅度的衰减有两个不同的度量——相干增益sum(w)/N（Hann 是 0.5，即 -6 dB）和功率增益sqrt(mean(w.^2))（Hann 是 0.6124，即 -4.26 dB）。算 RMS 电平时必须用功率增益，用相干增益会多扣 1.76 dB；算单频幅度谱时反过来要用相干增益。
解决：在代码里把两个量显式命名，winCoh和winPow各存一份，用在哪一步写清楚。这个坑最烦的地方在于它不会报错，只是让数字系统性偏移，如果正好和别处的误差部分抵消，可能很久都发现不了。
5.5 两个工具算出的 Leq 差 3 dB：双声道合并方式不一致
现象：同一段立体声录音，Matlab 算出 78.2 dB，另一个工具算出 81.1 dB。
原因：一边用mean(x, 2)做了幅度平均，另一边用sqrt(mean(x.^2, 2))做了能量合成。当两个通道信号不相关时，幅度平均会让能量降一半，也就是 3 dB；当两个通道完全相同时，两种做法结果一致。真实录音通常介于两者之间，所以差异在 0 到 3 dB 之间浮动。
解决：噪声测量统一用能量合成：
if size(x, 2) > 1 x = sqrt(mean(x.^2, 2)); % 能量域合成，保留总声压级 end如果是做波束形成之类的相位相关处理，那另当别论，但只要目的是"测这段录音有多响"，能量合成才是对的。
6. 让图能交付：坐标规范、自检信号与一次链路的回归
图能画出来和能交付是两件事。前者只要plot不报错，后者要求任何人拿到这张图，不看代码也能知道纵轴参考是什么、计权是什么、时间常数多少。我的做法是把这三项直接写进轴标签：声压级 (dB SPL, A 计权, Fast)，一个词都不省。色标图同理，在colorbar的标签里写明dB SPL/Hz。
导出用exportgraphics，300 dpi 起步，比saveas稳定：
exportgraphics(gcf, 'level_fast.png', 'Resolution', 300, 'BackgroundColor', 'white');真正省时间的不是画图，是自检。每次改完链路，我会跑一遍合成信号回归：生成一段 1 kHz、峰值 -20 dBFS 的正弦，走完全流程，看输出的 dB SPL 是否等于-20 + offset（A 计权在 1 kHz 处为 0，所以不用修正）。再生成一段白噪声，跑 A 计权后的 1/3 倍频程，和标准表比对：
| 中心频率 (Hz) | 63 | 125 | 250 | 500 | 1k | 2k | 4k | 8k |
|---|---|---|---|---|---|---|---|---|
| A 计权 (dB) | -26.2 | -16.1 | -8.6 | -3.2 | 0.0 | +1.2 | +1.0 | -1.1 |
白噪声在理想情况下各频带原始能量相同，所以计权后的柱子形状应该和这张表吻合，误差超过 0.5 dB 就说明滤波器或频带合成有偏差。这个自检比我见过的任何调试手段都有效，因为它把"看起来差不多"变成了可量化的比对。
血泪经验是：不要把校准和计权写成两段独立脚本。我早期就是分开的，改校准常数时忘了同步到计权脚本，结果一批图全部偏移 2.3 dB，等到发现时报告已经发出去了。现在我把offset、aWeightDb、分帧函数放进同一个工具箱目录，任何一处改动都会影响全部输出，回归测试跑一遍就知道有没有连带损伤。
另一个习惯是保留原始数据路径和参数结构体。图上好看不好看是次要的，半年后有人问"这个 78.6 是怎么来的"，能一键复现才算数。
希望帮到你。
本文还有配套的精品资源，点击获取
