ある音に対して、人間が感じるうるささと音圧レベルの大きさは必ずしも一致していません。
以下の等ラウドネス曲線のように、音圧レベルと周波数によって人間の聴覚の感度は異なっています。
https://ja.wikipedia.org/wiki/音の大きさ
そのため、音信号の分析、特に騒音の分析では、A特性、B特性、C特性といった周波数重み付けがよく利用されます。
すなわち、人間の聴覚の感度に合うように、対象となる音信号の音圧レベルを補正して分析を行います。
Wikipedia の A-weighting の記事では、以下の4つの特性とその算出方法が記載されています。
- A特性: 等ラウドネス曲線*1 の40phonの逆特性
- B特性: 〃 の70phonの逆特性
- C特性: 〃 の100phonの逆特性
- D特性:
- 等ノイジネス曲線 (ISO. R 507, 音のやかましさ) に基づく
- 高レベルの航空機騒音測定向け
騒音のうるささの実用的な測定・分析では、騒音源の音圧レベルや種類、騒音を聞く環境を考慮することが求められます。従って、A特性だけでなく様々な特性が開発されてきました。
今回は、Wikipedia の A-weighting の記事 に基づいて、4つの周波数重み付け特性を実装・プロットしてました。（記事冒頭プロット画像を出力）
サンプルコード
注意点として、下記コードのX_f_splの単位は音圧レベル[dB SPL]、つまり最小可聴値を基準にした絶対数値です。
一方で、オーディオファイルの波形信号の振幅値は基本的に相対値です。
従って、測定時の基準となる音圧レベル（実効値など）を騒音計で記録あるいはマイク感度の校正を行い、それに基づいて測定／録音データ（オーディオファイル）の時間波形の振幅値を補正する必要があります。
import numpy as np
import matplotlib.pyplot as plt
def get_weight(freqs, kind):
if kind == "A":
func = lambda f: (12194**2*f**4) \
/(np.sqrt((f**2+107.7**2)*(f**2+737.9**2))\
*(f**2+20.6**2)*(f**2+12194**2))
bias = 20*np.log10(func(1000))
elif kind == "B":
func = lambda f: (12194**2*f**3) \
/(np.sqrt((f**2+158.8**2))\
*(f**2+20.6**2)*(f**2+12194**2))
bias = 20*np.log10(func(1000))
elif kind == "C":
func = lambda f: (12194**2*f**2) \
/((f**2+20.6**2)*(f**2+12194**2))
bias = 20*np.log10(func(1000))
elif kind == "D":
func_sub = lambda f: ((1037918.48 - f**2)**2 + 1080768.16*f**2) \
/ ((9837328 - f**2)**2 + 11723776 *f**2)
func = lambda f: (f/6.8966888496476) \
* np.sqrt( func_sub(f) / ((f**2+79919.29)*(f**2+1345600)) )\
/ 10**(-5)
bias = 0.0
else:
func = lambda f: 0.0*f + 1.0
bias = 0.0
return 20*np.log10(func(freqs)) - bias
def add_weight(freqs, power_dB, kind):
weight = get_weight(freqs, kind)
return power_dB + weight
freqs = 20 ** np.arange(1.0, 4.0+1e-9, 0.01)
X_f_spl = np.ones(freqs.shape)
sig_power_dB = 20 * np.log10(X_f_spl)
plt.clf()
for kind in ["A", "B", "C", "D", "Without"]:
power_dB_weighted = add_weight(freqs, sig_power_dB, kind)
plt.plot(freqs, power_dB_weighted, label=kind+"-weighting")
plt.xlim(20, 20000)
plt.xscale("log")
plt.xlabel("Frequency")
plt.ylabel("Gain dB")
plt.grid()
plt.legend()
plt.show()
参考