まず、ITU-T G.227 https://www.itu.int/rec/T-REC-G.227-198811-I/en
記載の伝達関数を周波数特性としてグラフ化。普通のフィルタのグラフと違ってyはゲインではなくロスなので注意
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)
dB = 20 * np.log10(loss)
# adjust gain
dB -= np.min(dB)
loss = np.power(10, dB / 20)
return f, loss
sr = 44100
firlen = 513
xticks = [33, 50, 100, 500, 1000, 5000, 10000]
f_m, loss_m = g227(sr/2, 2**13)
f, loss = g227(sr/2, int(firlen/2))
print(np.min(20 * np.log10(loss)))
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, 20 * np.log10(loss), label="length=%d" % int(firlen/2))
plt.plot(f_m, 20 * np.log10(loss_m), label="expected")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
plt.show()
0.0
実数フィルタにしたいので左右対称にする
gain = 1.0/loss
gain = list(gain) + list(gain[::-1])
print(len(gain))
plt.figure(figsize=(10,5))
plt.title("Frequency response")
plt.grid(True)
plt.ylabel('Composite loss')
plt.xlabel('Frequency')
plt.plot(gain)
512
[<matplotlib.lines.Line2D at 0x12bd07850>]
IFFT して有限インパルス応答 (FIR) を得る。
fir = np.fft.fftshift(np.fft.ifft(gain))
print(len(fir))
x = np.linspace(0, len(fir), len(fir))
plt.figure(figsize=(15,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.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
plt.show()
fir = np.real(fir)
512
求めた FIR の周波数特性を評価する。freqz を使う。
w, h = signal.freqz(fir, worN=2**13)
x = w * sr * 1.0 / (2 * np.pi)
plt.figure(figsize=(10,10))
plt.title("Frequency response of impulse response (complex)")
plt.xscale('log')
plt.grid(True)
plt.ylabel('Gain')
plt.xlabel('Frequency')
plt.xlim((1, sr/2))
plt.ylim((-100, 10))
plt.xticks(xticks, xticks)
freqz_dB = 20 * np.log10(abs(h))
expected_dB = -20 * np.log10(loss_m)
plt.plot(x, 20 * np.log10(np.abs(h)), 'b', label="freqz")
plt.plot(f_m, -20 * np.log10(loss_m), 'r', label="expected")
plt.legend(bbox_to_anchor=(1, 1), loc='upper right', borderaxespad=1, fontsize=11)
print(len(h) / np.sum(np.abs(h)))
error = np.abs(freqz_dB- expected_dB)
print('Max Error', np.max(error), 'dB', 'at', x[np.argmax(error)], 'Hz')
print('MSE', 20*np.log10(np.sqrt(np.mean(np.power(10, error / 20)**2))), 'dB')
ax2 = plt.twinx()
angles = np.unwrap(np.angle(h))
ax2.plot(x, angles, 'g')
ax2.set_ylabel('Angle (radians)', color='g')
plt.show()
12.674339983109126 Max Error 6.3715686374777505 dB at 43.06640625 Hz MSE 0.1451089475236209 dB
print(json.dumps(list(map(lambda x: x.real, fir))))
[0.0, -4.959127400394614e-09, -1.984185008657713e-08, -4.465219102106621e-08, -7.942086416995453e-08, -1.2414527018542176e-07, -1.7889755019734332e-07, -2.4365250187244336e-07, -3.1853981026855327e-07, -4.034956754073518e-07, -4.987235258865907e-07, -6.04104936302943e-07, -7.199333812793213e-07, -8.460197299968095e-07, -9.82763928783579e-07, -1.1298899900745799e-06, -1.2879209015173587e-06, -1.4564775768222127e-06, -1.6362227797123846e-06, -1.8266579731376908e-06, -2.0286026184684197e-06, -2.241422246173838e-06, -2.4661101459249463e-06, -2.701879284430775e-06, -2.949914142313416e-06, -3.209258321477601e-06, -3.4813051117878885e-06, -3.7649117590110234e-06, -4.061698260547121e-06, -4.370318303494741e-06, -4.692636796376631e-06, -5.027086432512519e-06, -5.3757955666084226e-06, -5.73695820863928e-06, -6.11298505329063e-06, -6.501813460757424e-06, -6.906155746604949e-06, -7.323674354882317e-06, -7.757402919627947e-06, -8.204710378781777e-06, -8.668971830077653e-06, -9.147243767364515e-06, -9.64326337730666e-06, -1.0153755398442741e-05, -1.0682840245601535e-05, -1.1226891191535096e-05, -1.1790433568125216e-05, -1.2369469045634502e-05, -1.2968950149116858e-05, -1.3584486355499139e-05, -1.4221480285874031e-05, -1.4875128149902713e-05, -1.5551306236107152e-05, -1.6244775899739317e-05, -1.6961911380806126e-05, -1.7697017048562327e-05, -1.8456990137882953e-05, -1.923565532345765e-05, -2.0040458687275477e-05, -2.0864721890036892e-05, -2.171646657448084e-05, -2.2588487421751867e-05, -2.3489409266228022e-05, -2.4411475160947054e-05, -2.536394173953318e-05, -2.6338475057037535e-05, -2.7344993193962946e-05, -2.837455907607055e-05, -2.9437782986066387e-05, -3.052509778586176e-05, -3.164783789562299e-05, -3.2795778331889025e-05, -3.398101084484085e-05, -3.5192623931551035e-05, -3.644350120491019e-05, -3.772201502819713e-05, -3.904187683880547e-05, -4.039071226188886e-05, -4.178309804581245e-05, -4.3205881431311285e-05, -4.4674543591735995e-05, -4.617512064092235e-05, -4.772403902962558e-05, -4.9306489849605965e-05, -5.093988753946679e-05, -5.260854306214408e-05, -5.433090354192379e-05, -5.609036343328343e-05, -5.7906449371470953e-05, -5.976160158633649e-05, -6.16764753286691e-05, -6.363251748493781e-05, -6.565156347032271e-05, -6.771402623828788e-05, -6.984297554103291e-05, -7.201774826802129e-05, -7.426270550056342e-05, -7.655606431931352e-05, -7.892353715999775e-05, -8.134217586166961e-05, -8.383910750683476e-05, -8.639017149674702e-05, -8.902397637658617e-05, -9.171510007392595e-05, -9.449370321785316e-05, -9.733305131033501e-05, -0.00010026493180128194, -0.00010326124482353178, -0.00010635548384306482, -0.00010951812861477537, -0.00011278446265340029, -0.0001161234881916131, -0.00011957236808368872, -0.0001230985676953846, -0.0001267412242372491, -0.00013046620460585286, -0.00013431472163271082, -0.00013825097983880016, -0.00014231837577316416, -0.0001464793853387014, -0.00015077970438569102, -0.00015518001160764344, -0.0001597284259645241, -0.0001643837580128397, -0.00016919668268836193, -0.0001741240691918132, -0.00017921929130233424, -0.00018443720144192387, -0.00018983402617491624, -0.00019536252365793521, -0.00020108193948054572, -0.00020694285819144672, -0.00021300772434456868, -0.00021922486797551152, -0.00022566012784981947, -0.00023225949742219535, -0.00023909242207963662, -0.0002461024759991719, -0.0002533629429060038, -0.0002608148950738784, -0.00026853570807728374, -0.0002764638706390709, -0.00028468112838272715, -0.0002931233069731303, -0.00030187682834833846, -0.0003108747792266623, -0.00032020859614325097, -0.0003298085564652063, -0.00033977148625734797, -0.0003500247909549513, -0.00036067110317874275, -0.0003716349045942219, -0.0003830250999040192, -0.0003947632095296834, -0.00040696493182896566, -0.0004195488073384954, -0.0004326379145897345, -0.0004461478199147814, -0.0004602096439744591, -0.0004747360156012123, -0.0004898668473424183, -0.0005055119064009683, -0.0005218207493234381, -0.0005387004065441233, -0.000556311050173728, -0.0005745571595068789, -0.0005936106332608724, -0.0006133736599924196, -0.0006340311389062989, -0.0006554833194958214, -0.000677929565279607, -0.0007012686488595508, -0.0007257160830749518, -0.0007511697584014124, -0.0007778632788530487, -0.0008056944050980714, -0.0008349170712754429, -0.0008654298456843472, -0.0008975095733483262, -0.0009310567822648585, -0.0009663741995441969, -0.001003365709758295, -0.0010423633350888492, -0.0010832759870345252, -0.001126468889455813, -0.0011718579482510798, -0.0012198460378607735, -0.0012703583363853118, -0.0013238403997695848, -0.0013802292611085646, -0.0014400187928258776, -0.0015031607350125637, -0.0015702035069642276, -0.0016411165936085035, -0.0017165097288590133, -0.0017963732108215885, -0.0018813852588352186, -0.0019715598214604555, -0.0020676509290944536, -0.00216969837151045, -0.0022785390559375467, -0.002394239522564602, -0.0025177257078925746, -0.002649089585759948, -0.0027893503700483447, -0.0029386205499808216, -0.003098013498562655, -0.0032676517302780254, -0.0034487381806823933, -0.0036413865440421945, -0.003846876244409835, -0.004065281069918991, -0.004297931183878357, -0.004544811778141889, -0.004807259531427878, -0.005085097415542724, -0.005379598011402318, -0.005690313599881178, -0.006018344991135052, -0.006362817190679591, -0.006724500305447768, -0.0071018698222003315, -0.007495136332607807, -0.007901815016327852, -0.008321234275601063, -0.008749520309226853, -0.009184672621169019, -0.00962084511552893, -0.010054100939682292, -0.010475839221681606, -0.010879376145346589, -0.011252323265788036, -0.011584191483222876, -0.011857468357561326, -0.012056515755441801, -0.012157012994755315, -0.012136532899918333, -0.011961904332108548, -0.011602032884269971, -0.011012572898958151, -0.01015185955068889, -0.008962036205705385, -0.0073894848664154515, -0.005361172014796007, -0.0028119141471641936, 0.00034593471974791226, 0.004184191601468877, 0.00879886184843542, 0.014252671384170117, 0.02063024802919053, 0.0279437068941373, 0.03621433897973543, 0.04529488283578011, 0.055007627157714865, 0.06471525969702493, 0.07365598535160074, 0.07889769954583957, 0.07365598535160074, 0.06471525969702493, 0.055007627157714865, 0.04529488283578011, 0.036214338979735444, 0.027943706894137305, 0.02063024802919053, 0.014252671384170117, 0.00879886184843542, 0.0041841916014688785, 0.00034593471974791226, -0.0028119141471641927, -0.0053611720147960056, -0.0073894848664154515, -0.008962036205705383, -0.01015185955068889, -0.01101257289895815, -0.01160203288426997, -0.011961904332108548, -0.01213653289991833, -0.012157012994755316, -0.012056515755441801, -0.011857468357561327, -0.011584191483222876, -0.011252323265788034, -0.010879376145346587, -0.010475839221681608, -0.010054100939682294, -0.00962084511552893, -0.009184672621169017, -0.008749520309226853, -0.008321234275601063, -0.007901815016327852, -0.007495136332607806, -0.007101869822200331, -0.00672450030544777, -0.00636281719067959, -0.006018344991135053, -0.005690313599881178, -0.005379598011402316, -0.005085097415542724, -0.004807259531427877, -0.004544811778141889, -0.004297931183878357, -0.004065281069918991, -0.0038468762444098365, -0.003641386544042195, -0.0034487381806823933, -0.0032676517302780263, -0.003098013498562656, -0.002938620549980821, -0.0027893503700483434, -0.002649089585759948, -0.0025177257078925746, -0.0023942395225646012, -0.0022785390559375476, -0.0021696983715104525, -0.0020676509290944545, -0.001971559821460454, -0.0018813852588352201, -0.001796373210821591, -0.0017165097288590154, -0.0016411165936085026, -0.0015702035069642276, -0.0015031607350125629, -0.0014400187928258776, -0.0013802292611085638, -0.001323840399769585, -0.001270358336385315, -0.0012198460378607743, -0.0011718579482510776, -0.001126468889455814, -0.001083275987034525, -0.0010423633350888485, -0.0010033657097582938, -0.0009663741995441979, -0.0009310567822648593, -0.0008975095733483275, -0.0008654298456843464, -0.0008349170712754427, -0.0008056944050980744, -0.0007778632788530491, -0.000751169758401413, -0.0007257160830749526, -0.0007012686488595511, -0.000677929565279607, -0.000655483319495821, -0.0006340311389062985, -0.0006133736599924213, -0.0005936106332608737, -0.0005745571595068802, -0.000556311050173728, -0.0005387004065441219, -0.0005218207493234377, -0.0005055119064009672, -0.0004898668473424181, -0.00047473601560121337, -0.0004602096439744602, -0.0004461478199147812, -0.00043263791458973504, -0.00041954880733849505, -0.00040696493182896555, -0.00039476320952968355, -0.00038302509990401936, -0.0003716349045942217, -0.0003606711031787416, -0.00035002479095494975, -0.00033977148625734645, -0.00032980855646520466, -0.00032020859614325086, -0.0003108747792266624, -0.000301876828348339, -0.00029312330697312883, -0.00028468112838272476, -0.0002764638706390714, -0.0002685357080772822, -0.00026081489507387813, -0.0002533629429060039, -0.00024610247599917237, -0.0002390924220796362, -0.00023225949742219378, -0.0002256601278498179, -0.0002192248679755088, -0.00021300772434456683, -0.00020694285819144662, -0.00020108193948054423, -0.0001953625236579314, -0.00018983402617491624, -0.00018443720144192005, -0.00017921929130233275, -0.00017412406919181656, -0.00016919668268836356, -0.00016438375801284047, -0.00015972842596452427, -0.0001551800116076436, -0.00015077970438569062, -0.00014647938533870016, -0.00014231837577316492, -0.00013825097983879948, -0.00013431472163271034, -0.00013046620460585204, -0.000126741224237251, -0.000123098567695384, -0.00011957236808368837, -0.00011612348819161754, -0.00011278446265340105, -0.00010951812861477288, -0.00010635548384306504, -0.00010326124482353281, -0.00010026493180128254, -9.733305131033398e-05, -9.449370321785332e-05, -9.171510007392958e-05, -8.902397637658606e-05, -8.639017149674583e-05, -8.383910750683616e-05, -8.134217586167113e-05, -7.89235371599971e-05, -7.655606431931374e-05, -7.426270550056364e-05, -7.20177482680189e-05, -6.984297554103334e-05, -6.77140262382895e-05, -6.565156347032315e-05, -6.363251748493776e-05, -6.167647532866774e-05, -5.9761601586336437e-05, -5.790644937147106e-05, -5.609036343328126e-05, -5.4330903541922545e-05, -5.2608543062143426e-05, -5.0939887539465975e-05, -4.930648984960613e-05, -4.7724039029625146e-05, -4.6175120640923274e-05, -4.467454359173648e-05, -4.320588143130852e-05, -4.178309804581348e-05, -4.0390712261886906e-05, -3.904187683880422e-05, -3.7722015028197995e-05, -3.644350120491106e-05, -3.519262393155255e-05, -3.398101084484085e-05, -3.279577833188794e-05, -3.1647837895624944e-05, -3.052509778586035e-05, -2.943778298606628e-05, -2.8374559076068165e-05, -2.734499319396403e-05, -2.6338475057038728e-05, -2.536394173953307e-05, -2.4411475160948356e-05, -2.3489409266227914e-05, -2.2588487421753168e-05, -2.171646657448225e-05, -2.0864721890038193e-05, -2.0040458687274176e-05, -1.9235655323459385e-05, -1.8456990137882953e-05, -1.7697017048562977e-05, -1.6961911380806776e-05, -1.62447758997404e-05, -1.555130623610672e-05, -1.4875128149901629e-05, -1.422148028587338e-05, -1.3584486355497838e-05, -1.2968950149117075e-05, -1.236946904563602e-05, -1.1790433568125433e-05, -1.1226891191536397e-05, -1.0682840245601968e-05, -1.0153755398442741e-05, -9.64326337730536e-06, -9.147243767363648e-06, -8.668971830077653e-06, -8.204710378782645e-06, -7.757402919626212e-06, -7.323674354881883e-06, -6.906155746604949e-06, -6.501813460758725e-06, -6.1129850532888955e-06, -5.736958208636678e-06, -5.3757955666084226e-06, -5.027086432515121e-06, -4.692636796378366e-06, -4.370318303495609e-06, -4.0616982605479884e-06, -3.764911759010156e-06, -3.481305111789623e-06, -3.2092583214784684e-06, -2.949914142313416e-06, -2.70187928442904e-06, -2.4661101459249463e-06, -2.2414222461747055e-06, -2.0286026184675524e-06, -1.8266579731376908e-06, -1.6362227797115173e-06, -1.4564775768239474e-06, -1.2879209015173587e-06, -1.1298899900719778e-06, -9.827639287853136e-07, -8.460197299963758e-07, -7.199333812791045e-07, -6.041049363040272e-07, -4.987235258878918e-07, -4.034956754082192e-07, -3.1853981026855327e-07, -2.436525018741781e-07, -1.7889755019387388e-07, -1.2414527019236066e-07, -7.942086416995453e-08, -4.465219102106621e-08, -1.984185008657713e-08, -4.959127400394614e-09]