pythonで自作畳み込み演算とバンドストップフィルタ(BSF)を設計！

この記事では、

自力で1次元の畳み込み演算を実装して、バンドストップフィルタ(Band Stop Filter BSFもしくはBand Eliminate Filter BEF)を実装します!

Contents

畳み込み演算の実装
- 畳み込みの定義
ディジタルフィルタの設計
実際に関数を適応
最後に

畳み込み演算の実装

畳み込みの定義

入力信号と、システムの因果性を考慮して以下の式で実装します。

$y(n)= \sum^{n-1}_{k=0}{x(k)h(n-k)}$

これ元に参考になるべく繰り返しの回数が減るようにすると

def conv1d(x, h):
    """
    convolve x and h
    Parameters
    ----------
    x : array_like
        First one-dimensional input array.
    h : array_like
        Second one-dimensional input array.
    Returns
    -------
    out : ndarray
        Discrete, linear convolution of a and v.
    """
    N = len(x)
    out = np.zeros(N + len(h))
    for n in range(len(h)):
        out[n: n + N] += x * h[n]
    return out[:N]

こんな感じになります。

ディジタルフィルタの設計

バンドストップフィルタの設計

基本的にどんなフィルタでも考え方は同じです。

今回はバンドストップフィルタで行います。

Band Stop Filter を設計するためには周波数領域で、以下のようになるインパルス応答を実装すればよいです。（ただし$ω$は正規化された角周波数）

先の式を離散逆フーリエ変換した結果は以下のようになります。（ただし、 nは離散信号の1点を表す）

これがバンドストップフィルタのインパルス応答です。

なので、必要な分の窓幅$N$を定義して、そのサイズのバンドストップフィルタを返す関数を作成していきます。

バンドストップフィルタのインパルス応答

自作sinc関数を定義

numpyのsinc関数では定義が

$\frac{1}{πx}sin(πx)$となっており、式変形がわずらわしいので、自分で定義します。

ただ、自分で定義したことにより、入力が要素しか受け付けなくなっている点に注意が必要です。

def sinc(x):
    """
    sinc function
    Parameters
    ----------
    x : float
        ．input
    Returns
    -------
    y : float
        apply sinc function ．
    """
    if x == 0.0:
        return 1.0
    else:
        return np.sin(x) / x

入力のインデックスに対して、返答するインパルス応答を定義

from math import pi

def BSF_impulse(index, f_low, f_high, sr):
    """
    calculate BSF impulse at each point
    Parameters
    ----------
    idx : int
        index number
    f_low : float
        low edge frequency
    f_high : float
        high edge frequency
    sr : int
        sampling rate
    Returns
    -------
    y : float
        BSF impulse
    """
    # normalize sampling rate
    f_low = 2 * pi * f_low / sr
    f_high = 2 * pi * f_high / sr
    return sinc(pi * index) + (f_low * sinc(f_low * index) - f_high * sinc(f_high * index)) / pi

バンドストップフィルタを作成

今回は、窓関数法で適応をします。

フィルタは、0を中心に奇関数として定義をされているので、そうなるようにインデックスを入力します。

しかし、実際に使う際は、$(0, N+1)$の範囲で定義されていた方が使いやすいので、

b.append()のように逐次追加をしていき、最後にnumpy配列に変換をしています。

from scipy import hamming

