简介面向音频信号处理初学者与Matlab使用者的实操型文档围绕音频文件读取、快速傅立叶变换及分贝图绘制展开帮助读者理解时域波形与频域频谱的对应关系并掌握以对数尺度观察声音能量分布的方法。压缩包内共1个doc文件约154KB以图文结合的方式整理了音乐片段与安静环境敲击桌面两类音频的分析示例便于对照不同信号在时域复杂度和频域成分上的差异。相关资源已有2195人学习下载适合作为课程实验、信号分析入门或音频处理自学的参考材料。文档中的示例覆盖采样率、频率轴构建、FFT对称性处理以及20*log10归一化等关键环节读者可据此理解分贝图在噪声识别、频率成分检测和响度比较中的实际用途并迁移到音频编码、降噪与音乐制作等场景。1. 分贝图不是一条 plot 曲线先搞清楚你在画哪种 dBMatlab 实现音频的分贝图看上去一行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] 的 double1.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 dB24 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(idxEnd1) - cs(idxEnd-nWin1)) / nWin) offset;4. 频谱类分贝图spectrogram 的 dB 色标与倍频程谱4.1 spectrogram 从线性幅度到 dB参考值定错整张图就废spectrogram默认输出的是线性幅度或功率直接imagesc画出来是一片黑因为绝大多数能量集中在很小的动态范围内。转 dB 之前必须先确定参考量否则色标上的数字没有物理含义。参考量有三种常见选择对应三种图参考量计算方式图上的含义满量程20*log10(abs(S))dBFS用于数字域相对比较声压上式 offsetdB 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 Hz8192 点在 48 kHz 下只有 0.17 s20 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 dBdBFS 直接当成了 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)/NHann 是 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。导出用exportgraphics300 dpi 起步比saveas稳定exportgraphics(gcf, level_fast.png, Resolution, 300, BackgroundColor, white);真正省时间的不是画图是自检。每次改完链路我会跑一遍合成信号回归生成一段 1 kHz、峰值 -20 dBFS 的正弦走完全流程看输出的 dB SPL 是否等于-20 offsetA 计权在 1 kHz 处为 0所以不用修正。再生成一段白噪声跑 A 计权后的 1/3 倍频程和标准表比对中心频率 (Hz)631252505001k2k4k8kA 计权 (dB)-26.2-16.1-8.6-3.20.01.21.0-1.1白噪声在理想情况下各频带原始能量相同所以计权后的柱子形状应该和这张表吻合误差超过 0.5 dB 就说明滤波器或频带合成有偏差。这个自检比我见过的任何调试手段都有效因为它把看起来差不多变成了可量化的比对。血泪经验是不要把校准和计权写成两段独立脚本。我早期就是分开的改校准常数时忘了同步到计权脚本结果一批图全部偏移 2.3 dB等到发现时报告已经发出去了。现在我把offset、aWeightDb、分帧函数放进同一个工具箱目录任何一处改动都会影响全部输出回归测试跑一遍就知道有没有连带损伤。另一个习惯是保留原始数据路径和参数结构体。图上好看不好看是次要的半年后有人问这个 78.6 是怎么来的能一键复现才算数。希望帮到你。本文还有配套的精品资源点击获取