94

アップデート：

この質問に基づいた Scipy レシピを見つけました! したがって、興味のある方は、目次 » 信号処理 » バターワース バンドパスに直接アクセスしてください。

1-D numpy 配列 (時系列) にバターワース バンドパス フィルターを実装するという、当初は単純なタスクに見えたものを達成するのに苦労しています。

含める必要があるパラメーターは、sample_rate、ヘルツのカットオフ周波数、および場合によっては次数です (減衰、固有周波数などの他のパラメーターは私にはわかりにくいため、任意の「デフォルト」値で十分です)。

私が今持っているのはこれです。これはハイパスフィルターとして機能しているようですが、正しく行っているかどうかはわかりません:

def butter_highpass(interval, sampling_rate, cutoff, order=5):
    nyq = sampling_rate * 0.5

    stopfreq = float(cutoff)
    cornerfreq = 0.4 * stopfreq  # (?)

    ws = cornerfreq/nyq
    wp = stopfreq/nyq

    # for bandpass:
    # wp = [0.2, 0.5], ws = [0.1, 0.6]

    N, wn = scipy.signal.buttord(wp, ws, 3, 16)   # (?)

    # for hardcoded order:
    # N = order

    b, a = scipy.signal.butter(N, wn, btype='high')   # should 'high' be here for bandpass?
    sf = scipy.signal.lfilter(b, a, interval)
    return sf

ここに画像の説明を入力

ドキュメントと例は紛らわしくわかりにくいですが、「バンドパス用」とマークされたコマンドに示されているフォームを実装したいと思います。コメント内の疑問符は、何が起こっているのかを理解せずに例をコピーして貼り付けただけの場所を示しています。

私は電気工学者でも科学者でもありません。医療機器の設計者であり、EMG 信号に対してかなり単純なバンドパス フィルタリングを実行する必要があります。

4

3 に答える 3

129

buttord の使用をスキップして、代わりにフィルターの順序を選択して、フィルター基準を満たすかどうかを確認することができます。バンドパス フィルターのフィルター係数を生成するには、butter() にフィルター次数、カットオフ周波数Wn=[low, high](サンプリング周波数の半分であるナイキスト周波数の分数として表される)、およびバンド タイプを指定しbtype="band"ます。

以下は、バターワース バンドパス フィルターを操作するためのいくつかの便利な関数を定義するスクリプトです。スクリプトとして実行すると、2 つのプロットが作成されます。1 つは、同じサンプリング レートとカットオフ周波数に対する複数のフィルター次数での周波数応答を示しています。もう 1 つのプロットは、サンプル時系列に対するフィルター (次数 = 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()

このスクリプトによって生成されるプロットは次のとおりです。

複数のフィルター次数の周波数応答

ここに画像の説明を入力

于 2012-09-02T06:41:24.197 に答える
4

バンドパス フィルターの場合、 ws は下側と上側のコーナー周波数を含むタプルです。これらは、フィルター応答が通過帯域より 3 dB 小さいデジタル周波数を表します。

wp は、阻止帯域のデジタル周波数を含むタプルです。これらは、最大減衰が始まる位置を表します。

gpass は dB 単位の通過帯域の最大減衰であり、gstop は阻止帯域の減衰です。

たとえば、コーナー周波数が 300 Hz と 3100 Hz で、サンプリング レートが 8000 サンプル/秒のフィルターを設計したいとします。ナイキスト周波数は、サンプルレートを 2 で割ったもので、この例では 4000 Hz です。同等のデジタル周波数は 1.0 です。2 つのコーナー周波数は 300/4000 と 3100/4000 です。

ここで、阻止帯域をコーナー周波数から 30 dB +/- 100 Hz 下げたいとします。したがって、阻止帯域は 200 および 3200 Hz で始まり、200/4000 および 3200/4000 のデジタル周波数になります。

フィルターを作成するには、buttord を次のように呼び出します。

fs = 8000.0
fso2 = fs/2
N,wn = scipy.signal.buttord(ws=[300/fso2,3100/fso2], wp=[200/fs02,3200/fs02],
   gpass=0.0, gstop=30.0)

得られるフィルタの長さは、阻止帯域の深さと、コーナー周波数と阻止帯域周波数の差によって決まる応答曲線の急峻さに依存します。

于 2012-08-23T15:26:15.697 に答える