def BandStopFilter(f_low, f_high, sr, N=30):
    """
    make band stop filter
    Parameters
    ----------
    f_low : float
        low edge frequency
    f_high : float
        high edge frequency
    N : int, default 30
        the number of filter coefficients
    sr : float
        sampling rate
    Returns
    -------
    b : ndarray
        Filter coefficient
    """
    if N % 2 != 0:
        N += 1
    N = int(N)
    b = []
    # calculate filter coefficients
    for idx in range(-N // 2, N // 2 + 1):
        b.append(BSF_impulse(idx, f_low, f_high, sr))
    # Apply window function method
    b = np.array(b)
    window_func = hamming(N + 1)
    return b * window_func

上で定義したBFS_impulse()関数を変更すれば任意のフィルタを設計することができます。

今回作成したフィルタは、

窓幅：$N = 100$

阻止域周波数：$f_{low} = 8000, f_{high} = 12000$

として定義をしました。

インパルス応答を確認

こんなフィルタが作成できます。

位相特性を確認

振幅特性を確認

どちらも定義した阻止域で周波数が遮断されていることが分かります。

実際に関数を適応

下の関数では、バンドストップフィルタ適応前後の周波数特性と位相特性と、バンドストップフィルタ適応前後のスペクトログラムを表示しています。

import matplotlib.pyplot as plt
import librosa

def main():    
    N = args.N  # [step]
    f_low = args.f_low  # [Hz]
    f_high = args.f_high  # [Hz]
    wav, sr = librosa.load(fname, sr=48000)  # [sec], [Hz]
    x = wav
    F_x = np.fft.fft(x * hamming(len(x)))
    amp = 20 * np.log10(np.abs(F_x))
    arg = np.unwrap(np.angle(F_x)) * 180 / pi

    frequency_range = np.arange(0, sr // 2, (sr // 2) / (wav.size // 2 + 1))
    plt.figure(figsize=(15, 5))
    plt.subplot(1, 2, 1)
    plt.title('Amplitude characteristic', size=16)
    plt.plot(frequency_range, amp[: len(amp) // 2 + 1])
    plt.xlabel("Frequency [Hz]", size=14)
    plt.ylabel("Amplitude [dB]", size=14)

    plt.subplot(1, 2, 2)
    plt.title('Phase characteristic', size=16)
    plt.plot(frequency_range, arg[: len(arg) // 2 + 1])
    plt.xlabel("Frequency [Hz]", size=14)
    plt.ylabel("Phase [rad]", size=14)
    plt.savefig("Original_Amp_Phase_char.png")
    plt.close()

    b = BandStopFilter(f_low, f_high, sr, N=N)
    x_hat = conv1d(x, b)

    F_x = np.fft.fft(x_hat * hamming(len(x_hat)))
    amp = 20 * np.log10(np.abs(F_x))
    arg = np.unwrap(np.angle(F_x)) * 180 / pi
    frequency_range = np.arange(0, sr // 2, (sr // 2) / (wav.size // 2 + 1))
    plt.figure(figsize=(15, 5))
    plt.subplot(1, 2, 1)
    plt.title('Amplitude characteristic', size=16)
    plt.plot(frequency_range, amp[: len(amp) // 2 + 1])
    plt.xlabel("Frequency [Hz]", size=14)
    plt.ylabel("Amplitude [dB]", size=14)

    plt.subplot(1, 2, 2)
    plt.title('Phase characteristic', size=16)
    plt.plot(frequency_range, arg[: len(arg) // 2 + 1])
    plt.xlabel("Frequency [Hz]", size=14)
    plt.ylabel("Phase [rad]", size=14)
    plt.savefig("BST_Amp_Phase_char.png")
    plt.close()

    plt.figure(figsize=(12, 5))
    plt.title("Original signal spectrogram", size=16)
    plt.xlabel("Time [sec]", size=14)
    plt.ylabel("Frequency [Hz]", size=14)
    spectrum, freq, _, _ = plt.specgram(x, Fs=sr)
    cbar = plt.colorbar()
    cbar.set_label("magnitude [dB]", size=14)
    plt.savefig("Original_signal_spectrogram.png")
    plt.close()

    plt.figure(figsize=(12, 5))
    plt.title("Band Stop Filter signal spectrogram", size=16)
    plt.xlabel("Time [sec]", size=14)
    plt.ylabel("Frequency [Hz]", size=14)
    _=plt.specgram(x_hat, Fs=sr)
    cbar = plt.colorbar()
    cbar.set_label("magnitude [dB]", size=14)
    plt.savefig("BSF_signal_spectrogram.png")
    plt.close()

    librosa.output.write_wav('BSF_signal.wav', x_hat, sr)
    print("over")


if __name__ == "__main__":
    main()