如何在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 documentation 但我不明白。

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

来自维基百科：low-pass。

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 

来自维基百科：high-pass。

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)])

Answer 2

这是一个很好的答案，多利安，谢谢。经过多年的跟踪堆栈溢出，这将是我的第一个贡献。它是次要的，希望它有助于扩展用例。尽量保留原始符号，以便可追溯。

只需要两个导入

'''python
import numpy as np
import matplotlib.pyplot as plt

对 rc_high_pass alpha 值的小幅修正

def rc_low_pass(x_new, y_old, sample_rate_hz, lowpass_cutoff_hz):
    dt = 1/sample_rate_hz
    rc = 1/(2*np.pi*lowpass_cutoff_hz)
    alpha = dt/(rc + dt)
    y_new = x_new * alpha + (1 - alpha) * y_old
    return y_new


def rc_high_pass(x_new, x_old, y_old, sample_rate_hz, highpass_cutoff_hz):
    dt = 1/sample_rate_hz
    rc = 1/(2*np.pi*highpass_cutoff_hz)
    alpha = rc/(rc + dt)
    y_new = alpha * (y_old + x_new - x_old)
    return y_new


def rc_filters(xs, sample_rate_hz, 
               highpass_cutoff_hz, 
               lowpass_cutoff_hz):
    # Initialize. This can be improved to match wikipedia.
    x_prev = 0
    y_prev_high = 0
    y_prev_low = 0

    for x in xs:
        y_prev_high = rc_high_pass(x, x_prev, y_prev_high, sample_rate_hz, 
                                   highpass_cutoff_hz)
        y_prev_low = rc_low_pass(x, y_prev_low, sample_rate_hz, 
                                 lowpass_cutoff_hz)
        x_prev = x
        yield y_prev_high, y_prev_low

现在回到主要部分，我使用以赫兹和秒为单位的物理单位。请原谅以 2 为底的数字表示，它有助于为 ffts 奠定基础。

if __name__ == "__main__":
    """
    # RC filters for continuous signals
    """
    sample_rate = 2**13  # Close to 8 kHz
    duration_points = 2**10
    sec_duration = duration_points/sample_rate

    frequency_low = sample_rate/2**9
    frequency_high = sample_rate/2**3

    # Design the cutoff
    number_octaves = 3
    highpass_cutoff = frequency_high/2**number_octaves
    lowpass_cutoff = frequency_low*2**number_octaves

    print('Two-tone test')
    print('Sample rate, Hz:', sample_rate)
    print('Record duration, s:', sec_duration)
    print('Low, high tone frequency:', frequency_low, frequency_high)

    time_s = np.arange(duration_points)/sample_rate

    sig = np.sin(2*np.pi*frequency_low*time_s) + \
          np.sin(2*np.pi*frequency_high*time_s)

    filt_signals = np.array([[high, low]
                             for high, low in
                             rc_filters(sig, sample_rate,
                                        highpass_cutoff, lowpass_cutoff)])

输入信号不同，更接近我的典型用例，结果图显示通常的 RC 滤波器不是很快。这是一个经典的解决方案，在这里有它很好。

plt.plot(sig, label="Input signal")
plt.plot(filt_signals[:, 0], label="High-pass")
plt.plot(filt_signals[:, 1], label="Low-pass")
plt.title("RC Low-pass and High-pass Filter Response")
plt.legend()
plt.show()

Two-tone input, 6 octave separation; substantial spectral leakage.