Как реализовать полосовой фильтр Баттерворта с помощью Scipy.сигнал.масло

Question

Как реализовать полосовой фильтр Баттерворта с помощью Scipy.сигнал.масло

обновление:

Я нашел рецепт Scipy, основанный на этом вопросе! Итак, для всех, кто заинтересован, перейдите прямо к:

http://scipy-cookbook.readthedocs.io/items/ButterworthBandpass.html

Мне трудно достичь того, что изначально казалось простой задачей реализации полосового фильтра Баттерворта для 1-D массива numpy (time-series).

параметры, которые я должен включить являются 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

enter image description here

документы и примеры сбивают с толку и неясны, но я хотел бы реализовать форму, представленную в похвале, помеченной как "для бандпаса". Вопросительные знаки в комментариях покажите, где я просто скопировал-вставил какой-то пример, не понимая, что происходит.

Я не электротехника или ученый, просто разработчик медицинского оборудования, которому нужно выполнить довольно простую полосовую фильтрацию на ЭМГ-сигналах.

Спасибо за любую помощь!

1437 3

python scipy signal-processing digital-filter

3 ответов:

Comments

Ничего не найдено.

Warren Weckesser · Accepted Answer · 2012-09-02 09:41:24

вы можете пропустить использование 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()
вот графики, которые генерируются этим скриптом:

user13107 · Accepted Answer · 2018-03-27 03:23:41

метод проектирования фильтра в принято отвечать правильно, но у него есть недостаток. SciPy полосовые фильтры, разработанные с b, A являются нестабильная и может привести к ошибочное фильтров at более высокие заказы фильтра.

вместо этого используйте выход sos (секций второго порядка) конструкции фильтра.
from scipy.signal import butter, sosfilt, sosfreqz

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

def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
        sos = butter_bandpass(lowcut, highcut, fs, order=order)
        y = sosfilt(sos, data)
        return y
кроме того, вы можете построить частотную характеристику, изменив
b, a = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = freqz(b, a, worN=2000)
до
sos = butter_bandpass(lowcut, highcut, fs, order=order)
w, h = sosfreqz(sos, worN=2000)

sizzzzlerz · Accepted Answer · 2012-08-23 18:26:15

для полосового фильтра, WS-это кортеж, содержащий нижнюю и верхнюю граничные частоты. Они представляют цифровую частоту, где отклик фильтра на 3 дБ меньше полосы пропускания.

wp-это кортеж, содержащий стоп-полосу цифровых частот. Они представляют собой место, где начинается максимальное затухание.

gpass-это максимальное затухание в полосе пропускания в дБ, в то время как gstop-это аттенюация в полосах останова.

скажем, например, вы требуется разработать фильтр для частоты дискретизации 8000 сэмплов / сек, имеющий угловые частоты 300 и 3100 Гц. Частота Найквиста-это частота дискретизации, деленная на два, или в данном примере 4000 Гц. Эквивалентная цифровая частота равна 1,0. Затем эти две угловые частоты составляют 300/4000 и 3100/4000.

теперь предположим, что вы хотели, чтобы стоп - полосы были ниже 30 дБ + / - 100 Гц от угловых частот. Таким образом, ваши стоп-полосы будут начинаться с 200 и 3200 Гц, что приведет к цифровому частоты 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)
длина результирующего фильтра будет зависеть от глубины стоп-полос и крутизны кривой отклика, которая определяется разностью между угловой частотой и частотой стоп-полос.