简介三分之一倍频程程序解读是一份面向软件网络技术及信号处理学习者的PDF资料重点讲解基于MATLAB实现三分之一倍频程声压级分析的完整流程。整包为单个PDF文档体积仅46KB轻便易用已有686人学习。内容提供两种典型实现方法方法一通过汉宁窗滤波与快速傅里叶变换计算未计权频谱并叠加A计权值得到计权声压级方法二从时域瞬时声压曲线切入利用局部峰值检测提取频带有效值再累加获得总声压级。代码中给出了中心频率数组、A计权值、汉宁窗恢复系数等关键参数的详细定义并配有瞬时声压时程图、频谱图和倍频程声压级图谱的绘制示例能帮助读者理解三分之一倍频程的算法原理并快速迁移到自己的音频信号分析任务中。整体内容精炼适合作为程序解读笔记和参考实现。1. 三分之一倍频程程序里的常数表就是最容易被误读的部分如果只看一份三分之一倍频程程序的运行结果很容易以为它只是把时域信号做了一次傅里叶变换再画个柱状图。真正接触过实现的人会知道这套流程里最微妙、也最容易被拿着代码直接跑的人略过的是一张提前写死在程序里的中心频率常数表——比如 31.5、63、125、250 这些看起来毫无规律的数字。几乎所有工程上常见的“三分之一倍频程程序”核心都不是怎么算 FFT而是怎么把这些频率与数字滤波器的带宽准确对应上以及如何在离散频谱里把能量正确划分到各个频带。这篇文章要做的事就是拿着一份常见的“三分之一倍频程程序”逐段拆开讲先厘清三分之一倍频程本身是什么再解释程序里每一块在做什么之后贴出可以直接复现的最小实现最后落到参数调整和结果验证上。适合的场景包括设备噪声测试报告的复查、声学仿真后处理脚本的编写以及需要把旧程序从 MATLAB 移植到 Python 的场合。读完你会具备一件事拿到一个别人的倍频程程序能迅速找出它的频带划分逻辑、窗口设置和能量累加方式并判断结果是否可信。2. 三分之一倍频程的频带划分原则中心频率、截止频率与常数带宽比2.1 为什么是“三分之一”而不是别的分数倍频程Octave的本质是频率轴的等比划分每个频带的上限频率是下限频率的 2 倍。一倍频程的带宽比较粗在噪声源识别时往往分不清两个相邻的离散频率成分。三分之一倍频程就是把每一个倍频程再切成三段每一段的频率上限是下限的 2 的 1/3 次方倍也就是约 1.2599 倍。这个划分方式保证了在任何频率位置频带的相对宽度都恒定因此它既不是等带宽也不是等百分比带宽而是处于两者之间的一种对数坐标下的均匀分割。这个特性直接决定了程序里的做法在频域做能量统计时不能像窄带谱分析那样按固定 Hz 数去积分能量而必须按每个频带的上下限频率去查频谱索引。程序里的常数表本质上就是这些频带边界的一组预计算值。大部分程序会把中心频率写成数组再根据中心频率乘上 2 的 ±1/6 次方来得到上下限而不是直接硬编码上下限。这样做的好处是当采样率变化时只需要检查边界频率是否低于奈奎斯特频率而不需要重新计算整张表。2.2 中心频率与截止频率的换算关系三分之一倍频程的标准中心频率由 ISO 266 和 IEC 61260 定义常用的基准是 1000 Hz。程序里常见的一组中心频率从 1 Hz 到 20 kHz 按 1/3 倍频程递增。对任意一个中心频率 ( f_c )其下限 ( f_1 ) 和上限 ( f_2 ) 的计算方式是[ f_1 f_c \times 2^{-1/6}, \quad f_2 f_c \times 2^{1/6} ]频带宽度为 ( f_2 - f_1 )相对带宽恒为约 23.1%。很多程序在初始化时会生成这样一张表然后用它去频谱里圈定每个频带的起止位置。需要特别指出的是程序里如果出现的是“标称中心频率”而不是实际边界那它通常已经根据 IEC 标准做了舍入处理。比如 31.5 Hz 这个标称值对应的实际中心频率是 31.62 Hz程序里列出的数值可能是 31.5 也可能是 31.62这取决于标准版本。标称中心频率 (Hz)实际中心频率 (Hz)下限频率 (Hz)上限频率 (Hz)31.531.6228.1835.486363.1056.2370.79125125.89112.20141.25250251.19223.87281.84500501.19446.68562.342.3 程序里频带表的两种组织方式第一种方式是只存中心频率数组运行时统一计算上下限。这种写法的好处是表短、易读、便于调整频带范围。第二种方式是直接把上下限频率也存成数组省去每次运行时的浮点计算。工业级声学软件偏向用第二种因为它们的频带划分可能要遵循多个标准GB/T 3241、IEC 61260、ANSI S1.11而不同标准在高频段的中心频率舍入规则略有差异。解读程序时首先就要确认它用的是哪一种。如果是第一种下一步要去看它计算边界用的是 ( 2^{\pm 1/6} ) 还是 ( 10^{\pm 1/20} )。两者在低频段几乎没有区别但在 10 kHz 以上会有微小偏差。严格来说三分之一倍频程的带宽比应该用 2 的 1/3 次方而 10 的 1/10 次方是早期某些声学软件采用的近似值两者相差约 0.3%高频段某些窄带峰值可能会因此落进相邻频带。3. 从窄带频谱到三分之一倍频程谱的三种实现路线3.1 直接滤波法模拟滤波器组的数字实现严格意义上的三分之一倍频程分析是让信号通过一组带通滤波器每个滤波器对应一个中心频率然后对每个滤波器的输出做 RMS 计算。数字实现时每个频带需要一个 IIR 滤波器通常是二阶或四阶 Butterworth 滤波器级联。这种方法的优点是时域分辨率好可以逐帧输出频带电平随时间的变化适合做噪声事件分析。代价是计算量大而且滤波器系数需要专门设计不同采样率下要重新计算。程序里如果出现butter、sosfilter之类的函数名基本走的就是这条路线。解读这类代码时最需要检查的是滤波器的阶数和归一化方式。常见做法是使用 4 阶 Butterworth 滤波器阻带衰减在倍频程带宽外约为 24 dB/oct。但也要注意如果程序里使用了filtfilt做零相位滤波那么实际等效阶数翻倍且引入了非因果处理适合离线分析不适合实时系统。3.2 FFT 能量积分法程序里最常见的高效路线绝大多数可读程序采用的方法是先对整个信号或分段信号做 FFT获得功率谱密度然后按每个频带的频率索引做能量累加。中心频率对应的频带能量等于该频带内所有谱线功率之和再换算成声压级的话要除以参考值并取 10 倍对数。这种方法有一个关键前提FFT 的频率分辨率必须足够高保证每个三分之一倍频程频带内至少包含几条谱线否则低频段会出现很大的统计误差。例如采样率 48000 Hz、FFT 点数 8192 时频率分辨率为约 5.86 Hz。对于 31.5 Hz 中心频率的频带带宽只有约 7.3 Hz只够容纳 1 到 2 条谱线这时的频带估计非常不稳。程序里常见的缓解策略是低频段用较小的 FFT 点数做分段平均高频段用较大的 FFT 点数——但这会让程序复杂度上升不少。还有一些程序直接设定最低分析频率为 50 Hz 或 100 Hz避开低频分辨率不足的问题。3.3 1/n 倍频程通用函数一个函数覆盖多种分辨率稍微完善一点的程序会写一个通用函数输入n倍频程分母和频带范围输出所有频带的索引区间。n1时就是一倍频程n3时就是三分之一倍频程。这种设计的好处是便于扩展例如后续需要做 1/6 倍频程或 1/12 倍频程分析时不需要重写逻辑。import numpy as np def octave_bands(n, freq_min, freq_max, fs, nfft): 生成 1/n 倍频程频带的频率索引区间。 参数: n: 倍频程分母3 表示三分之一倍频程 freq_min: 最低分析频率 (Hz) freq_max: 最高分析频率 (Hz) fs: 采样率 (Hz) nfft: FFT 点数 返回: bands: 列表每个元素为 (center_freq, idx_start, idx_end, f_low, f_high) # 基准中心频率 1000 Hz按 2^(1/n) 向两侧扩展 center_ref 1000.0 step 2 ** (1.0 / n) # 寻找覆盖 freq_min 到 freq_max 的所有中心频率指数 exp_start int(np.floor(np.log(freq_min / center_ref) / np.log(step))) exp_end int(np.ceil(np.log(freq_max / center_ref) / np.log(step))) # 频率分辨率 df fs / nfft bands [] for k in range(exp_start, exp_end 1): fc center_ref * (step ** k) f_low fc * (2 ** (-1.0 / (2 * n))) f_high fc * (2 ** (1.0 / (2 * n))) if f_high fs / 2: continue idx_start int(np.floor(f_low / df)) idx_end int(np.ceil(f_high / df)) if idx_end nfft // 2: idx_end nfft // 2 if idx_end idx_start: continue bands.append((fc, idx_start, idx_end, f_low, f_high)) return bands这段代码的逻辑核心是以 1000 Hz 为基准通过指数 k 向低频和高频两侧扩展中心频率再用中心频率计算频带边界。索引换算部分把频率除以频率分辨率df得到 FFT 谱线序号。注意这里用了floor和ceil把边界向外扩了一条谱线而不是四舍五入到最近的谱线。这个细节直接影响能量统计的完整性——如果把边界频率向内收缩会丢掉频带边缘的能量导致低频段结果偏低。3.4 三种路线的选型对照路线时域分辨率计算量实时性程序复杂度适合场景IIR 滤波器组高高可实时中噪声事件监测、时变信号分析FFT 能量积分低依赖帧长低可实时低稳态噪声分析、报告生成小波包分解中中依赖实现高非平稳信号、瞬态分析如果程序的目标是输出一张稳态噪声的频带声压级表FFT 能量积分法是性价比最高的选择。如果程序还需要同时输出每个频带的时间历程曲线那么滤波器组方法更合适。对于只做数据后处理的人来说FFT 法也是最容易验证和调试的——你可以随手造一个白噪声信号检查每个频带的能量是否接近理论值。4. 逐步解读一份典型程序的完整实现4.1 信号预处理去直流、加窗、分帧拿到一个程序后我一般会先跳过前面的导入和绘图部分直接看信号进入分析函数之前做了什么。最容易被忽略的是直流分量去除。实际从采集卡读入的数据常常带有微小的直流偏置如果不做去除FFT 在 0 Hz 附近会出现一个很高的谱峰这个峰会通过频谱泄漏影响到最低的几个频带尤其是 31.5 Hz 和 63 Hz 频带。常见做法是在加窗之前减去信号的均值。更讲究一点的做法是先用高通滤波器滤除 10 Hz 以下的成分但这会带来滤波器瞬态。对于稳态信号直接减均值就够用。加窗这一步也值得注意矩形窗的旁瓣太高会让强线谱的能量泄漏到相邻频带Hanning 窗是通用选择平顶窗适合校准场景但主瓣太宽在低频段会让相邻频带互相污染。def analyze_third_octave(signal, fs, nfft8192, overlap0.5): 对信号进行三分之一倍频程分析。 参数: signal: 输入时域信号 (numpy 数组) fs: 采样率 (Hz) nfft: FFT 点数必须是 2 的幂 overlap: 帧重叠比例0~0.75 之间 返回: freqs: 中心频率数组 levels: 每个频带的 RMS 声压级 (dB) # 去直流 signal signal - np.mean(signal) # 加窗 win np.hanning(nfft) win_norm np.sum(win ** 2) # 用于能量归一化 # 分帧参数 hop int(nfft * (1 - overlap)) n_frames (len(signal) - nfft) // hop 1 # 获取频带索引 bands octave_bands(3, 20, 10000, fs, nfft) # 初始化能量累加器 band_energy np.zeros(len(bands)) frame_count 0 for i in range(n_frames): start i * hop frame signal[start:start nfft] # 加窗并做 FFT spectrum np.fft.rfft(frame * win) power np.abs(spectrum) ** 2 / (fs * nfft) # 单边功率谱密度 # 对每个频带累加能量 for j, (_, idx_start, idx_end, _, _) in enumerate(bands): band_energy[j] np.sum(power[idx_start:idx_end]) frame_count 1 # 多帧平均 band_energy / frame_count # 转换为声压级 (参考 20e-6 Pa) # 假设信号单位已校准为 Pa p_ref 20e-6 levels 10 * np.log10(band_energy / (p_ref ** 2) 1e-12) freqs [band[0] for band in bands] return np.array(freqs), levels这里有几个参数值得展开说。nfft8192是典型的默认值但如果你分析的是 44100 Hz 采样率的音频信号频率分辨率为约 5.38 Hz如果后续要对低频尤其是 20 Hz 到 50 Hz 频带做精确估计这个分辨率明显不够需要增大 nfft 到 32768 或更高。overlap0.5是最常见的选择它让相邻帧共享一半数据降低帧与帧之间的不连续性。4.2 频带能量累加的细节为什么要用功率谱密度而不是幅度谱程序里把np.abs(spectrum) ** 2除以fs * nfft这行代码很容易被误改或误删。它的作用是得到功率谱密度PSD单位是 Pa²/Hz。如果不除以fs * nfft累加出来的能量单位就是 Pa²数值会随 nfft 变化而变化程序在参数调整时结果不稳定。这个归一化系数的物理含义是FFT 输出的每个谱线值是时域信号的幅度乘以窗函数傅里叶变换的结果幅度谱平方后需要对窗函数的能量做补偿再除以频率分辨率即可得到 PSD。需要留意的是这里用的是单边谱所以不需要额外乘以 2 进行折半但如果程序的其它部分用的是双边谱就要注意补上这一步。做程序解读时这一行是最值得标记注释的地方。4.3 输出声压级与参考声压程序最后把能量除以(20e-6) ** 2再取对数20 微帕是空气中的参考声压。如果程序处理的是水下声信号或振动加速度信号这个参考值必须改。水下声学常用 1 微帕作为参考声压振动分析用的是加速度的参考值 1e-6 m/s²。很多人照搬程序跑水下声数据结果整体声压级偏低 26 dB——这正是因为参考值用错了。处理程序移植时我通常会把参考值和信号校准系数写成函数的显式参数而不是硬编码在函数里。这个习惯能避免在换数据源时返工也让程序更容易被同事检查和复用。4.4 边界检查奈奎斯特频率与频带截断octave_bands函数里有一句if f_high fs / 2: continue这是保护性逻辑。如果分析频带的上限超过了采样率的一半程序就无法正确计算该频带的能量。更微妙的情况是频带上限刚好略高于奈奎斯特频率此时 FFT 输出会缺少一部分谱线能量累加会偏低。程序里把idx_end限制在nfft // 2内保证了不越界但代价是该频带的能量结果是不完整的。高采样率下这种情况不常见但如果把 48 kHz 采样率的数据拿来跑一个分析上限为 20 kHz 的通用程序在 20 kHz 频带上就会遇到部分截断。一个稳妥的检查方法是在程序末尾把idx_end对应的实际频率打印出来看是否覆盖了完整的频带。如果发现截断可以选择降低最高分析频率到 16 kHz或者重采样后再分析。5. 程序调试与参数调优三个最容易翻车的点5.1 窗函数与频谱泄漏的折中Hanning 窗是三分之一倍频程分析的默认选择但它不是所有场景的最优解。当信号里有较强的离散线谱比如变压器 100 Hz 的电磁噪声叠加在宽频背景噪声上时Hanning 窗的旁瓣会把线谱能量泄漏到相邻频带。这时候有两个方向增加 FFT 点数让线谱位于一条谱线附近或者改用旁瓣衰减更大的 Blackman-Harris 窗但代价是主瓣变宽低频段的频带间分辨能力降低。一个实用的判断标准是如果信号中离散线谱和连续噪声谱的强度差超过 20 dB就需要警惕泄漏。可以在程序里对比采用 Hanning 和 Blackman-Harris 的频带结果如果 63 Hz 频带的声压级在两个窗函数下差异超过 0.5 dB说明泄漏影响显著此时应重新设计窗长。窗函数主瓣宽度 (bin)旁瓣衰减适用场景矩形2-13 dB瞬态信号、校准Hanning4-31 dB通用噪声分析Blackman-Harris8-92 dB强线谱存在时平顶约 5-70 dB幅度精确测量5.2 RMS 计算与分段平均的统计稳定性三分之一倍频程的输出值本质上是频带内的 RMS 声压级。如果只对一段 1 秒的信号做一次 FFT 分析结果的随机不确定度会很大尤其是高频段因为每个频带内的有效独立样本数太少。程序里常见做法是把信号切成多帧、逐帧计算后平均。帧数越多结果越稳定但时间分辨率会变差。通用经验是每一帧的时长应至少覆盖频带最低频率的 10 个周期。比如要可靠估计 31.5 Hz 频带单帧时长至少需要 0.32 秒所以 8192 点的 FFT 在 44100 Hz 采样率下远远不够。程序调整时我一般会从低分析频率反推所需 FFT 点数f_low 20 # 最低分析频率 fs 44100 min_frame_duration 10 / f_low # 推荐至少 10 个周期 min_nfft int(2 ** np.ceil(np.log2(min_frame_duration * fs))) print(f最低频带 {f_low} Hz建议 FFT 点数至少为 {min_nfft})输出会是 32768 或更高。实际中很多程序为了追求速度用了 4096 或 8192导致低频结果方差很大。如果你拿到一份程序跑出来的结果在低频段抖动明显优先检查 FFT 点数而不是去怀疑滤波器设计。5.3 有效频段截断模拟滤波器无法覆盖整个频带在滤波器组实现方式中一个常见误区是认为程序可以输出所有标准中心频率的结果。实际上很多 FIR 或 IIR 滤波器实现只覆盖了采样率允许范围内的频带。比如采样率 48 kHz 时20 kHz 以上确实没有数据但程序里可能仍有 20 kHz 中心频率这一项只是输出结果是噪声底或 0。程序解读时要在输出结果里过滤掉那些f_high fs/2或能量接近数值下溢的频带。常见做法是按 0.85 倍奈奎斯特频率作为上限截断这样给滤波器过渡带留出余量。如果程序没有做这个截断输出的最高两个频带通常不可信。6. 用白噪声自检验证程序正确性验证三分之一倍频程程序最有效的手段是使用白噪声。白噪声的功率谱密度是平坦的理论上每个三分之一倍频程频带内的总功率应当正比于该频带的带宽。由于三分之一倍频程的相对带宽恒定约 23.1%因此每个频带的声压级应该完全相等。跑完程序后如果输出结果呈一条水平线说明频带划分和能量累加逻辑正确如果低频段出现明显下滑说明 FFT 点数不足如果高频段下垂则需要检查抗混叠滤波器或频带截断。具体操作时用 NumPy 生成一段 10 秒的白噪声rng np.random.default_rng(42) noise rng.standard_normal(10 * fs) * 0.1 # 幅度约 0.1 Pa然后分别用 8192、16384、32768 三组 FFT 点数跑同一段数据对比 20 Hz 到 63 Hz 频带的结果。三个结果之间的差异如果超过 1 dB说明最小 FFT 点数不满足低频估计需求。如果三次结果基本一致但整体水平线与理论值有固定偏移那么问题出在校准系数或参考声压上而不是算法本身。另一个验证点是频带能量守恒将所有三分之一倍频程频带的能量相加应当约等于信号的总功率。程序中可以用信号时域的均方值来验证total_power_time np.mean(noise ** 2) total_power_bands np.sum(band_energy) ratio total_power_bands / total_power_time print(f频带能量与总能量比值: {ratio:.3f})这个比值理论上应接近 1。如果明显小于 1说明有频带被漏掉或 FFT 归一化有误如果明显大于 1常见原因是多帧平均时没有正确除以帧数或加窗后没有做能量归一化。该验证放在程序输出的最后一步是一份工程报告中最值得附上的自检数据。本文还有配套的精品资源点击获取