简介这套声品质分析工具包聚焦响度与尖锐度两类客观参数依据ISO 532-1:1991与Fastl方法实现适用于汽车、家电、工业设备等领域噪声主观感知量化研究。包内共14个文件以MATLAB核心算法与滤波器设计脚本为主配有说明文档、测试音频及运行辅助文件压缩包约484KB结构清晰可直接用于时域声压或频谱数据计算。目前已有15人学习下载。借助该工具读者可快速得到响度时间历程或单值sone与尖锐度acum结果用于声品质主客观关联建模、优化验证及基础教学是一套轻量实用的入门与进阶工具。 去年接了个新能源汽车电机啸叫的优化项目客户给的验收指标不是简单的“多少分贝”而是两个字——主观感受。坐在副驾听加速声A计权声压级大家都差不多但有些车就是尖锐刺耳有些车闷得让人心烦。研究声品质的人都知道这时候真正能说明问题的物理量是响度Loudness和尖锐度Sharpness。手头没有商用声品质测试软件授权我只好用MATLAB搭了一个快速计算工具包。做完之后发现同一段音频文件算出来的结果和实验室专业设备挺接近完全能用于工程预判和批量筛选。这篇博文就把这套工具包的实现思路、关键代码和踩过的坑完整梳理一遍适合正在做声品质分析、MATLAB做声音处理以及想搞懂响度尖锐度怎么落的同学。1. 先搞清楚响度和尖锐度到底在算什么很多人第一次接触响度会把它和声压级搞混。声压级是物理客观量用dB表示响度是人耳对声音“有多大”的主观量单位是sone。定义是这样的1kHz、40dB SPL的纯音响度是1 sone人耳感觉响2倍的声音响度就是2 sone。所以响度是和感受挂钩的不是简单把dB值取个指数就能画等号。尖锐度则描述声音“尖锐”或者“刺耳”的程度单位是acum。同样是60dB的噪声能量都集中在低频听起来闷闷的尖锐度就低能量集中在高频听起来刺刺的尖锐度就高。工程上常见的电机啸叫、齿轮啸叫、高压气动噪声尖锐度往往偏高这正是用户烦躁感的重要来源。计算响度和尖锐度目前工业界默认的底层模型是Zwicker模型对应标准是ISO 532B。这个模型把人耳看成一组带通滤波器滤波器中心频率不是线性排列而是按临界频带Critical Band分布。人耳能分辨的频率范围大概分成24个临界频带单位是Bark1 Bark代表一个临界频带宽度。Zwicker模型的做法是先把声音信号通过这24个频带的等效滤波网络再对每个频带计算“特定响度”Specific Loudness单位是sone/Bark最后把24个频带的特定响度积分起来就得到总响度。尖锐度的计算建立在响度基础上比较常用的是Zwicker和Aures两种加权法。我工具包里用的是Zwicker/Aures结合的思路公式可以写成[ S 0.11 \times \frac{\int_{0}^{24} N(z) \cdot g(z) \cdot z , dz}{N} ]其中 (N(z)) 是特定响度(z) 是以Bark为单位的临界频带序号(N) 是总响度(g(z)) 是加权函数。这个加权函数在低频段接近1到了高频段会急剧放大这样尖锐度才能突出高频成分的主观影响。理解到这个层面后面写代码就清楚多了先算响度分布再算加权积分整个过程就是一套滤波加积分。2. 工具包整体设计与计算流程这个工具包的设计目标很明确输入一段音频输出响度和尖锐度两个值同时尽可能快。市面上已经有商用的声品质分析软件但我需要能嵌入现有MATLAB处理流程里能批量跑上百个音频文件还能把中间结果拿出来画图分析。所以自己实现一套是可取的。整体流程分五步第一步读取音频并重采样。响度计算对不同采样率的处理略有不同我统一把输入音频重采样到48kHz因为Zwicker模型里的外耳中耳修正滤波器系数表本身就是按48kHz设计的。用MATLAB的resample函数即可重采样前最好先用anti-alias低通滤波否则高频伪影会影响尖锐度。第二步声压级校准。从音频文件读出来的数字信号没有绝对声压参考必须知道数字满幅对应的声压级。理想情况是录音时录过校准信号比如94dB SPL、1kHz的标准声源。如果没有校准信息只能按文件头的灵敏度或者默认值估这时算出的响度是相对值只能做同一批数据内的横向对比。第三步外耳中耳修正。Zwicker模型里面包含人体外耳道和中耳对声音的传递修正这一环节直接影响高频频段的权重尤其是尖锐度。实现上可以用一组固定系数的IIR/FIR滤波器完成也可以直接在频域做加权我用的是MATLAB里的filtfilt做零相位滤波避免相位失真影响瞬时响度。第四步24个临界频带分解。工程上做快速计算不需要严格复刻Zwicker的专用滤波器组可以用1/3倍频程滤波器组近似。1/3倍频程的中心频率和Bark临界频带大致有对应关系虽然每个Bark不一定正好等于一个1/3倍频程但对于工程筛选来说误差可以接受。MATLAB的Audio Toolbox提供了octaveFilter和filterBank可以直接生成1/3倍频程滤波器组但如果没装Audio Toolbox也可以用designfilt自己设计一组带通滤波器。第五步响度积分与锐度加权。每个频带分解出来后先按Zwicker的特定响度转换公式把频带能量转成特定响度然后对频带轴积分得到总响度再按尖锐度加权公式计算acum值。工具包模块划分大概是这样的loadAudioFile.m负责读取、重采样、校准。outerMiddleEarFilter.m外耳中耳修正。thirdOctaveFilterBank.m生成1/3倍频程滤波器组。computeSpecificLoudness.m频带能量到特定响度的转换。computeLoudnessSharpness.m主函数整合上面所有模块。这种设计的核心思路是模块化每个环节都可以单独验证。比如先单独跑thirdOctaveFilterBank用扫频信号确认每个滤波器中心频率和增益正确再去跑后面的积分不然后面算出的数字有问题时根本定位不到是哪一步出的错。3. 核心函数与关键代码的实现细节下面给出主函数骨架这个结构基本可以直接搬用function [N, S] loudnessSharpness(x, fs, splCalib) % x: 单声道音频信号double类型范围[-1, 1] % fs: 采样率 % splCalib: 数字满幅对应的声压级(dB SPL)没有校准信息时可设0此时输出为相对值 % N: 总响度(sone)S: 尖锐度(acum) % 1. 重采样到48kHz fsTarget 48000; if fs ~ fsTarget x resample(x, fsTarget, fs); fs fsTarget; end % 2. 声压级校准 if nargin 3 || isempty(splCalib) splCalib 0; end fullScaleSPL splCalib; xPa x * 10^((fullScaleSPL - 94) / 20) * 2e-5; % 粗略换算到Pa % 3. 外耳中耳修正 xFiltered outerMiddleEarFilter(xPa, fs); % 4. 生成1/3倍频程滤波器组并分解 bands thirdOctaveBands(); Nbands length(bands); specLoudness zeros(1, Nbands); for i 1:Nbands [b, a] thirdOctaveFilter(bands(i), fs); y filtfilt(b, a, xFiltered); % 频带有效值转换成特定响度 rmsVal sqrt(mean(y.^2)); specLoudness(i) specificLoudness(rmsVal, bands(i)); end % 5. 总响度把特定响度在Bark轴上积分 barkAxis 0.5:1:23.5; % 24个频带的近似Bark中心 N trapz(barkAxis, specLoudness); % 6. 尖锐度加权积分 g zeros(1, Nbands); for i 1:Nbands z barkAxis(i); if z 16 g(i) 1; else g(i) 0.15 * exp(0.42 * (z - 16)) 0.85; end end S 0.11 * trapz(barkAxis, specLoudness .* g .* barkAxis) / N; end注意这里面specificLoudness函数要按Zwicker的模型做非线性压缩。低频段和高频段的特定响度增益不一样具体可以使用经典Loudness模型的查表法或者用近似公式。实操中比较稳妥的做法是直接采用标准里给出的自由场/扩散场外耳中耳修正表以及特定响度压缩系数按频带插值。还有一个细节如果输入音频时间很长比如几十秒的噪声样本直接用filtfilt跑全部数据会占用大量内存而且逐频带跑循环会有点慢。我实际优化的时候是先把信号切成若干帧每帧2秒加50%重叠逐帧计算响度再做时域平均。这样既省内存又方便分析响度随时间的变化画成“响度-时间曲线”更直观。代码里需要解释的另一个重点是xPa的换算。这个写法是用参考声压2e-5 Pa算绝对声压如果只是相对比较可以把校准值设为0这样算出来的响度值不会准确落在sone标尺上但横向排序完全没问题。我做批量筛选时常用相对模式谁比谁吵谁比谁刺耳排序结果可信。4. 参数设置、校准与容易踩的坑这个工具包写完之后最难的不是算法本身而是参数校准和边界条件处理。有几个坑我印象特别深。第一个坑滤波器组的中心频率和Bark轴对不上。最初我直接用标准的1/3倍频程中心频率表从25Hz到20kHz一共31个频带然后硬映射到24个Bark结果边界频带的响度计算偏大。后来我改成按Bark轴重新排列频带低频段部分Bark对应多个1/3倍频程频带高频段则一个频带对应一个Bark积分起来才比较稳。建议在做积分之前先把频带-临界频带对应表画出来确认没有错位。第二个坑校准环节最容易被忽略。如果输入信号是用声卡直接录的数字满幅对应的声压级不是随便猜的必须用标准声级校准器标定。我用的是94dB SPL/1kHz校准声源先把校准信号录进MATLAB计算满幅下的RMS反推fullScaleSPL。不同声卡的输入增益不同每次换设备都要重新标定否则精度根本谈不上。第三个坑filtfilt虽然能避免相位失真但不能直接用于高频段滤波器因为高频段带通滤波后信号经常包含很大的瞬态filtfilt的边缘效应会把响度算高。处理办法是给信号做前后各0.5秒的镜像延拓滤波后再截掉延拓段这样边缘效应基本消除。第四个坑尖锐度对高频成分极其敏感。同样一段音频如果录音设备本底噪声很大高频噪声会把尖锐度值整体抬高。实际使用时要先看频谱确认被分析频段是目标声源主导不是录音设备底噪主导。对于底噪超标的文件建议先做250Hz到10kHz之间的带通预滤波再进响度尖锐度计算流程。这些参数问题我用一个表整理过写在这里可以直接对照排查现象可能原因处理方式总响度整体偏大校准值设置错误用标准声源重新标定响度有波动但不稳定滤波器组中心频率/Bark映射错误打印频带与Bark对应表检查错位尖锐度异常高录音底噪/预滤波未做加带通滤波器处理高频段能量过高外耳中耳修正未生效检查滤波器系数是否按48kHz设计长音频计算过慢循环过多/无预延拓分帧计算并行化5. 实测对比与使用心得工具包跑通之后我用三组数据做了验证。第一组是1kHz、40dB SPL纯音理论响度应该在1 sone左右。我的程序跑出来是1.08 sone误差8%在工程允许范围内。第二组是一段白噪声响度应该在约11 sone附近程序结果落在10.5到11.2之间。第三组是真实电机啸叫音频响度17.3 sone尖锐度2.9 acum和商用声品质终端测出的17.1 sone、2.8 acum很接近。说明这套流程虽然做了1/3倍频程近似但结果依然可参考。批量处理的时候我习惯把每个文件的响度和尖锐度打成一个表格同时输出特定响度曲线。这样在筛选候选方案时看的不是一堆孤立的数字而是每个频带对响度和尖锐度的贡献量。有一次做空调风噪优化对比两个设计方案的响度相近但尖锐度一个1.8 acum一个2.3 acum明显后者更刺耳主观试听也验证了这一点。这就是单独看分贝值发现不了的问题响度和尖锐度组合起来才完整。我个人在实际使用中的体会是声品质分析工具包这种项目如果只是做单次工程评估安装现成软件最快。但如果要做批量筛选、算法对比、研究不同频带对主观感受的贡献自己用MATLAB写一套反而更灵活。关键是校准和数据预处理不能省否则计算出花来结果也不可靠。针对后续扩展这套工具包还能加粗糙度、抖动度等指标核心框架不变只是多接几个积分分支而已值得持续维护。本文还有配套的精品资源点击获取