ОБНОВИТЬ:
Я нашел 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
Документы и примеры сбивают с толку и непонятны, но я хотел бы реализовать форму, представленную в благодарности, помеченной как «для пропускной способности». Знаки вопроса в комментариях показывают, что я просто скопировал какой-то пример, не понимая, что происходит.
Я не инженер-электрик или ученый, я просто разработчик медицинского оборудования, которому необходимо выполнить довольно простую полосовую фильтрацию сигналов ЭМГ.
источник
Ответы:
Вы можете пропустить использование батторда и вместо этого просто выбрать порядок для фильтра и посмотреть, соответствует ли он вашему критерию фильтрации. Чтобы сгенерировать коэффициенты фильтра для полосового фильтра, задайте 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()
Вот графики, которые генерирует этот скрипт:
источник
x[0]
? Я пробовал аналогичные вещи с фильтром нижних частот Cheby1, и у меня возникла та же проблема.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])
. (Но это может не иметь значения с полосовым фильтром или фильтром высоких частот.)scipy.signal.lfiter()
к исходному сигналу иsignal.filtfilt()
выходу, почему это? Должен ли я использоватьfiltfilt()
вместо этого, если для меня важно время?filtfilt()
. Мой ответ здесь включает пример использования,filtfilt()
чтобы избежать задержки, вызванной фильтром.Метод построения фильтра в принятом ответе верен, но имеет недостаток. Полосовые фильтры 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)
источник
sosfilt
наsosfiltfilt
.sos
выход, потому что это всегда позволит избежать нестабильности. И если вам не нужна обработка в реальном времени, вы должны всегда использоватьsosfiltfilt
.fs
аргумент был добавлен ко многим функциям вscipy.signal
. Ответ давно пора обновить.Для полосового фильтра 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)
Длина результирующего фильтра будет зависеть от глубины полос задерживания и крутизны кривой отклика, которая определяется разницей между частотой среза и частотой полосы задерживания.
источник
gpass=0.0
возникает ошибка деления на ноль, поэтому я изменил ее на 0,1, и ошибка прекратилась. Кроме того, документыbutter
говорят:Passband and stopband edge frequencies, normalized from 0 to 1 (1 corresponds to pi radians / sample).
я сомневаюсь, правильно ли ваш ответ сделал расчеты, поэтому я все еще работаю над этим и скоро дам некоторые отзывы.ws
иwp
есть по два элемента, фильтр выполняет только низкие или высокиеbtype