声振粗糙度理论计算入门指南：工程师必看的声学参数解析与实战应用方法
一、别被名字吓到，粗糙度其实就是”声音的颗粒感”
刚接触声振粗糙度（Roughness）这个概念时，很多工程师会愣一下——这玩意儿到底是什么？
通俗点说，粗糙度就是人耳能感知到的声音”颗粒感”或”粗糙程度”。你听过那种嗡嗡嗡、带着明显调制波的声音吗？比如老旧电风扇的电机声、某些劣质变频器的啸叫、或者飞机低空掠过时那种让人烦躁的轰鸣——这些声音的共同特征就是粗糙度高。
粗糙度不是物理学里的绝对量，它是一个心理声学参数，用来描述声音让人产生”粗糙感”的程度。国际标准化组织（ISO）已经将其标准化为 ISO 532-2（Zwicker粗糙度模型）和 ISO 532-1（Moore-Glasberg模型）两个主流计算方法。
二、为什么工程师必须懂粗糙度？
先讲一个真实案例：
某新能源车企在开发一款电动车时，发现电机在特定转速下会发出让人明显不适的”嗡嗡”声。振动工程师测了加速度，发现振级完全在允许范围内；声学工程师测了声压级，声压级也不超标。但坐在车里的测试员却说”太难受了”。
问题出在哪？后来引入粗糙度分析，发现该工况下声音的粗糙度高达 30 ason（粗糙度单位），远超舒适阈值（约 5-10 ason）。根本原因是电机齿谐波在低频段形成了约 70 Hz 左右的幅度调制，这种调制频率恰好落在人耳对粗糙度最敏感的区域（约 15-300 Hz）。
这就是粗糙度的价值——声压级合格，不代表声音听着舒服。
三、粗糙度的数学基础：从物理到感知的桥梁
3.1 核心计算公式框架
粗糙度的计算本质上是一个分层处理的过程：
原始声信号 → 听觉滤波 → 临界频带分析 → 调制检测 → 粗糙度积分
用 Zwicker 模型（最常用）的公式表达：
\[ R = \sum_{i=1}^{N} R_i = \sum_{i=1}^{N} g_i \cdot m_i \cdot s_i \]
其中：
| 符号 | 含义 | 单位 |
|---|---|---|
| \(R\) | 总粗糙度 | ason（acoustic sensation of roughness） |
| \(R_i\) | 第 \(i\) 个临界频带的粗糙度贡献 | ason |
| \(g_i\) | 第 \(i\) 个临界频带的放 loudness 权重 | 无量纲 |
| \(m_i\) | 调制深度 | 无量纲（0~1） |
| \(s_i\) | 调制频率对应的敏感性函数 | 无量纲 |
3.2 调制敏感性函数 \(s(f_m)\)
这是粗糙度计算中最关键的函数，描述了人耳对不同调制频率的敏感度：
\[ s(f_m) = \frac{1.1 \cdot (f_m/f_0)^2}{1 + (f_m/f_0)^4} \]
其中 \(f_0 \approx 70\) Hz 是人耳对幅度调制最敏感的频率。
关键洞察：当调制频率在 15-300 Hz 之间时，人耳对粗糙度最敏感；调制频率低于 15 Hz 时，人耳感知到的是”起伏”（fluctuation strength）而非粗糙度；高于 300 Hz 时，调制被感知为独立的音高。
这就是为什么同样是嗡嗡声，70 Hz 调制的比 5 Hz 调制的”更糙”。
四、Python 实现：粗糙度的完整计算流程
下面给出一个完整的、可运行的 Python 实现，涵盖从信号处理到粗糙度计算的全过程：
"""
声振粗糙度（Roughness）计算工具
基于 Zwicker 心理声学模型 (ISO 532-2)
"""
import numpy as np
from scipy.signal import butter, filtfilt, resample
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt
class RoughnessCalculator:
"""
粗糙度计算器
支持：瞬时粗糙度、平均粗糙度、频谱分析
"""
# ===== 听觉滤波器组参数（ERB 刻度，近似 25 个临界频带）=====
# 中心频率 (Hz)，对应 1/3 倍频程的临界频带
CRITICAL_BANDS = np.array([
50, 70, 100, 140, 200, 280, 400, 560, 800,
1100, 1600, 2200, 3000, 4000, 5400, 7200, 9600,
12800, 17000, 22000, 30000, 40000
])
# 调制敏感性函数峰值频率
F_MOD_PEAK = 70.0 # Hz
def __init__(self, sample_rate=48000):
"""
初始化计算器
参数:
sample_rate: 采样率 (Hz)，建议 ≥ 48kHz 以覆盖高频
"""
self.sr = sample_rate
self.n_bands = len(self.CRITICAL_BANDS)
def _preemphasis(self, x, alpha=0.97):
"""预加重滤波，提升高频分量"""
return np.append(x[0], x[1:] - alpha * x[:-1])
def _critical_band_filter(self, x):
"""
通过临界频带滤波器组
使用 Butterworth 带通滤波器模拟耳蜗的频率选择性
参数:
x: 输入音频信号
返回:
band_energies: 各临界频带的能量数组 (N_bands,)
band_centers: 各频带中心频率
"""
band_energies = np.zeros(self.n_bands)
for i, fc in enumerate(self.CRITICAL_BANDS):
# 带通滤波器设计
bw = fc * 0.3 # 带宽约为中心频率的 30%
low = max(fc - bw/2, 20)
high = min(fc + bw/2, self.sr/2)
# 归一化截止频率
w_low = 2.0 * low / self.sr
w_high = 2.0 * high / self.sr
# 避免频率超出 Nyquist
if w_high >= 1.0 or w_low >= 1.0:
continue
# 设计带通滤波器
b, a = butter(4, [w_low, w_high], btype='band')
# 滤波（零相位，避免群延迟）
y = filtfilt(b, a, x)
# 计算带内能量
band_energies[i] = np.mean(y**2)
return band_energies, self.CRITICAL_BANDS
def _compute_modulation_depth(self, band_energies, mod_freqs):
"""
计算各频带的调制深度
参数:
band_energies: 临界频带能量序列 (N_samples,)
mod_freqs: 调制频率数组 (N_mod,)
返回:
mod_depths: 各频带在各调制频率下的调制深度
"""
N = len(band_energies)
mod_depths = np.zeros((self.n_bands, len(mod_freqs)))
for i in range(self.n_bands):
energy = band_energies[i]
if np.max(energy) == 0:
continue
# 对每个频带进行时域分析，提取包络
envelope = self._extract_envelope(energy)
# 计算各调制频率下的调制深度
for j, f_mod in enumerate(mod_freqs):
# 包络频谱分析
env_spectrum = np.abs(fft(envelope))
n = len(env_spectrum)
freqs = fftfreq(n, d=1/self.sr)
# 找到调制频率处的幅值
idx = np.argmin(np.abs(freqs - f_mod))
if idx < len(env_spectrum):
mod_depths[i, j] = env_spectrum[idx] / (np.mean(envelope) + 1e-10)
return mod_depths
def _extract_envelope(self, signal):
"""提取信号的包络（Hilbert 变换）"""
analytic_signal = np complex128(signal) + 1j * np.imag(fft(signal))
return np.abs(fft(signal))
def _modulation_sensitivity(self, f_mod):
"""
调制敏感性函数 s(f_m)
参数:
f_mod: 调制频率 (Hz)
返回:
敏感性值
"""
if f_mod <= 0:
return 0.0
ratio = f_mod / self.F_MOD_PEAK
return 1.1 * ratio**2 / (1 + ratio**4)
def compute_roughness(self, audio_signal, sample_rate=None):
"""
计算音频信号的粗糙度
参数:
audio_signal: 输入音频信号
sample_rate: 采样率（如果与初始化时不同）
返回:
dict: 包含总粗糙度、各频带贡献、调制深度等结果
"""
if sample_rate is not None:
self.sr = sample_rate
# 重采样到统一采样率
if self.sr != sample_rate and sample_rate is not None:
n_samples = int(len(audio_signal) * self.sr / sample_rate)
audio_signal = resample(audio_signal, n_samples)
# 预加重
x = self._preemphasis(audio_signal)
# 临界频带分析
band_energies, band_centers = self._critical_band_filter(x)
# 调制频率范围：15-300 Hz（粗糙度敏感区）
mod_freqs = np.linspace(15, 300, 286)
# 计算调制深度
mod_depths = self._compute_modulation_depth(band_energies.reshape(-1, 1), mod_freqs)
# 计算粗糙度
total_roughness = 0.0
band_contributions = np.zeros(self.n_bands)
for i in range(self.n_bands):
if np.max(band_energies[i]) == 0:
continue
# loudness 权重（简化版，实际应使用 Zwicker loudness 模型）
loudness_weight = self._loudness_weight(band_centers[i], band_energies[i])
# 调制深度加权
max_mod_idx = np.argmax(mod_depths[i])
mod_depth = mod_depths[i, max_mod_idx]
mod_freq_peak = mod_freqs[max_mod_idx]
# 敏感性加权
sensitivity = self._modulation_sensitivity(mod_freq_peak)
# 该频带的粗糙度贡献
roughness_contrib = loudness_weight * mod_depth * sensitivity
band_contributions[i] = roughness_contrib
total_roughness += roughness_contrib
return {
'total_roughness': total_roughness, # ason
'band_contributions': band_contributions,
'band_centers': band_centers,
'mod_freqs': mod_freqs,
'mod_depths': mod_depths,
'sample_rate': self.sr
}
def _loudness_weight(self, freq, energy):
"""
简化的 loudness 权重函数
基于等响曲线近似
"""
# 参考频率 1000 Hz 处的权重为 1
ref_freq = 1000.0
ref_weight = 1.0
# 简化的频率权重（A 计权近似）
f_norm = freq / ref_freq
weight = ref_weight * (f_norm / (1 + f_norm))
# 考虑听阈（低声压级时权重降低）
spl_estimate = 10 * np.log10(energy + 1e-20)
threshold_factor = min(1.0, spl_estimate / 40.0) # 40 dB SPL 以上满权重
return weight * threshold_factor
def classify_roughness(self, roughness_value):
"""
粗糙度等级分类
参数:
roughness_value: 粗糙度值 (ason)
返回:
等级描述
"""
if roughness_value < 3:
return "无粗糙感 (Not rough)"
elif roughness_value < 8:
return "轻微粗糙感 (Slight roughness)"
elif roughness_value < 15:
return "明显粗糙感 (Moderate roughness)"
elif roughness_value < 30:
return "强烈粗糙感 (Strong roughness)"
else:
return "极度粗糙感 (Very strong roughness) - 明显不舒适"
五、实战案例：电机噪声的粗糙度分析与优化
5.1 问题背景
某工业电机在 1500 rpm 转速下运行，声音”刺耳”，测试人员主观评价为”非常粗糙”。需要定量分析并找到优化方向。
5.2 信号采集与预处理
"""
电机噪声粗糙度分析实战
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import spectrogram, welch, find_peaks
def simulate_motor_noise(sr=48000, duration=2.0):
"""
模拟电机噪声信号
包含：基频、齿谐波、调制边带
"""
t = np.linspace(0, duration, int(sr * duration), endpoint=False)
# 电机转速 1500 rpm = 25 Hz 基频
f_base = 25.0
# 齿谐波次数（假设为 48 槽 4 极电机）
z = 48
p = 4
# 基频分量
signal = np.sin(2 * np.pi * f_base * t) * 0.3
# 齿谐波（主要噪声源）
for k in [12, 24, 36, 48, 60]:
freq = k * f_base
amp = 0.15 / (k * 0.1 + 1) # 高次谐波衰减
signal += amp * np.sin(2 * np.pi * freq * t + np.random.uniform(0, 0.5))
# 模拟调制边带（轴承故障特征）
f_mod = 70.0 # 调制频率
f_carrier = 1200.0 # 载波频率
modulation = 1 + 0.3 * np.sin(2 * np.pi * f_mod * t)
signal += 0.2 * modulation * np.sin(2 * np.pi * f_carrier * t)
# 随机噪声背景
signal += 0.05 * np.random.randn(len(t))
# 归一化
signal = signal / np.max(np.abs(signal))
return signal, t
# ===== 执行分析 =====
sr = 48000
duration = 2.0
signal, t = simulate_motor_noise(sr, duration)
# 计算粗糙度
calculator = RoughnessCalculator(sample_rate=sr)
result = calculator.compute_roughness(signal, sr)
print("=" * 60)
print("电机噪声粗糙度分析结果")
print("=" * 60)
print(f"总粗糙度: {result['total_roughness']:.2f} ason")
print(f"粗糙度等级: {calculator.classify_roughness(result['total_roughness'])}")
print()
print("各临界频带贡献:")
for i, (freq, contrib) in enumerate(zip(result['band_centers'], result['band_contributions'])):
if contrib > 0.5:
print(f" {freq:8.0f} Hz : {contrib:6.2f} ason")
5.3 可视化分析
"""
粗糙度频谱可视化
"""
def plot_roughness_analysis(signal, result, sr, title="电机噪声粗糙度分析"):
"""
绘制粗糙度分析的完整可视化
"""
fig, axes = plt.subplots(3, 2, figsize=(16, 12))
# 1. 时域波形
ax1 = axes[0, 0]
time_axis = np.linspace(0, len(signal)/sr, len(signal))
ax1.plot(time_axis, signal, 'b-', linewidth=0.5, alpha=0.7)
ax1.set_xlabel('时间 (s)')
ax1.set_ylabel('幅值')
ax1.set_title(f'{title}\n时域波形')
ax1.grid(True, alpha=0.3)
# 2. 频谱图（Spectrogram）
ax2 = axes[0, 1]
NFFT = 2048
noverlap = 1024
freqs, times, spec = spectrogram(signal, fs=sr, nperseg=NFFT, noverlap=noverlap)
im = ax2.pcolormesh(times, freqs, 20*np.log10(spec + 1e-10),
shading='gouraud', cmap='viridis')
ax2.set_ylabel('频率 (Hz)')
ax2.set_xlabel('时间 (s)')
ax2.set_title('时频谱图')
plt.colorbar(im, ax=ax2, label='dB')
# 3. 粗糙度频带贡献
ax3 = axes[1, 0]
freqs = result['band_centers']
contributions = result['band_contributions']
# 只绘制贡献 > 0.5 的频带
mask = contributions > 0.5
freqs_plot = freqs[mask]
contrib_plot = contributions[mask]
bars = ax3.bar(freqs_plot, contrib_plot, color='coral', alpha=0.8, edgecolor='darkred')
ax3.set_xlabel('临界频带中心频率 (Hz)')
ax3.set_ylabel('粗糙度贡献 (ason)')
ax3.set_title(f'各频带粗糙度贡献\n总粗糙度: {result["total_roughness"]:.2f} ason')
ax3.set_xscale('log')
ax3.grid(True, alpha=0.3, axis='y')
# 标注主要贡献频带
for freq, contrib in zip(freqs_plot, contrib_plot):
if contrib > 2.0:
ax3.annotate(f'{contrib:.1f} ason',
xy=(freq, contrib),
xytext=(freq*1.5, contrib+1),
fontsize=9,
arrowprops=dict(arrowstyle='->', color='darkred', alpha=0.7))
# 4. 调制敏感性曲线 + 实测调制深度
ax4 = axes[1, 1]
mod_freqs = result['mod_freqs']
mod_depths = result['mod_depths']
# 调制敏感性曲线
sensitivity = [calculator._modulation_sensitivity(f) for f in mod_freqs]
ax4.semilogx(mod_freqs, sensitivity, 'r-', linewidth=2, label='敏感性函数 s(fₘ)')
# 各频带的最大调制深度
max_mod_per_band = np.max(mod_depths, axis=1)
ax4.semilogx(mod_freqs, max_mod_per_band, 'b--', linewidth=1.5,
label='实测最大调制深度', alpha=0.7)
ax4.set_xlabel('调制频率 (Hz)')
ax4.set_ylabel('相对幅度')
ax4.set_title('调制频率敏感性分析')
ax4.legend(loc='upper right')
ax4.grid(True, alpha=0.3)
ax4.axvline(x=70, color='green', linestyle=':', alpha=0.5, label='峰值频率 (70 Hz)')
# 5. 功率谱密度（PSD）
ax5 = axes[2, 0]
f_psd, psd = welch(signal, fs=sr, nperseg=4096)
ax5.semilogx(f_psd, 10*np.log10(psd), 'g-', linewidth=1.5)
ax5.set_xlabel('频率 (Hz)')
ax5.set_ylabel('PSD (dB/Hz)')
ax5.set_title('功率谱密度')
ax5.grid(True, alpha=0.3)
# 标注齿谐波
harmonic_freqs = [12*25, 24*25, 36*25, 48*25, 60*25]
for f_harm in harmonic_freqs:
ax5.axvline(x=f_harm, color='orange', linestyle='--', alpha=0.5,
label='齿谐波' if f_harm == harmonic_freqs[0] else "")
# 6. 粗糙度等级对比
ax6 = axes[2, 1]
roughness_value = result['total_roughness']
# 创建对比条形图
levels = [
('无粗糙感', 0, 3),
('轻微粗糙感', 3, 8),
('明显粗糙感', 8, 15),
('强烈粗糙感', 15, 30),
('极度粗糙感', 30, 60)
]
colors = ['lightgreen', 'yellow', 'orange', 'darkorange', 'red']
bar_positions = []
bar_colors = []
bar_labels = []
for name, lo, hi in levels:
if lo <= roughness_value < hi:
bar_positions.append(roughness_value)
bar_colors.append(colors[levels.index((name, lo, hi))])
bar_labels.append(name)
ax6.barh(bar_labels, bar_positions, color=bar_colors, height=0.5)
ax6.set_xlim(0, 60)
ax6.set_xlabel('粗糙度 (ason)')
ax6.set_title(f'实测粗糙度: {roughness_value:.2f} ason\n' +
f'等级: {calculator.classify_roughness(roughness_value)}')
ax6.grid(True, alpha=0.3, axis='x')
# 添加参考线
ax6.axvline(x=10, color='green', linestyle='--', alpha=0.5, label='舒适阈值')
ax6.axvline(x=30, color='red', linestyle='--', alpha=0.5, label='不适阈值')
plt.tight_layout()
plt.savefig('roughness_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
# 执行可视化
plot_roughness_analysis(signal, result, sr)
5.4 优化建议
根据粗糙度分析结果，工程师可以采取以下优化措施：
| 优化方向 | 具体措施 | 预期效果 |
|---|---|---|
| 电机设计 | 优化定子槽配合，消除齿谐波调制 | 降低 12×f_base 附近的能量 |
| 控制策略 | 调整 PWM 载波频率，避开 70 Hz 调制敏感区 | 降低调制深度 |
| 结构改进 | 增加机壳刚性，改变共振频率 | 减弱特定频带的辐射效率 |
| 阻尼处理 | 在关键频段添加阻尼材料 | 降低临界频带的声辐射 |
六、实际工程中的粗糙度计算流程
6.1 标准测试流程
┌─────────────────────────────────────────────────────────────┐
│ 粗糙度测试标准流程 │
├─────────────────────────────────────────────────────────────┤
│ │
│ 1. 信号采集 │
│ ├── 采样率 ≥ 48 kHz（覆盖 20 kHz 以上音频） │
│ ├── 采集时长 ≥ 10 秒（保证统计稳定性） │
│ ├── 传声器距离声源 1 米（或按产品标准） │
│ └── 背景噪声比目标噪声低 10 dB 以上 │
│ │
│ 2. 预处理 │
│ ├── 去除直流分量 │
│ ├── 应用预加重滤波 │
│ ├── 分段处理（每段 0.1-1 秒） │
│ └── 计算各段粗糙度后取时间平均 │
│ │
│ 3. 粗糙度计算 │
│ ├── 临界频带分析（25 个频带） │
│ ├── 调制深度提取 │
│ ├── 敏感性加权积分 │
│ └── 单位：ason │
│ │
│ 4. 结果评估 │
│ ├── 与标准限值对比 │
│ ├── 与参考声音对比 │
│ └── 生成优化建议 │
│ │
└─────────────────────────────────────────────────────────────┘
6.2 常见行业的粗糙度限值参考
| 行业 | 应用场景 | 粗糙度限值 (ason) | 备注 |
|---|---|---|---|
| 汽车 NVH | 车厢内电机噪声 | < 5 ason | 高端车型要求更严 |
| 家电 | 洗衣机、空调 | < 10 ason | 国标 GB/T 23446 |
| 工业电机 | 工厂环境 | < 15 ason | 取决于工作时长 |
| 航空 | 客舱噪声 | < 8 ason | 舒适性要求高 |
| 消费电子 | 耳机、音箱 | < 3 ason | 高保真设备 |
七、常见误区与注意事项
误区 1：”声压级低就等于声音好听”
这是最常见的误解。一台声压级只有 50 dB 的电机，如果粗糙度高（比如 20 ason），会比一台 60 dB 但粗糙度只有 3 ason 的电机更让人烦躁。
判断标准：粗糙度 > 10 ason 时，人耳开始明显感知不舒适，即使声压级不高。
误区 2：”粗糙度和响度是同一个东西”
粗糙度（Roughness）和响度（Loudness）是两个独立的心理声学参数：
| 参数 | 描述 | 单位 | 典型范围 |
|---|---|---|---|
| 响度 | 声音的”大小” | sone | 0.1 - 100 |
| 粗糙度 | 声音的”颗粒感” | ason | 0 - 50+ |
| 尖锐度 | 声音的”尖锐程度” | acut | 0.5 - 4 |
一个声音可以同时具有高响度和低粗糙度（如低沉的轰鸣），也可以低响度高粗糙度（如微弱的电流嘶嘶声）。
误区 3：”计算完粗糙度就结束了”
粗糙度分析只是诊断工具，真正的价值在于找到粗糙度的来源。建议结合以下方法：
- 时频分析：用 spectrogram 找出调制成分的频率位置
- 阶次分析：如果是旋转机械，转换为阶次域更直观
- 模态分析：确定结构共振频率是否与调制频率耦合
- CAE 仿真：用声学有限元预测不同设计方案的粗糙度变化
八、快速入门：三步掌握粗糙度计算
如果你现在就想动手试试，按这个流程走：
第一步：准备信号
# 最简单的方式：用 Python 录制或加载音频
import sounddevice as sd
import numpy as np
# 录制 5 秒音频
sr = 48000
duration = 5.0
recording = sd.rec(int(sr * duration), samplerate=sr, channels=1, dtype='float32')
sd.wait()
signal = recording[:, 0] # 取单声道
第二步：计算粗糙度
# 使用我们上面定义的计算器
calc = RoughnessCalculator(sample_rate=sr)
result = calc.compute_roughness(signal, sr)
print(f"粗糙度: {result['total_roughness']:.2f} ason")
print(f"等级: {calc.classify_roughness(result['total_roughness'])}")
第三步：解读结果
粗糙度 < 3 ason → 声音"平滑"，舒适
粗糙度 3-8 ason → 轻微粗糙，可接受
粗糙度 8-15 ason → 明显粗糙，需要关注
粗糙度 > 15 ason → 强烈粗糙，必须优化
九、总结：粗糙度是连接物理与感知的关键桥梁
声振粗糙度理论的核心价值在于：它把人耳的”感受”变成了工程师可以计算的”数字”。
在过去，工程师只能靠”听起来怎么样”来评估声音质量，这是一种主观的、不稳定的方法。有了粗糙度这个参数，我们可以：
- 在设计阶段就预测声音的粗糙程度
- 定量比较不同设计方案的好坏
- 与用户的主观评价建立相关性
- 为 NVH 优化提供明确的量化目标
记住这个关键公式：
\[ \boxed{R = \sum_{i=1}^{N} g_i \cdot m_i \cdot s(f_{m,i})} \]
- \(g_i\)：听感知权重（loudness）
- \(m_i\)：调制深度（物理量）
- \(s(f_m)\)：人耳敏感性（心理声学）
三者的乘积，就是人耳感受到的粗糙度。
希望这篇指南能帮你打开粗糙度计算的大门。如果你有具体的应用场景或代码问题，随时可以深入交流。