ITU-T G.227 https://www.itu.int/rec/T-REC-G.227-198811-I/en
import numpy as np
from matplotlib import pyplot as plt
from scipy import signal
import json
def g227(niq, size):
f = np.linspace(0, niq, size)
p = 1j * f / 1000
numerator = 18400 + 91238*p**2 + 11638*p**4 + p*(67280 + 54050*p**2)
denominator = 400 + 4001*p**2 + p**4 + p*(36040 + 130*p**2)
loss = np.abs(numerator / denominator)
return f, loss
sr = 44100
xticks = [33, 50, 100, 500, 1000, 5000, 10000]
f_m, loss_m = g227(sr/2, 2**15)
plt.figure(figsize=(8,10))
plt.title("Frequency response")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Composite loss')
plt.xlabel('Frequency')
plt.ylim((0, 70))
#plt.xlim((33, 10000))
plt.xlim((1, sr/2))
plt.xticks(xticks, xticks)
plt.plot(f_m, 20 * np.log10(loss_m), label="G.227")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
plt.show()
アナログフィルタの係数を計算する。以下は規格書に掲載されているアナログフィルタの式だが、p=f/1000 (Hz) になっているのを角周波数 ω をとる関数にする
$$ H(p) = \frac{ 11638\ p^4 + 54050\ p^3 + 91238\ p^2 + 67280\ p + 18400 }{ p^4 + 130\ p^3 + 4001\ p^2 + 36040\ p + 400 } $$$$ p = j \frac{\omega}{2000 \pi} $$以下の形にして係数を計算しなおす
$$ H(\omega) = \frac{ 11638\ (2000\pi)^{-4} ({j\omega})^4 + 54050\ (2000\pi)^{-3} ({j\omega})^3 + 91238\ (2000\pi)^{-2} ({j\omega})^2 + 67280\ (2000\pi)^{-1} ({j\omega}) + 18400 }{ (2000\pi)^{-4} ({j\omega})^4 + 130\ (2000\pi)^{-3} ({j\omega})^3 + 4001\ (2000\pi)^{-2} ({j\omega})^2 + 36040\ (2000\pi)^{-1} ({j\omega}) + 400 } $$biquad フィルタのみで構成する場合は、同じく G.227 に書いてあるネットワークを利用する
$$ \begin{eqnarray} H(p) = \frac{ 46 p^2 + 90 p + 46 }{ p^2 + 90 p + 1 } \\ H(p) = \frac{ 11 p + 20 }{ p + 20 } \\ H(p) = \frac{ 23 p + 20 }{ p + 20 } \\ \end{eqnarray} $$いずれにしてもこれはゲインの関数になっていて、フィルタとしては分母と分子が逆なので気をつける。
これを双一次変換してデジタルフィルタを得る filtz = signal.lti(*signal.bilinear(num, den, sr))
の部分
# biquad フィルタのみで構成する場合
network = [
[
[1, 90, 1],
[46, 90, 46]
],
[
[0, 1, 20],
[0, 11, 20]
],
[
[0, 1, 20],
[0, 23, 20]
]
]
# 4次のIIRフィルタで構成する場合
network = [
[
[1, 130, 4001, 36040, 400],
[11638, 54050, 91238, 67280, 18400]
]
]
filters = []
for (num, den) in network:
num = [ x * ( (2 * np.pi * 1000)**-(4-i) ) for i, x in enumerate(num) ]
den = [ x * ( (2 * np.pi * 1000)**-(4-i) ) for i, x in enumerate(den) ]
# 双一次変換
filtz = signal.lti(*signal.bilinear(num, den, sr))
wz, hz = signal.freqz(filtz.num, filtz.den, worN=2**15)
ws, hs = signal.freqs(num, den, worN=sr*wz)
filters.append({
"filter": filtz,
"wz": wz,
"hz": hz,
})
plt.figure(figsize=(8,10))
plt.title("Frequency response")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Gain')
plt.xlabel('Frequency')
#plt.ylim((0, 70))
plt.xlim((33, 10000))
plt.xlim((1, sr/2))
plt.xticks(xticks, xticks)
freqz_dB = 20 * np.log10(np.abs(hz))
expected_dB = -20 * np.log10(loss_m)
plt.plot(wz*sr/(2*np.pi), 20 * np.log10(np.abs(hz)), 'b', label="IIR")
plt.plot(wz*sr/(2*np.pi), 20 * np.log10(np.abs(hs)), label="Expected")
#plt.plot(f_m, -20 * np.log10(loss_m), label="expected")
error = np.abs(freqz_dB - expected_dB)
print('最大誤差とその周波数: Max Error', np.max(error), 'dB', 'at', wz[np.argmax(error)]*sr/(2*np.pi), 'Hz')
print('二乗平均平方根誤差: RMSE', 20*np.log10(np.sqrt(np.mean(np.power(10, error / 20)**2))), 'dB')
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
plt.show()
最大誤差とその周波数: Max Error 17.766762355103175 dB at 21666.439819335934 Hz 二乗平均平方根誤差: RMSE 10.434663310271596 dB
双一次変換の IIR フィルタには誤差(歪み)が残る。特にナイキスト周波数に近いほど大きくなる。これを FIR フィルタで補正する。 FIR フィルタは高域の補正は短い長さですむ。
まず補正する周波数特性を算出する
h = np.multiply.reduce([ np.abs(f["hz"]) for f in filters])
x = filters[0]["wz"]*sr/(2*np.pi)
diff = -20 * np.log10(loss_m) - 20 * np.log10(np.abs(h))
diff = np.power(10, diff / 20)
plt.figure(figsize=(8,10))
plt.title("Frequency response")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Gain')
plt.xlabel('Frequency')
#plt.ylim((-100, 0))
#plt.xlim((33, 10000))
plt.xlim((1, sr/2))
plt.xticks(xticks, xticks)
plt.plot(x, 20 * np.log10(np.abs(h)), 'b', label="IIR")
plt.plot(f_m, -20 * np.log10(loss_m), label="Expected")
plt.plot(f_m, 20 * np.log10(diff), label="Error to be compensated")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
plt.show()
実数フィルタにするため、左右対称にして IFFT してインパルスレスポンスを得る
gain = list(diff) + list(diff[::-1])
plt.figure(figsize=(10,5))
plt.title("IFFT target frequency response")
plt.grid(True)
plt.ylabel('Gain')
plt.xlabel('Frequency')
plt.plot( 20 * np.log10(gain))
plt.xticks([0, len(gain)/2, len(gain)-1], [0, "+%d/-%d" % (sr/2, sr/2), 0])
plt.show()
fir = np.fft.ifft(gain)
fir = list(fir[:65])
fir = list(fir[::-1]) + list(fir[1:])
print(len(fir))
x = np.linspace(0, len(fir)-1, len(fir))
plt.figure(figsize=(10,5))
plt.title("Impulse response")
plt.grid(True)
plt.ylabel("Coeffs")
plt.xlabel("Samples")
plt.plot(x, np.real(fir), label="real")
plt.plot(x, np.imag(fir), label="imag")
plt.xticks([64], [64])
#plt.xlim((60,68))
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
plt.show()
fir = np.real(fir)
129
IIR と補正 FIR を組合せて、誤差を評価する。
w, h = signal.freqz(fir, worN=2**15)
x = w * sr * 1.0 / (2 * np.pi)
plt.figure(figsize=(10,10))
plt.title("Frequency response (IIR+FIR)")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Gain')
plt.xlabel('Frequency')
plt.xlim((1, sr/2))
plt.ylim((-70, 10))
plt.xticks(xticks, xticks)
h = np.multiply.reduce([ np.abs(f["hz"]) for f in filters]) * np.abs(h)
freqz_dB = 20 * np.log10(h)
expected_dB = -20 * np.log10(loss_m)
plt.plot(x, freqz_dB, 'b', label="IIR+FIR")
plt.plot(f_m, expected_dB, 'r', label="expected")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
error = np.abs(freqz_dB- expected_dB)
print('最大誤差とその周波数: Max Error', np.max(error), 'dB', 'at', wz[np.argmax(error)]*sr/(2*np.pi), 'Hz')
print('二乗平均平方根誤差: RMSE', 20*np.log10(np.sqrt(np.mean(np.power(10, error / 20)**2))), 'dB')
最大誤差とその周波数: Max Error 0.04974724191783508 dB at 22049.327087402344 Hz 二乗平均平方根誤差: RMSE 0.0019852001730218663 dB
# 係数を出力
data = [
{
"num": list(f["filter"].num),
"den": list(f["filter"].den),
}
for f in filters
]
data.append(list(fir))
print(json.dumps(data, indent=2))
[ { "num": [ 0.0027293620390779223, 0.0002239451280681405, -0.002162405616471103, -0.0007195307151076326, -6.105389511747936e-05 ], "den": [ 1.0, -3.3926835929532446, 4.312959033230201, -2.4347384558596867, 0.5149375948434236 ] }, [ -0.00034784079278272323, 0.00035910574883141135, -0.0003708374974680221, 0.00038327793286990866, -0.00039619927956737475, 0.0004099876593157417, -0.00042426286213959426, 0.0004396050209241575, -0.00045542725478700096, 0.0004725704195613511, -0.0004901681025685535, 0.0005094111811268232, -0.0005290559939055743, 0.0005507628426983046, -0.0005727799686419476, 0.0005973965867091121, -0.0006221778700522442, 0.0006502547671927741, -0.0006782756914777671, 0.000710497051269187, -0.0007423386638185929, 0.0007795603458064955, -0.0008159374296345047, 0.0008592362492915217, -0.0009010330034797562, 0.0009517698709320109, -0.0010000837043960351, 0.0010599825539323803, -0.0011161744286636785, 0.0011874162003435063, -0.0012531603211528694, 0.0013384838490677756, -0.0014157961845169013, 0.0015185798228837788, -0.0016097750047847116, 0.0017340315561666319, -0.001841492121691065, 0.001991618072479125, -0.002117118266651802, 0.002297038366106365, -0.002440059903865149, 0.0026509738593397487, -0.00280479452656793, 0.0030397876383340817, -0.0031826987966766585, 0.003414403583244268, -0.0034902765501177574, 0.003643140784362165, -0.0035185056606318052, 0.003406610076082206, -0.002775004443246457, 0.001961087743010665, -0.00012539059974168645, -0.0024071680300668166, 0.007048899286148081, -0.013696453162376332, 0.024947163573301463, -0.041616742389665465, 0.06918097408309692, -0.11210433316434221, 0.18560341816508352, -0.31478694745543906, 0.5796783480514935, -1.2164627712267784, 2.66715630727313, -1.2164627712267784, 0.5796783480514935, -0.31478694745543906, 0.18560341816508352, -0.11210433316434221, 0.06918097408309692, -0.041616742389665465, 0.024947163573301463, -0.013696453162376332, 0.007048899286148081, -0.0024071680300668166, -0.00012539059974168645, 0.001961087743010665, -0.002775004443246457, 0.003406610076082206, -0.0035185056606318052, 0.003643140784362165, -0.0034902765501177574, 0.003414403583244268, -0.0031826987966766585, 0.0030397876383340817, -0.00280479452656793, 0.0026509738593397487, -0.002440059903865149, 0.002297038366106365, -0.002117118266651802, 0.001991618072479125, -0.001841492121691065, 0.0017340315561666319, -0.0016097750047847116, 0.0015185798228837788, -0.0014157961845169013, 0.0013384838490677756, -0.0012531603211528694, 0.0011874162003435063, -0.0011161744286636785, 0.0010599825539323803, -0.0010000837043960351, 0.0009517698709320109, -0.0009010330034797562, 0.0008592362492915217, -0.0008159374296345047, 0.0007795603458064955, -0.0007423386638185929, 0.000710497051269187, -0.0006782756914777671, 0.0006502547671927741, -0.0006221778700522442, 0.0005973965867091121, -0.0005727799686419476, 0.0005507628426983046, -0.0005290559939055743, 0.0005094111811268232, -0.0004901681025685535, 0.0004725704195613511, -0.00045542725478700096, 0.0004396050209241575, -0.00042426286213959426, 0.0004099876593157417, -0.00039619927956737475, 0.00038327793286990866, -0.0003708374974680221, 0.00035910574883141135, -0.00034784079278272323 ] ]