在办公室做过一次盲听评价同一台冰箱放在不同的减振垫上用声级计测出来的A声级几乎没差别但房间里的每个人都能毫不犹豫地告诉你有一台明显更烦人。有人形容那个声音沙沙的毛躁的有人说是像有什么东西一直在摩擦。这种仪器读数和主观感受对不上的情况在噪声工程里太常见了——也正是声品质存在的意义。声品质Sound Quality不是简单的声音好坏它把人对声音的主观感受拆解成可计算的客观参数其中**粗糙度Roughness**专门描述声音中快速起伏带来的毛刺感。配合MATLAB我们可以把这种感受变成一条曲线、一组数值甚至进一步定位到具体频段反推是哪个零部件在捣乱。这篇文章我会把粗糙度的原理、MATLAB实现方法、参数调优经验以及如何把它放进声品质综合分析框架里从头到尾讲清楚。适合看这篇文章的主要是三类人做汽车NVH、家电噪声整改的工程师想用客观指标量化听感的音频算法工程师以及声学方向做项目需要写代码处理数据的学生。不管你属于哪一类按下面的链路走一遍基本能搭出自己的一套粗糙度分析工具。1. 粗糙度到底在衡量耳朵的什么感觉1.1 把毛躁感拆成两个物理量人耳听到的粗糙本质上来源于声音在时间轴上的快速调制。想象一个纯音如果它的响度在短时间内反复起伏耳朵就会感觉这个声音不干净像打磨不光滑的表面。这种起伏有两个关键参数一是多久起伏一次也就是调制频率二是起伏的幅度有多大也就是调制深度。调制频率太低时人耳能分辨出一声一声的波动比如风扇叶片转动时呜——呜——的喘振感那叫波动强度Fluctuation Strength调制频率太高时人耳已经跟不上这种变化主观上反而感觉声音变宽、变亮。只有在中间的某个频段人耳会强烈感知到沙沙咔啦咔啦的粗糙感这就是粗糙度的核心敏感区。在实际工程里这个敏感区大致对应15Hz到300Hz的调制频率。其中20Hz到200Hz附近的感觉最明显调制频率再往上粗糙感会逐渐减弱但这不意味着可以完全忽略。1.2 那个1 asper到底是怎么定的粗糙度的单位是asper这个单位是Zwicker和Fastl通过大量主观听音实验定出来的。基准定义是一个1kHz的纯音声压级60dB SPL用70Hz的正弦波做100%调幅此时人耳感知到的粗糙度规定为1 asper。这个基准有两层含义。第一它选了一个中频、中响度的信号这是日常生活中最常出现的普通噪声范畴第二100%调制深度意味着声音的包络从最小到最大完全起伏对应一个明确的物理边界。有了这个基准粗糙度计算出来就不仅仅是相对大小而是可以用asper为单位做绝对判断。很多初学者会忽略一个细节粗糙度的感知是非线性的。调制深度从10%增加到20%带来的粗糙感提升和从80%增加到90%不一样。所以计算时不能简单地把调制深度做线性放大而是要经过压缩映射。这也是为什么网上有些粗糙度MATLAB代码算出来数值一会儿夸张一会儿偏低很多人就是在这里处理得太粗暴。1.3 粗糙度和波动强度同一件事的两个频段粗糙度和波动强度虽然在主观感受上一个偏沙沙的、一个偏喘气的但在算法上它们高度相似都是对信号包络做调制分析。区别在于调制频率范围不同。我习惯把这两个参数看作一对搭档波动强度管0.5到20Hz左右的慢起伏粗糙度管15到300Hz的快起伏中间的过渡区有重叠但不算冲突。在分析空调压缩机这类旋转机械噪声时低频喘振可能是波动强度主导而齿轮啮合的高频抖动则表现为粗糙度。如果只算粗糙度低频喘振的贡献会被高通滤波器滤掉一部分最后结果就是吵但粗糙度不高容易误导整改方向。1.4 响度、尖锐度、粗糙度声品质参数的各自分工刚接触声品质的人最常犯的毛病是把所有问题都往响度上推。响度Loudness描述的是多大音量单位sone尖锐度Sharpness描述的是刺不刺耳单位acum而粗糙度描述的是时间维度的细腻质感单位asper。举个例子一条平坦的粉红噪声和一条在3kHz处有尖锐峰值的等响度噪声响度可以做到很接近但尖锐度差异明显一个稳定纯音和一个做70Hz调制的声音响度和尖锐度都能保持一致但粗糙度天差地别。这就解释了为什么A声级测不出来的差异声品质参数可以敏感地捕捉到。在做产品噪声评价时这三个参数加上波动强度共同构成基础指标体系缺一个都不完整。2. 在MATLAB里实现粗糙度计算的完整链路2.1 程序主流程从WAV文件到asper值要经过哪六步我在MATLAB里实现粗糙度计算时大致把流程分成六步信号读取与前端校正、Bark尺度分频带、频带内包络提取、包络调制参数估计、分频带粗糙度加权、总粗糙度合成。这个链路里最核心的思路是先分频带再在每一个频带内独立做包络→调制谱→调制深度→加权的处理。为什么要先分频带因为人耳对粗糙度的感知不是全频谱一次性总体处理的而是先由耳蜗在不同特征频率处做滤波再对每个通道的时域波形进行包络感知最后在中枢做整合。所以MATLAB代码也必须模拟这种并行处理。2.2 前端处理采样率、分帧与Bark滤波器组第一步是采样率检查。粗糙度关心的调制频率最高到300Hz按Nyquist定理至少要600Hz采样率但实际信号的载波频率可能远高于此所以建议输入信号采样率不低于48kHz。我经常遇到有人拿8kHz采样率的电话录音去算粗糙度结果包络高频部分严重失真数值毫无意义。第二步是分帧。虽然粗糙度是准稳态参数但工程信号往往时变不能把整段信号揉在一起算一个值。我习惯把信号切成512ms的帧、75%重叠这样既能覆盖最低调制频率15Hz的两个完整周期又能在时间上平滑追踪粗糙度的变化。窗口函数一般选Hann窗减少帧边界泄漏。第三步是Bark尺度滤波器组。Bark尺度是模仿人耳临界频带划分出来的频率刻度全频谱一般分成24个临界频带。下面是我常用的前几个Bark频带边界Bark带编号频率范围Hz带宽中心趋势10–100低频区很窄2100–200线性增长3200–300线性增长4300–400线性增长5400–510开始变宽6510–630变宽7630–770变宽152000–2320进入中高频2412000–15500最高临界带在MATLAB里实现这组滤波器可以直接用butter或cheby2设计带状通滤波器组。要注意的是Bark尺度带通滤波器不建议用FIR直接实现因为低频带0-100Hz要求非常窄的带宽FIR需要极高阶数才能达到足够陡峭的过渡带计算量很大。IIR滤波器更高效相位问题通过filtfilt零相位滤波解决。2.3 包络提取与调制参数估计希尔伯特不是唯一选择得到每个Bark带的时域信号后下一步是提取包络。最常用的手段是希尔伯特变换MATLAB里就是abs(hilbert(xb))。希尔伯特包络在平稳信号下表现很好但对于瞬态冲击较多的噪声比如敲击声、点火声包络会带有明显的高频残迹必须先做平滑或低通滤波再分析。我实测下来一个更稳的做法是对包络再做一次15Hz到300Hz的带通滤波把缓变的趋势成分和更高频的杂散都滤掉。这样处理后的包络信号能量主要集中在调制频率范围后续做频谱分析才干净。调制度的估计有两种思路。一种是在时域直接用包络的均值与方差来估算比如[ m \frac{E_{\max} - E_{\min}}{E_{\max} E_{\min}} ]另一种更精确的做法是取包络做FFT找到调制频率处的谱峰幅度与频谱直流分量的比值。两种方法在小调制深度时结果接近但深度大、波形畸变时会有偏差。我在代码里通常先做FFT找主峰值频率f_mod再用时域峰值法算调制指数m同时得到频率和深度两个参数。2.4 分频带加权与总粗糙度合成算法公式与代码骨架每个Bark频带的粗糙度贡献可以按简化Zwicker模型来算。核心关系是该频带的粗糙度大致正比于调制频率f_mod和该频带内的调制深度以dB计再乘一个该频带的权重系数。常见参考公式形式为[ R C \cdot \sum_{i1}^{24} g_i \cdot \Delta L_i \cdot \frac{f_{\text{mod},i}}{1000} ]其中ΔL_i是第i个Bark带内时变响应的调制深度单位dBf_mod,i是该频带内的调制频率单位Hzg_i是考虑临界频带间掩蔽效应的权重。系数C用于把基准信号1kHz、70Hz纯音调制、100%深度校准到1 asper。这里的关键是ΔL要尽量在响度域计算而不是直接拿声压包络算。因为人耳对声压存在压缩非线性同样20dB的声压起伏在高声级下的主观起伏感比低声级下小。用特征响度代替声压包络得到的ΔL能更好地反映真实感知。MATLAB中可以用Zwicker响度模型需要外中耳传递函数修正分频带计算特征响度再做包络分析。如果项目周期紧也可以先用声压包络近似但最后校准一定要做。下面给一个精简的算法骨架重点是让你理解每一步怎么衔接不是直接抄了就跑。真正工程化时你需要把滤波器组、包络处理等细节做扎实function R_total rough_estimate(x, fs) % 1. 分帧参数 frameLen round(0.512 * fs); overlap round(frameLen * 0.75); nFrame floor((length(x) - frameLen) / overlap) 1; % 2. Bark频带边界(Hz)这里列出部分示意 barkEdges [0 100 200 300 400 510 630 770 920 1080 ... 1270 1480 1720 2000 2320 2700 3150 3700 ... 4400 5300 6400 7700 9500 12000 15500]; nBark length(barkEdges) - 1; % 3. 逐帧计算 R_frame zeros(nFrame, 1); for k 1:nBark [b, a] butter(4, [barkEdges(k) barkEdges(k1)]/(fs/2), bandpass); xb filtfilt(b, a, x); env abs(hilbert(xb)); % 包络 env filtfilt(bp15_300b, bp15_300a, env); % 15-300Hz带通 for n 1:nFrame seg env(n*overlap1 : n*overlapframeLen); % 主调制频率提取 [f_mod, m_depth] modulation_para(seg, fs); % 分频带粗糙度加权(权重gi需根据频带和模型预先设定) R_frame(n) R_frame(n) g(k) * m_depth * f_mod / 1000; end end R_total C * R_frame; % 需要根据基准信号校准C end需要注意上面代码里的g(k)和C不能随便填。标准做法是用一组已知粗糙度的标准信号1kHz、70Hz调制的基准声去反向标定C这样整条链路才具备横向对比的意义。3. 调参与误差控制为什么不同工具算出来总差一截3.1 采样率和滤波器阶数算得准的第一道关粗糙度计算对前端参数极其敏感。采样率不够时频带滤波器在接近Nyquist频率处会产生畸变滤波器阶数不够时相邻Bark带之间发生能量泄漏尤其是低频带窄带宽部分泄漏会直接污染包络调制度。我自己的经验是每个Bark带滤波器的阻带衰减至少做到60dB以上过渡带尽量收窄。如果一个频带的滤波器阶数已经很高还是不够干净与其继续加阶数不如考虑把滤波器拆成两段级联或者用designfilt做等纹波IIR设计。实测中4阶Butterworth在某些窄带上不够陡换成6阶Chebyshev II型效果更好代价是相位非线性但filtfilt能很大程度缓解。3.2 包络算法的差异检波、希尔伯特还是响度域网上能搜到的粗糙度MATLAB代码至少有三种包络处理方式简单半波整流低通、希尔伯特包络、响度域包络。这三种方式计算出的调制指数差异可能超过20%尤其对于调幅深度比较大、波形上下不对称的信号。我的建议是做科研对比或产品评价时统一采用响度域包络。虽然计算量更大但这是最接近人耳真实路径的做法。如果只是为了快速排查哪个频带在产生粗糙感用希尔伯特包络就够了但要在报告里注明算法版本避免和其他工具直接比绝对值。另外要特别小心低频Bark带比如0-100Hz的包络。这个频带里的载波本身就只有几十赫兹和调制频率范围高度重叠包络提取后几乎分不清载波和包络算出来的调制频率和深度都不可靠。我通常会对前几个低频带的粗糙度贡献单独查看如果发现异常高先确认是不是滤波器或包络泄漏再做结论。3.3 与商用软件的数值对齐别追求逐位相同很多工程师拿着自己的MATLAB代码去和Head Acoustics ArtemiS、Siemens LMS等商业软件对数据发现R值总有20%到30%的偏差然后就慌了。其实这不一定是代码错了而是商业软件往往实现了完整的Zwicker响度层计算包含外中耳传递函数、特征响度压缩、时变响度的特定平均方式等你的简化模型没有这些环节差异自然存在。正确的对齐方式不是追求绝对值一致而是先固定一组基准信号做线性校准。比如拿5个已知粗糙度范围的信号算出自己代码的结果和商业软件结果的比值如果比值基本恒定说明算法趋势一致只是增益不同标定一个偏移系数就能对齐。如果比值忽大忽小才需要回头检查调制深度估计或权重系数是否有结构性错误。3.4 边界频率附近的人工伪迹处理15Hz高通和300Hz低通是粗糙度分析的边界。IIR滤波器的相位特性在边界频率附近会带来瞬态畸变如果直接作用于包络会伪造出假调制。解决方法是一用filtfilt零相位滤波但这个函数对长信号内存占用较高分段处理时一定要有足够的前后重叠二滤波之前先对包络做边缘延拓如对称延拓100个采样点算完再裁掉。这两个技巧细节虽然简单但能让边界频率附近的粗糙度曲线稳定很多。4. 声品质应用分析当粗糙度成为产品决策依据4.1 典型场景判据从感觉不对到数值超标粗糙度最常用于产品噪声的横向对比和合格判定。比如我评估过的一款车载空调压缩机怠速工况下正常状态粗糙度在0.3 asper左右更换一个供应商后同一工况直接跳到0.9 asper主观听感确实出现了明显的嗒嗒声。这就是客观参数和主观感受的一致性体现。不过不要把这些数值当金科玉律。不同算法版本、不同麦克风测量位置同一台机器的粗糙度会有偏差。比较科学的做法是企业先定义自己的标准测量工况标准算法版本标准听音人员采集一批可接受和不可接受的样件数据用统计方法定出内部阈值。4.2 多维声品质参数如何一起用粗糙度很少单独作为决策依据我通常会把响度、尖锐度、粗糙度、波动强度放在一起看。主观抱怨最可能相关的客观参数太吵响度Loudness刺耳啸叫尖锐度Sharpness沙沙的磨砂感粗糙度Roughness一阵一阵喘波动强度Fluctuation Strength在MATLAB里可以依次算出这四个参数然后做一个雷达图或者二维散点图。比如用fitlm把主观打分和这四个参数做多元线性回归看哪些参数对主观感受的贡献大。不过要注意响度和粗糙度在某些机械噪声中相关性很高一起放进回归模型会产生共线性问题。遇到这种情况先用corrplot看看参数间的相关系数或者用主成分分析降维后再回归。4.3 贡献谱分析定位粗糙感来自哪个频带粗糙度计算过程中每个Bark频带的贡献值其实都保留了这给了我们一个非常实用的分析维度画出各频带粗糙度贡献柱状图或伪彩图一眼就能看出粗糙感主要来自哪个频段。有一次我分析一个电动工具的噪声总粗糙度不算特别高但主观上明显觉得中频很毛。一看贡献谱4到8 Bark约300到1270Hz一段的贡献显著高于其他频带。结合谐波阶次分析最后定位到电机换向器和碳刷摩擦的激励频率上。这就是贡献谱的价值它把粗糙感这条笼统的抱怨变成了可执行的整改线索。5. 我踩过的几个坑以及给新手的复现建议5.1 时窗参数与统计量选择刚开始做时我用256ms帧长算整个文件发现结果波动很大同一段录音不同帧之间粗糙度能从0.2跳到1.2。后来把帧长加大到512ms重叠提高到75%结果曲线平滑很多但代价是时间分辨率下降。对非稳态噪声我建议输出时间-粗糙度曲线同时统计整个时段的N5值即超过95%时间的噪声数值作为代表值不要只取平均值。5.2 用标准测试信号自检写粗糙度算法时一定要用基准信号自检。最简单的方法生成一段1kHz、60dB、70Hz调制深度100%的正弦信号期望算出1 asper。如果算出来不是1就调整C系数。另一个自检手段是对一个恒定纯音做0%、10%、50%、100%四种调制深度看粗糙度是否单调上升且上升幅度符合主观趋势。这两关过了算法才算能拿去处理真实信号。5.3 不要把不同算法的绝对值直接横向对比我在项目里吃过一次亏早期用简化算法得出某产品粗糙度0.5 asper后期改用响度域算法变成0.75 asper如果不知道前后差异很容易得出噪声变差了的错误结论。所以在声品质分析项目里算法版本和参数设置必须冻结任何对比都要在同一套流程下完成。我习惯把算法版本、滤波器参数、帧长、重叠率都写进结果Excel的元数据里防止隔了几个月自己都忘了当初怎么算的。5.4 复现建议先跑流程再调精度给新手的复现路径是第一步先用希尔伯特包络分段加权跑通全流程得到一个趋势正确的粗糙度值第二步在4个Bark频带比如低频带、中低频、中高频、最高频上做详细包络和基准确认第三步引入响度域处理和外中耳修正校准绝对值。不要一上来就照着商业软件实现全套Zwicker模型容易淹没在细节里。推荐顺手把这几个算法也看一眼Sottek的听觉模型Hearing Model、Moore的响度模型、以及ISO 532-1对响度计算的规定。粗糙度虽然没有对应ISO标准但响度域的处理方式基本承接自这些标准理解之后再看文献里的粗糙度公式就不会发懵。最后再分享一个我在实际项目里的体会粗糙度数值不要脱离信号本身去解读。同一个0.6 asper出现在空调吹风噪声里可能觉得轻微出现在方向盘后面的异响里可能就非常恼人。做工程判断时一定要把粗糙度、贡献谱、主观回放三者放在一起看。算法给你的是一个线索而不是结论。带上这套MATLAB分析工具你至少能少走我当初绕过的那些弯路——从感觉这声音怪到知道它怪在哪、怪在哪个频带这个跨越不值钱但很实用。