Python 高通滤波器

Posted 2023-03-12

【中文标题】Python 高通滤波器

【英文标题】：Python High Pass Filter

【发布时间】：2016-12-26 05:14:25

【问题描述】：

我使用以下代码在 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.

【问题讨论】：

你在哪里设置你的通带？你有你可以粘贴的公式吗？和你的filtfilt() 函数。和你的butter() 函数

@Daniel Lee：我编辑了我的问题。 butter() 和 filtfilt() 函数来自 scipy.signal。

你能举一些例子，你是如何使用它的，观察到的行为是什么样的，以及预期的行为会是什么？寻求调试帮助的问题需要明确的问题陈述。 在这里工作得很好 vs. 在那里不工作不够清楚。

@moooeeeep：我添加了一个样本信号和一个如何应用过滤器的示例。

【参考方案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()

【讨论】：

感谢您的帮助！现在我有点困惑。 cutoff和order是什么关系？

没有关系。 “截止”是指截止频率。 “顺序”决定过滤器的精度。

我就是这么想的。但是，如果您将订单设置为“30”。结果变化很大。对于高阶，低截止频率不起作用。

这取决于很多因素。您可以使用 Matlab 滤波器设计工具包设计和尝试您想要的东西。我认为 LMS 滤波器是去除极低频的最佳方法 (***.com/questions/18252160/…)

嗨。你知道我如何使用 mp3 文件而不是正弦波。我想从 mp3 和 wav 文件中删除低于某个阈值的所有频率

【参考方案2】：

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

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

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

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

【参考方案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()