Как реализовать непрерывный фильтр высоких / низких частот в python?

Question

У меня есть приложение, которое принимает входящие сигналы с фиксированной частотой. Я пытался реализовать фильтр для этого входящего сигнала без необходимости сохранять N временных шагов и выполнять на нем функцию фильтрации. То, что я хотел бы сделать, похоже на одномерный фильтр Калмана, в котором я обновляю свое текущее состояние с новым наблюдением. Я не математик, и поэтому страницы Википедии о том, возможно ли это, для меня совершенно непонятны. Кроме того, ответы StackOverflow (я нашел) в этом домене касаются только того, есть ли у вас доступный раздел из N временных шагов сигнала и как выполнить фильтрацию для этого, и такую ​​фильтрацию я могу выполнить без проблем.

Ниже я привел некоторый фиктивный код, который иллюстрирует типы функций, которые я пытался написать.

import numpy as np
import matplotlib.pyplot as plt

def high_pass(new, previous, cutoff=20):
    # Howto?
    return 0

def low_pass(new, previous, cutofff=20):
    # Howto?
    return 0

def continuous_filter(vs):
    prev_highpass, prev_lowpass = 0, 0
    for v in vs:
        prev_highpass = high_pass(v, prev_highpass)
        prev_lowpass = low_pass(v, prev_lowpass)
        yield prev_highpass, prev_lowpass


np.random.seed(21)
sec_duration = 10.0
time_resolution = 1e3
dt = 1/time_resolution
steps = int(sec_duration * time_resolution)
signal = np.sum([np.sin(np.linspace(0, np.random.randint(10, 100), steps)) * np.random.random() for _ in range(3)], axis=0)

filt_signals = np.array([[high, low] for high, low in continuous_filter(signal)])

plt.plot(signal, label="Input signal")
plt.plot(filt_signals[:, 0], label="High-pass")
plt.plot(filt_signals[:, 1], label="Low pass")
plt.legend()
plt.show()

Кто-нибудь, кто может сказать мне, возможно ли это? Я смотрел документацию SciPy , но не понимаю.

Answer 1

Во-первых, вам нужна постоянная времени для вашего фильтра из отсечки: alpha = dt / (RC + dt) и cutoff = 1 / (2 * pi * RC).

Этот коэффициент нужен вам для вычисления следующего отфильтрованного значения:

def low_pass(x_new, y_old, cutoff=20):
    alpha = dt / (dt + 1 / (2 * np.pi * cutoff))
    y_new = x_new * alpha + (1 - alpha) * y_old
    return y_new

Из википедии: низкочастотный .

def high_pass(x_new, x_old, y_old, cutoff=20):
    alpha = dt / (dt + 1 / (2 * np.pi * cutoff))
    y_new = alpha * (y_old + x_new - x_old)
    return y_new 

Из википедии: высокочастотный .

def continuous_filter(xs):
    prev_highpass, prev_lowpass = 0, 0
    x_prev = 0  # need initializatoin for highpass
    y_prev_high = 0  # initialization
    y_prev_low = 0  # initialization

    for x in xs:
        y_prev_high = high_pass(x, x_prev, y_prev_high)
        y_prev_low = low_pass(x, y_prev_low)
        x_prev = x
        yield y_prev_high, y_prev_low


np.random.seed(21)
sec_duration = 10.0
time_resolution = 1e3
dt = 1/time_resolution

steps = int(sec_duration * time_resolution)
signal = np.sum([np.sin(np.linspace(0, np.random.randint(10, 100), steps)) * np.random.random() for _ in range(3)], axis=0)

filt_signals = np.array([[high, low] for high, low in continuous_filter(signal)])