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

84

ОБНОВИТЬ:

Я нашел Scipy рецепт, основанный на этом вопросе! Итак, всем, кто интересуется, сразу переходите к: Содержание »Обработка сигналов» Полоса пропускания Баттерворта


Мне трудно достичь того, что изначально казалось простой задачей реализации полосового фильтра Баттерворта для одномерного массива 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

введите описание изображения здесь

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

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

Heltonbiker
источник
Я пробовал что-то в dsp.stackexchange, но они слишком много (больше, чем я могу справиться) сосредоточены на концептуальных проблемах разработки, а не на использовании scipy-функций.
heltonbiker

Ответы:

119

Вы можете пропустить использование батторда и вместо этого просто выбрать порядок для фильтра и посмотреть, соответствует ли он вашему критерию фильтрации. Чтобы сгенерировать коэффициенты фильтра для полосового фильтра, задайте Butter () порядок фильтра, частоты среза Wn=[low, high](выраженные как часть частоты Найквиста, которая составляет половину частоты дискретизации) и тип полосы btype="band".

Вот сценарий, который определяет несколько удобных функций для работы с полосовым фильтром Баттерворта. При запуске в виде сценария он создает два графика. Один показывает частотную характеристику для нескольких порядков фильтрации для одной и той же частоты дискретизации и частот среза. Другой график демонстрирует влияние фильтра (с порядком = 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()

Вот графики, которые генерирует этот скрипт:

Частотная характеристика для нескольких порядков фильтров

введите описание изображения здесь

Уоррен Векессер
источник
1
Вы знаете, почему отфильтрованный вывод всегда начинается с нуля? Можно ли сопоставить его с фактическим входным значением x[0]? Я пробовал аналогичные вещи с фильтром нижних частот Cheby1, и у меня возникла та же проблема.
LWZ 03
2
@LWZ: используйте функцию scipy.signal.lfilter_ziи ziаргумент для lfilter. Подробнее см. В строке документации для lfilter_zi. TL; DR? Просто смени y = lfilter(b, a, data)на zi = lfilter_zi(b, a); y, zo = lfilter(b, a, data, zi=zi*data[0]). (Но это может не иметь значения с полосовым фильтром или фильтром высоких частот.)
Уоррен Векессер
1
Я заметил, что есть сдвиг фазы на 180 градусов на выходе по отношению scipy.signal.lfiter()к исходному сигналу и signal.filtfilt()выходу, почему это? Должен ли я использовать filtfilt()вместо этого, если для меня важно время?
Джейсон
1
Это фазовая задержка фильтра на этой частоте. Фазовая задержка синусоиды через фильтр Баттерворта нелинейно зависит от частоты. Для нулевой задержки фазы да, можно использовать filtfilt(). Мой ответ здесь включает пример использования, filtfilt()чтобы избежать задержки, вызванной фильтром.
Уоррен Векессер
1
Привет, Джейсон, я рекомендую задавать вопросы о теории обработки сигналов на dsp.stackexchange.com . Если у вас есть вопрос о написанном вами коде, который не работает должным образом, вы можете задать новый вопрос здесь, в stackoverflow.
Уоррен Векессер
41

Метод построения фильтра в принятом ответе верен, но имеет недостаток. Полосовые фильтры SciPy, разработанные с использованием b, a, нестабильны и могут приводить к ошибочным фильтрам на более высоких порядках фильтров .

Вместо этого используйте вывод конструкции фильтра 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)
user13107
источник
+1, потому что теперь это лучший способ во многих случаях. Как и в комментариях к принятому ответу, также можно устранить фазовую задержку с помощью прямой-обратной фильтрации. Просто замените sosfiltна sosfiltfilt.
Майк,
@Mike и user13107. Эта же ошибка влияет на фильтры Баттерворта высоких и низких частот? И решение такое же?
dewarrn1
3
@ dewarrn1 На самом деле неправильно называть это "ошибкой"; алгоритм реализован правильно, но по своей сути нестабилен, так что это просто плохой выбор алгоритма. Но да, это влияет любой фильтр более высокого порядка - не только на верхние или нижние частоты и не только на фильтры Баттерворта, но и на другие фильтры, такие как Чебышева и так далее. В любом случае, в общем, лучше всегда выбирать sosвыход, потому что это всегда позволит избежать нестабильности. И если вам не нужна обработка в реальном времени, вы должны всегда использовать sosfiltfilt.
Майк
Извините, я давно не заметил этого ответа! @ user13107, да, представление передаточной функции (или ba) линейного фильтра имеет некоторые серьезные числовые проблемы, когда порядок фильтра большой. Даже фильтры относительно низкого порядка могут иметь проблемы, когда желаемая полоса пропускания мала по сравнению с частотой дискретизации. Мой первоначальный ответ был написан до того, как представление SOS было добавлено в SciPy, и доfs аргумент был добавлен ко многим функциям в scipy.signal. Ответ давно пора обновить.
Уоррен Векессер,
любая помощь для этого? stackoverflow.com/q/60866193/5025009
seralouk
4

Для полосового фильтра 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)

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

Sizzzzlerz
источник
Пытался реализовать, но чего-то все равно не хватает. Одна вещь заключается в том, что gpass=0.0возникает ошибка деления на ноль, поэтому я изменил ее на 0,1, и ошибка прекратилась. Кроме того, документы butterговорят: Passband and stopband edge frequencies, normalized from 0 to 1 (1 corresponds to pi radians / sample).я сомневаюсь, правильно ли ваш ответ сделал расчеты, поэтому я все еще работаю над этим и скоро дам некоторые отзывы.
heltonbiker
(также, хотя у меня wsи wpесть по два элемента, фильтр выполняет только низкие или высокие btype
частоты
1
Согласно документации на docs.scipy.org/doc/scipy/reference/generated/… , buttord разрабатывает фильтры низких, высоких и полосовых частот. Что касается gpass, я предполагаю, что buttord не допускает ослабления 0 дБ в полосе пропускания. Затем установите его на какое-нибудь ненулевое значение.
sizzzzlerz