ITU-T G.227 の周波数特性

ITU-T G.227 https://www.itu.int/rec/T-REC-G.227-198811-I/en

In [1]:
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)) の部分

In [2]:
# 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

FIR フィルタで補正する

双一次変換の IIR フィルタには誤差(歪み)が残る。特にナイキスト周波数に近いほど大きくなる。これを FIR フィルタで補正する。 FIR フィルタは高域の補正は短い長さですむ。

まず補正する周波数特性を算出する

In [3]:
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 してインパルスレスポンスを得る

In [4]:
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 の合成特性を評価する

IIR と補正 FIR を組合せて、誤差を評価する。

In [5]:
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
In [6]:
# 係数を出力

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
  ]
]
In [ ]: