巴特沃斯滤波器——python实现

butter（）函数是求Butterworth数字滤波器的系数向量，在求出系数后对信号进行滤波时需要用scipy.signal.filtfilt()。
需要安装scipy包。

函数butter()

设计滤波器就是设计滤波器系数[B,A]。
[b,a]=butter(n,Wn)，根据阶数n和归一化截止频率Wn计算ButterWorth滤波器分子分母系数（b为分子系数的矢量形式，a为分母系数的矢量形式）。

[B,A] = BUTTER(N,Wn,‘high’) 高通滤波器
[B,A] = BUTTER(N,Wn,‘low’) 低通滤波器
[B,A] = BUTTER(N,Wn,‘stop’) 带阻滤波器 Wn = [W1 W2].
[B,A] = BUTTER(N,Wn) 带通滤波器
带通滤波器：它允许一定频段的信号通过，抑制低于或高于该频段的信号、干扰和噪声；
带阻滤波器：它抑制一定频段内的信号，允许该频段以外的信号通过。

其中Wn是归一化频率，具体计算方法是（2*截止频率)/采样频率（也就是除以fs/2）。
当设计低通和高通时，Wn是一个值，表示截止频率；
当设计带通和带阻时，Wn是一个二个元素的数组，表示通带或阻带的上下截止频率。频率的归一化是对fs/2进行归一。
对于原始信号x。
比如说你的采样频率fs=1000Hz，设计一个8阶、通带为fc1=100,fc2=200Hz的带通滤波器：
[b,a]=butter(8,[0.2 0.4])=butter(8,fc1/fa fc2/fa])
这里fa=fs/2，fa是分析频率
得到滤波器系数后，就可以直接用了。
y=filter(B,A,x)

scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘ba’)

输入参数：
N:滤波器的阶数
Wn：归一化截止频率。计算公式Wn=2*截止频率/采样频率。（注意：根据采样定理，采样频率要大于两倍的信号本身最大的频率，才能还原信号。截止频率一定小于信号本身最大的频率，所以Wn一定在0和1之间）。当构造带通滤波器或者带阻滤波器时，Wn为长度为2的列表。

btype : 滤波器类型{‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’},

output : 输出类型{‘ba’, ‘zpk’, ‘sos’},

输出参数：

b，a: IIR滤波器的分子（b）和分母（a）多项式系数向量。output=‘ba’

z,p,k: IIR滤波器传递函数的零点、极点和系统增益. output= ‘zpk’

sos: IIR滤波器的二阶截面表示。output= ‘sos’

函数scipy.signal.filtfilt()

scipy.signal.filtfilt(b, a, x, axis=-1, padtype=‘odd’, padlen=None, method=‘pad’, irlen=None**)**

输入参数：
b: 滤波器的分子系数向量
a: 滤波器的分母系数向量
x: 要过滤的数据数组。（array型）

（以下参数可以省略）
axis: 指定要过滤的数据数组x的轴

padtype: 必须是“奇数”、“偶数”、“常数”或“无”。这决定了用于过滤器应用的填充信号的扩展类型。{‘odd’, ‘even’, ‘constant’, None}

padlen：在应用滤波器之前在轴两端延伸X的元素数目。此值必须小于要滤波元素个数- 1。（int型或None）

method：确定处理信号边缘的方法。当method为“pad”时，填充信号；填充类型padtype和padlen决定，irlen被忽略。当method为“gust”时，使用古斯塔夫森方法，而忽略padtype和padlen。{“pad” ，“gust”}

irlen：当method为“gust”时，irlen指定滤波器的脉冲响应的长度。如果irlen是None，则脉冲响应的任何部分都被忽略。对于长信号，指定irlen可以显著改善滤波器的性能。（int型或None）

输出参数：
y:滤波后的数据数组

带通滤波代码：

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


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    fa = 0.5 * fs
    low = lowcut / fa
    high = highcut / fa
    b, a = butter(order, [low, high], btype='band')
    y = filtfilt(b, a, data)
    return y


# 生成噪声信号
fs = 1000  # 采样率
t = np.arange(0, 1, 1/fs)  # 时间序列

x = 5*np.sin(2*np.pi*50*t) + 2*np.sin(2*np.pi*120*t) + 0.5*np.sin(2*np.pi*200*t)
noise = 0.5*np.random.randn(len(t))
y = x + noise

# 设计巴特沃斯滤波器，滤除50 Hz以下和200 Hz以上的信号
filtered_signal = butter_bandpass_filter(x, 50, 200, fs, order=6)





# 绘制原始信号和滤波后的信号的频谱图和波形图
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq

# 频谱图
fig, axs = plt.subplots(2, 1)
axs[0].plot(t, y, 'b', label='raw signal')
axs[0].plot(t, filtered_signal, 'g', linewidth=2, label='filtered signal')
axs[0].set_xlabel('Time [s]')
axs[0].set_ylabel('Amplitude')
axs[0].legend()

n = len(y)
f = fftfreq(n, 1/fs)
Y = fft(y)/n
filtered_Y = fft(filtered_signal)/n
axs[1].semilogy(f[:n//2], np.abs(Y[:n//2]), 'b', label='raw signal')
axs[1].semilogy(f[:n//2], np.abs(filtered_Y[:n//2]), 'g', linewidth=2, label='filtered signal')
axs[1].set_xlabel('Frequency [Hz]')
axs[1].set_ylabel('Amplitude')
axs[1].legend()
plt.show()

# 波形图
plt.plot(t, y, 'b', label='raw signal')
plt.plot(t, filtered_signal, 'g', linewidth=2, label='filtered signal')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()
plt.show()

低通滤波代码：

from scipy.signal import butter, filtfilt
import pandas as pd
from pylab import *
import mplcursors
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator
#显示中文
mpl.rcParams['font.sans-serif'] = ['SimHei']


def filter_data(WK):
    N = 4  # Filter order
    Wn = 0.11  # Cutoff frequency
    B, A = butter(N, Wn, output='ba')
    smooth_data = filtfilt(B, A, WK)
    return smooth_data

df = pd.read_csv('PL19-3-E35ST11.csv')

plt.subplot(1, 2, 1)
x = df['PROD_TIME']
y = df['PUMP_CURRENT']
plt.title("滤波前", fontsize=8)
plt.plot(x, y, c="g", label="CURRENT")
plt.ylabel("PUMP_CURRENT")
plt.xlabel("PROD_TIME")

plt.gca().xaxis.set_major_locator(MultipleLocator(df.shape[0]/10))
plt.tick_params(axis='x', rotation=45)
plt.legend()


plt.subplot(1, 2, 2)
x = df['PROD_TIME']
y = df['PUMP_CURRENT']
y2 = filter_data(y)
plt.title("滤波后", fontsize=8)
plt.plot(x, y2, c="r", label="CURRENT")
plt.ylabel("PUMP_CURRENT")
plt.xlabel("PROD_TIME")

plt.gca().xaxis.set_major_locator(MultipleLocator(df.shape[0]/10))
plt.tick_params(axis='x', rotation=45)
plt.legend()


# 放大拖动功能
# enable zooming on both subplots
mplcursors.cursor(hover=True).connect("add", lambda sel: dir(sel))
plt.tight_layout(pad=1.08)
# 显示图像
plt.show()