Python 高通滤波器

Question

我使用以下代码在 python 中实现了一个高通滤波器：

from scipy.signal import butter, filtfilt
import numpy as np

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = filtfilt(b, a, data)
    return y

rawdata = np.loadtxt('sampleSignal.txt', skiprows=0)
signal = rawdata
fs = 100000.0

cutoff = 100
order = 6
conditioned_signal = butter_highpass_filter(signal, cutoff, fs, order)

我将此滤波器应用于 100 kHz 电压信号，它适用于 >= 60 Hz 的截止频率。但它在下面不起作用。我想切断所有低于 10 Hz 的频率。任何提示我的错误在哪里？我观察到的是滤波器的阶数越低，截止频率就越低。

The sample Signal can be found here.

Answer 1

希望对你有帮助：

import numpy as np
import pandas as pd
from scipy import signal
import matplotlib.pyplot as plt
def sine_generator(fs, sinefreq, duration):
    T = duration
    nsamples = fs * T
    w = 2. * np.pi * sinefreq
    t_sine = np.linspace(0, T, nsamples, endpoint=False)
    y_sine = np.sin(w * t_sine)
    result = pd.DataFrame({ 
        'data' : y_sine} ,index=t_sine)
    return result

def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y

fps = 30
sine_fq = 10 #Hz
duration = 10 #seconds
sine_5Hz = sine_generator(fps,sine_fq,duration)
sine_fq = 1 #Hz
duration = 10 #seconds
sine_1Hz = sine_generator(fps,sine_fq,duration)

sine = sine_5Hz + sine_1Hz

filtered_sine = butter_highpass_filter(sine.data,10,fps)

plt.figure(figsize=(20,10))
plt.subplot(211)
plt.plot(range(len(sine)),sine)
plt.title('generated signal')
plt.subplot(212)
plt.plot(range(len(filtered_sine)),filtered_sine)
plt.title('filtered signal')
plt.show()

Answer 2

由于我的声誉很低，我无法评论您的问题 - “截止和过滤顺序之间的关系是什么？”这不是您原来问题的答案。

对于 FIR 滤波器，对于给定的截止频率，脉冲响应图的斜率（|H(f)| vs f）对于更高阶的滤波器更陡峭。因此，为了在不需要的频率范围内实现更高的衰减，您需要增加滤波器阶数。但是当滤波器阶数如此之高以至于脉冲响应是一个理想的盒函数时会发生什么？您将看到类似符号间干扰（数字通信中的 ISI）的效果。当截止频率与采样频率的比率变小时，这种效应的强度会增加（想想频域中盒函数的宽度与 sinc 函数的主瓣宽度之间的关系）。

当我尝试在 TI DSP 微控制器上实现非常窄带的低通 IIR 滤波器时，我首先观察到了这一点。 TI 库将滤波器实现为级联双四边形结构，以处理众所周知的截断效应。这仍然没有解决问题，因为问题不仅仅是由于截断。我解决这个问题的方法是使用抗混叠滤波器，然后对输入信号进行下采样，然后使用我想要的低通 IIR 滤波器。

我了解到您正在实施 HPF，它是在频域中转换的 LPF。希望这能回答你的一些问题。让我知道下采样是否适合您。

Answer 3

我在Konstantin Purtov的代码上方添加了验证功能，因此您可以看到顺序和截止频率之间的关系。代码大部分来自Warren Weckesser。

import scipy
import sys  
from scipy import signal
from scipy import pi
from scipy.io.wavfile import write
import matplotlib.pyplot as plt
import numpy as np    
from scipy.signal import butter, lfilter, freqz   


# ----- ----- ----- -----    
def butter_highpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = signal.butter(order, normal_cutoff, btype='high', analog=False)
    return b, a

def butter_highpass_filter(data, cutoff, fs, order=5):
    b, a = butter_highpass(cutoff, fs, order=order)
    y = signal.filtfilt(b, a, data)
    return y


# ----- -----    
# (1)
def foo(sel):
    if (sel == 1):
        # Filter requirements.
        order = 6
        fs = 300.0  # sample rate, Hz
        cutoff = 10  # desired cutoff frequency of the filter, Hz

        # Get the filter coefficients so we can check its frequency response.
        b, a = butter_highpass(cutoff, fs, order)

        # Plot the frequency response.
        w, h = freqz(b, a, worN=8000)
        plt.subplot(2, 1, 1)
        plt.plot(0.5 * fs * w / np.pi, np.abs(h), 'b')
        plt.plot(cutoff, 0.5 * np.sqrt(2), 'ko')
        plt.axvline(cutoff, color='k')
        plt.xlim(0, 0.5 * fs)
        plt.title("High Filter Frequency Response")
        plt.xlabel('Frequency [Hz]')
        plt.grid()

        # Demonstrate the use of the filter.
        # First make some data to be filtered.
        T = 0.5  # seconds
        n = int(T * fs)  # total number of samples
        t = np.linspace(0, T, n, endpoint=False)
        # "Noisy" data.  We want to recover the 20 Hz signal from this.
        data = np.sin(1.2 * 2 * np.pi * t) + 1.5 * np.cos(5 * 2 * np.pi * t) + 0.5 * np.sin(20.0 * 2 * np.pi * t)

        # Filter the data, and plot both the original and filtered signals.
        y = butter_highpass_filter(data, cutoff, fs, order)

        plt.subplot(2, 1, 2)
        plt.plot(t, data, 'b-', label='data')
        plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
        plt.xlabel('Time [sec]')
        plt.grid()
        plt.legend()

        plt.subplots_adjust(hspace=0.35)
        plt.show()
    else:
        print ('Please, choose among choices, thanks.')


# ----- -----
def main():
    sel = int (sys.argv[1])
    foo(sel)


# ----- ----- ----- ----- ----- -----
if __name__ == '__main__':
    main()