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)

我正在对一个100kHz的电压信号应用此滤波器，对于截止频率>=60Hz的情况效果很好。但是它在低于这个值时不起作用。我想切掉所有低于10 Hz的频率。请问我的错误在哪里？我观察到滤波器阶数越低，则截止频率越低。

样本信号可以在这里找到。

- ChrisG

---

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

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

你能举一些例子来说明你如何使用它，观察到的行为是什么样的，期望的行为又是什么样的吗？寻求调试帮助的问题需要清晰的问题陈述。"这里可以工作"与"那里不能工作"不足够清晰明了。 - moooeeeep

@moooeeeep：我添加了一个示例信号和一个应用滤波器的示例。 - ChrisG

Answer 1

I hope this can help you:

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

- Konstantin Purtov

---

谢谢你的帮助！现在我有点困惑。截止频率和阶数之间有什么关系？ - ChrisG

没有关系。 "截止频率" 意味着剪切频率。 "阶数" 确定滤波器的精度。 - Konstantin Purtov

这就是我想的。但如果你将顺序设置为“30”，结果会发生很大变化。而对于高阶，低截止频率不起作用。 - ChrisG

这取决于许多因素。您可以使用Matlab滤波器设计工具进行设计和尝试。我认为，使用LMS滤波器是去除非常低频率的最佳方法（https://dev59.com/t2Ml5IYBdhLWcg3woYIv）。 - Konstantin Purtov

你好。你有没有想过如何使用mp3文件代替正弦波？我想从mp3和wav文件中删除所有低于某个阈值的频率。 - It's a trap

@It'satrap 你好。是的，你可以将mp3文件读取为信号，使用过滤器并保存回去。在此处有一个很好的例子，介绍如何读取mp3文件并使用过滤器：https://dev59.com/cW_Xa4cB1Zd3GeqP0266 - Konstantin Purtov

Answer 2

由于我的声望低，无法对您的问题进行评论 - “截止频率和滤波器阶数之间的关系是什么？”这不是对您最初问题的回答。 对于一个FIR滤波器，对于给定的截止频率，脉冲响应图（| H（f）| vs f）的斜率对于更高阶滤波器更陡峭。因此，为了实现对不需要的频率范围的更高衰减，应增加滤波器阶数。但是，当滤波器阶数非常高时，脉冲响应是理想的矩形函数吗？您将看到类似于数字通信中的干扰（ISI）的效果。当截止频率与采样频率的比例变小时，这种效应的强度会增加（考虑在频域中盒函数的宽度与sinc函数的主瓣宽度之间的关系）。 我第一次观察到这一点是当我试图在TI DSP微控制器上实现非常窄的带低通IIR滤波器时。 TI库将滤波器实现为级联双二次结构，以处理众所周知的截断效应。这仍然没有解决问题，因为问题不仅仅是由于截断造成的。我解决这个问题的方法是使用一个抗混叠滤波器，然后对输入信号进行下采样，然后进行所需的低通IIR滤波器。 我知道您正在实现HPF，这是在频域中翻译的LPF。希望这回答了您的一些问题。如果下采样适用于您，请告诉我。

- Sriharsha Madala

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

- Cloud Cho