[python] Scipy.signal.butter로 대역 통과 버터 워스 필터를 구현하는 방법

buttord의 사용을 건너 뛰고 대신 필터 주문을 선택하고 필터링 기준을 충족하는지 확인할 수 있습니다. 대역 통과 필터에 대한 필터 계수를 생성하려면 butter ()에 필터 차수, 차단 주파수 Wn=[low, high](샘플링 주파수의 절반 인 나이 퀴 스트 주파수의 일부로 표현됨) 및 대역 유형을 제공하십시오 btype="band".

다음은 Butterworth 대역 통과 필터로 작업하기위한 몇 가지 편의 함수를 정의하는 스크립트입니다. 스크립트로 실행하면 두 개의 플롯이 생성됩니다. 하나는 동일한 샘플링 속도 및 차단 주파수에 대해 여러 필터 차수의 주파수 응답을 보여줍니다. 다른 그림은 샘플 시계열에 대한 필터 (차수 = 6)의 효과를 보여줍니다.

from scipy.signal import butter, lfilter


def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y


if __name__ == "__main__":
    import numpy as np
    import matplotlib.pyplot as plt
    from scipy.signal import freqz

    # Sample rate and desired cutoff frequencies (in Hz).
    fs = 5000.0
    lowcut = 500.0
    highcut = 1250.0

    # Plot the frequency response for a few different orders.
    plt.figure(1)
    plt.clf()
    for order in [3, 6, 9]:
        b, a = butter_bandpass(lowcut, highcut, fs, order=order)
        w, h = freqz(b, a, worN=2000)
        plt.plot((fs * 0.5 / np.pi) * w, abs(h), label="order = %d" % order)

    plt.plot([0, 0.5 * fs], [np.sqrt(0.5), np.sqrt(0.5)],
             '--', label='sqrt(0.5)')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain')
    plt.grid(True)
    plt.legend(loc='best')

    # Filter a noisy signal.
    T = 0.05
    nsamples = T * fs
    t = np.linspace(0, T, nsamples, endpoint=False)
    a = 0.02
    f0 = 600.0
    x = 0.1 * np.sin(2 * np.pi * 1.2 * np.sqrt(t))
    x += 0.01 * np.cos(2 * np.pi * 312 * t + 0.1)
    x += a * np.cos(2 * np.pi * f0 * t + .11)
    x += 0.03 * np.cos(2 * np.pi * 2000 * t)
    plt.figure(2)
    plt.clf()
    plt.plot(t, x, label='Noisy signal')

    y = butter_bandpass_filter(x, lowcut, highcut, fs, order=6)
    plt.plot(t, y, label='Filtered signal (%g Hz)' % f0)
    plt.xlabel('time (seconds)')
    plt.hlines([-a, a], 0, T, linestyles='--')
    plt.grid(True)
    plt.axis('tight')
    plt.legend(loc='upper left')

    plt.show()

다음은이 스크립트에 의해 생성 된 플롯입니다.