ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

MATLAB实现A计权1/3倍频程声压级分析完整指南

MATLAB实现A计权1/3倍频程声压级分析完整指南 简介面向环境噪声、工业设备噪声及声学测试数据的合规性评估这套MATLAB脚本实现A计权处理与1/3倍频程声压级计算。脚本内置ISO 266和IEC 61672推荐的滤波器系数支持自定义采样率与频带范围即使没有额外工具箱也可在MATLAB R2015b及以上版本直接运行适合声学测试人员与相关专业学生快速掌握频带级噪声分析方法。压缩包内共8个文件除MATLAB主程序外还附带了Python辅助脚本、txt说明文档以及png曲线图等便于理解A计权曲线与计权网络特性整体大小约250KB结构精简、查用方便。已有102人学习参考。输出结果以结构化数据集呈现逐一给出各1/3倍频程中心频率对应的A加权声压级可直接用于后续绘图、导出或与限值比对在环境噪声监测、工业设备噪声评价及声学测试数据复现中具有较高实用价值。 我最近在整理设备噪声测试数据时正好需要把一批时域录音转换成 A 计权 1/3 倍频程声压级频谱。翻了一圈现成工具要么是商业软件收费不低要么是自带算法和声级计对不上干脆自己写了一个 MATLAB 脚本。这个脚本的核心功能就一句话输入一段声压时域信号输出 1/3 倍频程各频带的 A 计权声压级附带频谱图。如果你也在做环境噪声评估、产品噪声排查或者想跟声级计实测结果互相对标这篇文章里从算法原理到代码实现、从校准到避坑应该都能省你不少时间。这里先给结论用 MATLAB 做 A 计权 1/3 倍频程声压级分析完全可行而且如果按标准公式正确实现误差可以控制在零点几分贝以内。但前提是你得明白两条技术路线的区别以及在低频段和滤波器设计上会踩哪些坑。下面按我实际开发的思路逐步拆开讲。1. 为什么噪声分析要用 A 计权和 1/3 倍频程1.1 人耳听觉特性与 A 计权的由来噪声评价里最常出现的指标就是 dB(A)也就是 A 计权声压级。人耳不是对所有频率都同样敏感同样是 60 dB 的声压级1 kHz 声音听起来往往比 100 Hz 的响得多。A 计权曲线就是模拟人耳在中等响度下的频率响应特性对低频做了大幅度衰减对高频略有提升在 1 kHz 位置修正量为 0 dB。工程上无论环境噪声标准、家电产品噪声限值还是职业健康评估绝大多数都直接采用 A 计权声压级作为最终评价量。所以我们在做声学测试时不能只看宽带总声压级更要看频率分布而 A 计权恰好把“人耳听到的大小”这个维度放进了每一段频谱里。1.2 1/3 倍频程的分辨率为什么够用用 FFT 做频谱分析可以得到几百上千条谱线但人耳对频率的分辨能力是有限的更接近对数频率尺度下的恒定比例分辨率。1/3 倍频程把整个可听频段划分为一系列频带每个频带的带宽是中心频率的 23% 左右相邻中心频率之比是 2 的 1/3 次方。20 Hz 到 20 kHz 范围内大约分布 30 个频带这已经能很好地识别噪声的主要特征比如低频嗡嗡声、中频风扇噪声、高频气动噪声。相比窄带 FFT1/3 倍频程谱更贴近听觉感知也方便和标准限值表中的频带声压级直接对比。如果你最后要出噪声检测报告用的几乎都是 1/3 倍频程谱而不是几百根的 FFT 谱线。1.3 脚本要解决的实际问题我需要的工具具体包括读取一段校准后的声压信号单位 Pa、按标准 1/3 倍频程中心频率划分频带、在每个频带内计算 A 计权声压级、以 dB(A) 为单位输出频谱表和图形。这个脚本不是要替代专业声学软件而是用来在开发阶段快速评估噪声源、对比不同方案的降噪效果顺便验证第三方测试报告的数据。2. 脚本整体架构两条技术路线怎么选2.1 路线一时域滤波链路标准但重经典声级计内部处理流程是这样的原始信号先经过带通滤波器组分成 1/3 倍频程各频带信号每一路分别经过 A 计权滤波器再求 RMS 值得到声压级。这种方式的优点是严格符合 IEC 61672 标准能处理非稳态、瞬态信号还能做实时连续监测。缺点也很明显需要为每个频带设计带通滤波器低频段滤波器阶数高了容易数值不稳阶数低了衰减不陡峭频率边界容易出问题。加上还要做 A 计权时域滤波代码量和调试成本都不小。如果是长时连续噪声监测这条路是必须的但只是为了分析一段几分钟的录音有点杀鸡用牛刀。2.2 路线二FFT 频带能量积分轻快但有前提另一条思路是对整段信号做 FFT得到功率谱密度然后按 1/3 倍频程频带边界把功率谱积分。这种方法的原理是帕塞瓦尔定理时域 RMS 平方等于频域能量积分。由于 FFT 本身已经包含了全部频率信息就可以在频域直接做 A 计权加权不需要设计任何滤波器。它的局限在于要求信号在分析时长内近似平稳同时低频段频率分辨率受制于数据长度。比如采样率 48 kHz、分析时长 1 秒FFT 频率分辨率为 1 Hz在 20 Hz 频带内只有 7 个频点积分误差会比较大。但对于常规产品噪声测试采样几秒甚至几十秒完全可行实测结果跟声级计对比差值通常在 0.5 dB 以内。2.3 我的最终选型我的需求是快速出频谱图、方便批量跑不同工况数据不需要实时所以我选择 FFT 频带积分法做主体脚本同时留一个滤波法接口用于小信号比对。再强调一句如果要出正式报告或用来做型式试验建议还是老老实实用符合 IEC 61260 的滤波链路或者直接用专业声学前端。FFT 积分法更适合研发自测和前期排查。3. A 计权滤波器工程实现从模拟公式到 MATLAB 代码3.1 标准 A 计权频率响应公式A 计权在频域有明确的解析式。IEC 61672 给出的模拟频率响应是R_A(f) (12194^2 * f^4) / ((f^2 20.6^2) * sqrt((f^2 107.7^2) * (f^2 737.9^2)) * (f^2 12194^2))然后换算成以 dB 为单位的 A 计权修正量A_db(f) 20 * log10(R_A(f)) 2.00这个 2.00 dB 的常数是为了让 1 kHz 处的修正量正好等于 0 dB。注意有些老资料里没加这个 2 dB直接套用会发现 1 kHz 附近总差 2 dB这是一个很容易翻车的点。MATLAB 里写一个简单的函数function A A_weight(f) f1 20.6; f2 107.7; f3 737.9; f4 12194; num f4^2 * f.^4; den (f.^2 f1^2) .* sqrt((f.^2 f2^2) .* (f.^2 f3^2)) .* (f.^2 f4^2); A 20 * log10(num ./ den) 2.00; end3.2 频域加权实现方式FFT 积分法中最直接的 A 计权应用是把频带内每一条谱线的功率密度乘以对应频率的线性加权系数再做频带积分。线性加权系数是 10^(A_db/10)因为 A_db 本身是 20log10 表示的电压/声压比功率比就是平方关系。对应代码片段W_linear 10.^(A_weight(f_band) / 10); energy_A trapz(f_band, pxx_band .* W_linear);其中f_band是当前频带内的频率向量pxx_band是功率谱密度。这样做的精度比单纯取中心频率加权高很多尤其是 200 Hz 以下A 计权曲线斜率变化很快用中心频率一点代替整个频带会带来明显误差。3.3 时域 IIR 滤波需要实时运行时用如果你的应用是实时测量或信号很长不能整段 FFT那需要设计 A 计权时域滤波器。工程上常用二阶级联滤波器逼近模拟曲线。MATLAB 里可以用designfilt尝试d designfilt(weighting, FilterOrder, 4, ... FrequencyWeights, A, SampleRate, fs); x_a filtfilt(d, x);不过designfilt的加权滤波器功能依赖 DSP System Toolbox如果没有这个工具箱可以用频域法做替代或者把模拟 A 计权传递函数通过双线性变换离散化自己写双二阶滤波。我个人建议非实时场景优先用频域加权简单且不容易出错。4. 1/3 倍频程带通滤波与声压级计算4.1 中心频率与频带边界生成1/3 倍频程中心频率以 1000 Hz 为基准按 2^(1/3) 等比排列。标准列出的常用频率从 1 Hz 到 20 kHz 共有几十个但工程测试一般只关心 20 Hz 到 20 kHz。可以用一行代码生成fc 1000 * 2.^((-30:30)/3); % 生成以 1000 Hz 为基准的系列 fc fc(fc 20 fc 20000); % 截取需要的范围 fc [20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 ... 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000 12500 16000 20000];建议直接用上面这个列表因为标准标称频率是规定好的自己算容易把 31.5 写成 31.25 这类误差带进去。每个频带的下限和上限分别为fl fc / 2^(1/6)fh fc * 2^(1/6)4.2 如果采用滤波法带通滤波器怎么设计滤波法需要为每个频带设计带通滤波器。这里我推荐用 4 阶 Butterworth双边零相位滤波。设计一个频带的例子fc_i 1000; fl fc_i / 2^(1/6); fh fc_i * 2^(1/6); [b, a] butter(4, [fl, fh] / (fs / 2), bandpass); y filtfilt(b, a, x); Lp_i 20 * log10(rms(y) / 2e-5);注意butter的归一化频率以 Nyquist 频率为 1所以截止频率要除以 fs/2。所有频带共用一段原始信号x逐频带滤波后计算 RMS 即可。但滤波法计算量较大低频段滤波器阶数不足时通带边界精度差需要小心验证。如果只是想让结果接近标准建议把 Butterworth 阶数提高到 6 到 8并用freqz观察通带纹波。4.3 声压级计算、校准与总声压级合成无论哪种方法最后都回到声压级公式L_p 20 * log10(p_rms / p_ref)其中 p_ref 20 μPa 2e-5 Pa如果用periodogram计算功率谱密度频带内 A 计权均方声压为energy_A那么该频带声压级就是L_pA 10 * log10(energy_A / p_ref^2)这里用 10log10 是因为 energy_A 本身是平方量Pa^2。校准是另一个容易忽略的点。如果你的信号已经是校准好的 Pa 值直接算没问题。但大多数采集系统导出的只是数字码值需要先乘以校准系数。标准做法是用声校准器在 1 kHz 位置输出 94 dB 声压级对应有效声压为 p_ref * 10^(94/20) ≈ 1 Pa。采集到这段信号后计算数字 RMS校准系数就是目标 RMS 除以实测数字 RMS。有了这个系数后面所有时域数据乘上去就能转成实际声压。总 A 计权声压级可以用频带能量叠加L_total 10 * log10(sum(10.^(L_pA / 10)))但要注意频带叠加得到的总级和直接对全频带 A 计权信号求 RMS 的结果应当非常接近。如果二者差超过 0.5 dB多半是频带覆盖不完整或者某个频带积分有问题。我最终写的核心计算循环大致是这样的p_ref 2e-5; N length(x); [pxx, f] periodogram(x, hanning(N), N, fs); % 单边功率谱密度 LpA zeros(size(fc)); for i 1:length(fc) fl fc(i) / 2^(1/6); fh fc(i) * 2^(1/6); idx f fl f fh; f_band f(idx); pxx_band pxx(idx); W 10.^(A_weight(f_band) / 10); energy_A trapz(f_band, pxx_band .* W); LpA(i) 10 * log10(energy_A / p_ref^2); endtrapz用梯形积分是为了适应频带内频率点不均匀的问题虽然periodogram返回的频率是均匀分布的但用trapz更通用例如从pwelch输出非均匀频率时也能用。绘图可以直接用barfigure; bar(1:length(fc), LpA); set(gca, XTick, 1:4:length(fc), XTickLabel, fc(1:4:end)); xlabel(中心频率 / Hz); ylabel(A计权声压级 / dB(A)); grid on;如果频带太多横轴用对数坐标更清晰可以用semilogx画折线或者直接把fc映射到等距位置方便观察刻度。5. 实测校验与避坑我踩过的几个坑5.1 低频段结果抖得厉害第一次用 1 秒数据跑 31.5 Hz 频带发现声压级比声级计低 3 个多 dB。原因是 1 秒数据在 31.5 Hz 频带内只有大约 31 个周期而且 Hann 窗对低频分量的能量泄露影响很大。FFT 频率分辨率只有 1 Hz31.5 Hz 频带带宽才 9 Hz 左右频点太少积分误差自然大。解决办法有三个一是把分析时长加到至少 5 到 10 秒频率分辨率到 0.2 Hz 以下二是用periodogram时做重叠平均降低方差三是对低频段单独使用滤波法频带 RMS 更稳定。我后来统一把测试录音时长加长到 10 秒低频结果就稳了和声级计对上了。5.2 和声级计对不上排查出 2 dB 偏差有个朋友跑我脚本发现所有频段都比声级计高 2 dB问我怎么回事。我一看他用的 A 计权公式没加 2.00导致 1 kHz 处不是 0 dB而是 -2 dB等于给全部结果加了 2 dB。这类问题在网上一堆代码里都有注意核对 1 kHz 处的修正量是否等于 0。另一个常见偏差是校准系数没对声级计显示 94 dB采集系统数字 RMS 可能不是 1必须用校准系数统一到 Pa 单位。否则所有频带整体平移差值固定。5.3 带通滤波器边界频点归属问题FFT 积分法里频带边界的频率点不能简单算到某一个频带里否则带间会重叠或泄漏。我建议用严格的大小判断例如idx f fl f fh这样相邻频带在边界处只归入一边。为了更严谨可以把边界频点对半拆给两个频带但工程上影响很小一般不做也行。5.4 批量处理时的脚本组织建议如果有一段多工况数据我建议把核心处理封装成函数输入是x、fs输出是fc和LpA。再写一个外层循环读取不同文件最后把结果汇总到表格里。这样可以避免每个文件都重跑一遍 FFT 和频带循环。另外.mat文件记得同时保存采样率fs很多丢失校准信息的问题都是采样率对不上导致的。个人体会是自己写的脚本做研发分析非常够用但不要盲目挑战标准合规。如果项目需要正式计量报告还是找专业声学仪器对照验证用手持声级计在实验室里打一遍同样信号结果在误差范围内再大规模使用脚本。这也是我目前的工作流先用脚本快速扫一遍所有数据选出超标的几个工况再用标准声学系统复测确认。这样既省时间又不至于在关键数据上翻车。本文还有配套的精品资源点击获取
返回列表