Wizard Notes

Pythonを使った音楽信号分析の技術録、作曲活動に関する雑記

Python:FIR/IIRフィルタの直列/並列接続の周波数応答(伝達関数)プロット (scipy.signal.freqz)

f:id:Kurene:20210925151016p:plain

複雑な信号処理プラグインを開発する場合、複数のFIR/IIRフィルタを組み合わせることがよくあります。

そこで,FIR/IIRフィルタの直列/並列接続時の周波数応答をプロットするサンプルコードを作成しました。

回路全体の伝達関数を確認するだけでなく、フィルタ単体や部分的な回路の伝達関数も確認できるように、複数の伝達関数を入力として与えられるように設計しました。

また、カットオフ周波数を周波数応答とともに確認できるように垂直線プロットしています。

なお、直列回路全体の伝達関数各フィルタの伝達関数の積並列回路全体の伝達関数各フィルタの伝達関数の和として計算することができます。

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


def bode_plot(
        w, h, 
        sr=None,
        labels=None,
        ylim=[-25, 5], 
        to_dB = lambda h : 20 * np.log10(np.abs(h)),
        xscale="linear",
        fc=[],
        ):
    
    h_list = [h] if type(h) is not list else h
    labels = [None] * len(h_list) if type(labels) is not list else labels
    fc_list = [fc] if type(fc) is not list else fc
    
    n_h = len(h_list)
    alpha = (1/n_h)**0.5
    
    freqs = w/np.pi*(sr/2) if sr is not None else w
    
    plt.figure(figsize=(6.0, 8.0))
    plt.subplot(2,1,1)
    for k in range(n_h):
        plt.plot(freqs, to_dB(h_list[k]), label=labels[k], alpha=alpha)
    plt.ylabel('Magnitude (dB)')
    plt.xlabel('Frequency')
    plt.xscale(xscale)
    plt.grid()
    plt.ylim(ylim[0], ylim[1])
    for tmp_fc in fc_list:
        if tmp_fc < freqs[-1]:
            plt.vlines(tmp_fc, ylim[0], ylim[1], linestyles="dotted", colors ="k", alpha=0.5)

    if len(labels) > 0:
        plt.legend(loc='lower right')

    
    plt.subplot(2,1,2)
    for k in range(n_h):
        plt.plot(freqs, np.angle(h_list[k]) * 180 / np.pi, label=labels[k], alpha=alpha)
    plt.ylabel('Phase (deg)')
    plt.xlabel('Frequency')
    plt.xscale(xscale)
    plt.grid()
    plt.ylim(-180, 180)
    for tmp_fc in fc_list:
        if tmp_fc < freqs[-1]:
            plt.vlines(tmp_fc, -180, 180, linestyles="dotted", colors ="k", alpha=0.5)

    if len(labels) > 0:
        plt.legend(loc='lower right')


    plt.tight_layout()
    plt.show()


# パラメタ
sr = 44100
fc_iir = sr//10 # カットオフ周波数
fc_fir = sr//20 # カットオフ周波数
order_iir = 7   # フィルタ次数
order_fir = 9 # フィルタ次数

# フィルタ係数算出(バターワースフィルタ)
b, a         = scipy.signal.butter(order_iir, fc_iir/(sr/2), btype='low')
w, h_iir_low = scipy.signal.freqz(b, a)

# フィルタ係数算出(FIR)
b = scipy.signal.firwin(order_fir, fc_fir/(sr/2), pass_zero=False)
a = 1.0
w, h_fir_high = scipy.signal.freqz(b, a)

# プロット
bode_plot(w, 
    [
        h_iir_low, 
        h_fir_high,
        h_iir_low * h_fir_high, # 直列接続
        h_iir_low + h_fir_high, # 並列接続
    ], 
    labels=[
        "IIR-Low", 
        "FIR-High", 
        "Cascade(IIR-Low, FIR-High)", 
        "Parallel(IIR-Low, FIR-High)"
   ],
   sr=sr,
   fc=[fc_iir, fc_fir], 
)


np.set_printoptions(precision=3)

import code
console = code.InteractiveConsole(locals=locals())
console.interact()