がれすたさんのDIY日記

電子回路、Python、組み込みシステム開発、自作エフェクターを語るblog

Python:scipy.signal.freqzのテスト

PythonMATLABみたいなfreqzが使えるみたいなのでためしてみた

RBJ cookbookの公式を用いてPEAKフィルタを構成した。

fc1:150Hz カット
fc2:400Hz ブースト

みたいな感じにした。
カットオフ帯域ゲインとbwあるけどパラメータはご自由に


wがそのままだと正規化周波数のまま帰ってくるので周波数に変えたければナイキスト周波数(サンプリング周波数/2.0*pi)をかけてあげればいい
オプションパラメータでfsを指定すればfsスケールで計算することも可能(Scipy ver1.2必須)

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np


def calc_coef_PF(fc,g,bw):
    a = [0.0] * 3
    b = [0.0] * 3

    omega = 2.0 * np.pi * fc / 44100 #正規化カットオフ
    alpha = np.sin(omega) * np.sinh(np.log(2.0) / 2.0 * bw * omega / np.sin(omega))
    A = np.power(10.0,(g/40.0))

    a[0] =  1.0 + alpha / A;
    a[1] = -2.0 * np.cos(omega);
    a[2] =  1.0 - alpha / A;
    b[0] =  1.0 + alpha * A;
    b[1] = -2.0 * np.cos(omega);
    b[2] =  1.0 - alpha * A;
    return a, b

a1, b1 = calc_coef_PF(150,-15,0.1)#[cutoff -> 150Hz][gain -> -15][bw -> 0.1]
a2, b2 = calc_coef_PF(400,2,1.0)  #[cutoff -> 400Hz][gain -> 2.0][bw -> 1.0]


w, h = signal.freqz(b1,a1,worN=1024*5)
w1, h1 = signal.freqz(b2,a2,worN=1024*5)

h = h*h1
f = w*(44100/2.0*np.pi)



fig, ax1 = plt.subplots()
plt.xticks([10,100,1000,10000], ["10","100","1k","10k"])
plt.xscale("log")
plt.grid(which="both")
ax1.set_title('Digital filter frequency response')
ax1.plot(f[0:100000], 20 * np.log10(abs(h)), 'b')
ax1.set_ylabel('Amplitude [dB]', color='b')
ax1.set_xlabel('Frequency [Hz]')

plt.show()

f:id:gsmcustomeffects:20190108043005p:plain
周波数応答


アナログではこういう回路組むけどパラメータ調整次第で同じような特性も出せそうなので頑張ろうと思う

f:id:gsmcustomeffects:20190108043242p:plain
Spiceでのシミュレーション結果

個別にしたやつ

畳み込み使わなくてもフィルタかける前に係数合わせてしまって高次フィルタにしちゃえばもっと短くも書けると思う(解釈が間違ってたらすいません

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np


def calc_coef_PF(fc,g,bw):
    a = [0.0] * 3
    b = [0.0] * 3

    omega = 2.0 * np.pi * fc / 44100 #正規化カットオフ
    alpha = np.sin(omega) * np.sinh(np.log(2.0) / 2.0 * bw * omega / np.sin(omega))
    A = np.power(10.0,(g/40.0))

    a[0] =  1.0 + alpha / A;
    a[1] = -2.0 * np.cos(omega);
    a[2] =  1.0 - alpha / A;
    b[0] =  1.0 + alpha * A;
    b[1] = -2.0 * np.cos(omega);
    b[2] =  1.0 - alpha * A;
    return a, b

a1, b1 = calc_coef_PF(150,-2,1.0)#[cutoff -> 150Hz][gain -> -15][bw -> 0.1]
a2, b2 = calc_coef_PF(400,2,1.0)  #[cutoff -> 400Hz][gain -> 2.0][bw -> 1.0]
a3, b3 = calc_coef_PF(1000,2,1.0)  #[cutoff -> 400Hz][gain -> 2.0][bw -> 1.0]

w1, h1 = signal.freqz(b1,a1,worN=1024*5)
w1, h2 = signal.freqz(b2,a2,worN=1024*5)
w1, h3 = signal.freqz(b3,a3,worN=1024*5)

h = h1*h2*h3
f = w*(44100/2.0*np.pi)



fig, (ax1,ax2) = plt.subplots(nrows=2,figsize=(10,15), sharex=True)


plt.title('Digital filter frequency response')
plt.xscale("log")

ax1.plot(f[0:100000], 20 * np.log10(abs(h)), 'b',label='conv all filter')
ax1.legend()
ax1.set_ylim([-40,40])
ax1.set_ylabel('Amplitude [dB]', color='b')
ax1.set_xlabel('Frequency [Hz]')
ax1.grid(which='minor')

ax2.plot(f[0:100000], 20 * np.log10(abs(h1)),label='filter1' )
ax2.plot(f[0:100000], 20 * np.log10(abs(h2)),label='filter2' )
ax2.plot(f[0:100000], 20 * np.log10(abs(h3)),label='filter3' )
ax2.legend()
ax2.set_ylim([-40,40])
ax2.set_ylabel('Amplitude [dB]', color='b')
ax2.set_xlabel('Frequency [Hz]')
ax2.grid(which='minor')
plt.show()

f:id:gsmcustomeffects:20190108065834p:plain
frequency response