プログラミング

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

この記事では、

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

畳み込み演算の実装

畳み込みの定義

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

$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()

 

元の波の振幅特性と位相特性

BSF適応後の振幅特性と位相特性

位相特性はわかりにくいですが、振幅特性は削られていることが分かります。

元の波のスペクトログラム

BSF後のスペクトログラム

阻止域周波数:$f_{low} = 8000, f_{high} = 12000$の部分がごっそり取り除かれていることが分かります。

目的のフィルタが完成したと考えて良いのではないでしょうか!

最後に

今回はバンドストップフィルタでしたが、同じ方法でどんなフィルタでも作成できることが分かりました。

手元でフィルタを作成しているとCNNがどういった働きをしているのかがよく分かります。

すごく荒く説明をすると、目的の周波数を強調するフィルタを作成して、その周波数成分をもとに回帰問題を解いているのです。

何事も基礎から学んでいくことは大切だとしみじみ感じました!

pythonで自作短時間フーリエ変換を用いたスペクトログラム関数作成!こちらの記事で作成をした短時間フーリエ変換の関数を用いて、スペクトログラム関数を自作していきます! https://ibuki-s...
pythonのnumpyで自作のSTFT関数とISTFT関数を実装!今年の4月に某大学のデジタル信号を扱う研究室に配属されました! 希望が通り一安心です。 そんな思いで研究をしています...
オススメのプログラミングスクールをご紹介

タイピングもままならない完全にプログラミング初心者から

アホいぶきんぐ
アホいぶきんぐ
プログラミングってどこの国の言語なの~?

たった二ヶ月で

いぶきんぐ
いぶきんぐ
え!?人工知能めっちゃ簡単にできるじゃん!

応用も簡単にできる…!!

という状態になるまで、一気に成長させてくれたオススメのプログラミングスクールをご紹介します!

テックアカデミーのPython+AIコースを受講した僕が本音のレビュー・割引あり! というプログラミング完全初心者だった僕が Tech Academy(テックアカデミー)のPython×AIコース を二ヶ月間...

COMMENT

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です