如何使用 Scipy – Python实现 IIR 带通巴特沃斯滤波器？

在数字信号处理中，滤波器是必不可少的工具之一，它可以帮我们对信号进行去噪、频谱分析、信号重构等操作。本文将介绍如何使用 Scipy 库中的 butter 函数，实现 IIR（Infinite Impulse Response，无限脉冲响应）带通巴特沃斯滤波器。

步骤一：导入必要的库

Scipy 是 Python 中的一个科学计算库，其中包含了一些常用的数字信号处理函数。在使用之前，需要先导入 Scipy 库，以及一些常用的 Python 库，如 numpy 和 matplotlib。

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter, freqz

步骤二：定义滤波器参数

在定义滤波器参数之前，需要先确定所需的滤波器类型、截止频率、通带和阻带的增益等参数。这里我们以 4 阶带通巴特沃斯滤波器为例。

# 滤波器阶数
order = 4

# 通带和阻带增益
gain = 10

# 截止频率
fs = 1000  # 采样频率
lowcut = 20  # 通带下限频率
highcut = 100  # 通带上限频率
Wn = [lowcut/fs*2, highcut/fs*2]  # 截止频率归一化（Nyquist 频率为采样频率的一半）

步骤三：计算滤波器系数

接下来，使用 Scipy 库中的 butter 函数计算滤波器的系数。

# 计算滤波器系数
b, a = butter(order, Wn, btype='bandpass', analog=False, output='ba')

其中，order 是滤波器的阶数，Wn 是归一化截止频率，btype 是滤波器类型，analog 表示是否为模拟滤波器，output 表示返回的系数类型。

步骤四：绘制滤波器频率响应

可以使用 freqz 函数绘制滤波器的频率响应曲线。

# 绘制频率响应
w, h = freqz(b, a)
fig = plt.figure(figsize=(8, 6))
plt.plot(fs*w/(2*np.pi), 20*np.log10(np.abs(h)), 'b')
plt.axvline(lowcut, color='k')
plt.axvline(highcut, color='k')
plt.ylim(-gain*2, gain)
plt.xlim(0, fs/2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain (dB)')
plt.grid(True)
plt.show()

其中，freqz 函数返回滤波器的频率响应，w 是频率数组，h 是复数数组。

步骤五：应用滤波器

最后，使用 lfilter 函数将滤波器应用到信号上。这里我们以一个正弦波作为例子。

# 生成正弦波信号
t = np.linspace(0, 1, fs, endpoint=False)
sig = np.sin(2*np.pi*80*t) + 0.5*np.sin(2*np.pi*200*t)

# 应用滤波器
filtered = lfilter(b, a, sig)

# 绘制原始信号和滤波后的信号
fig = plt.figure(figsize=(8, 6))
plt.plot(t, sig, 'b-', label='Signal')
plt.plot(t, filtered, 'g-', linewidth=2, label='Filtered')
plt.xlim(0, 0.01)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend(loc='best')
plt.show()

其中，lfilter 函数将滤波器系数 b 和 a 应用到信号 sig 上，返回滤波后的信号 filtered。

这样就完成了一个 IIR 带通巴特沃斯滤波器的设计和实现。可以根据需要修改滤波器类型、阶数、截止频率等参数，以满足不同的应用需求。

完整代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, lfilter, freqz

# 滤波器阶数
order = 4

# 通带和阻带增益
gain = 10

# 截止频率
fs = 1000  # 采样频率
lowcut = 20  # 通带下限频率
highcut = 100  # 通带上限频率
Wn = [lowcut/fs*2, highcut/fs*2]  # 截止频率归一化（Nyquist 频率为采样频率的一半）

# 计算滤波器系数
b, a = butter(order, Wn, btype='bandpass', analog=False, output='ba')

# 绘制频率响应
w, h = freqz(b, a)
fig = plt.figure(figsize=(8, 6))
plt.plot(fs*w/(2*np.pi), 20*np.log10(np.abs(h)), 'b')
plt.axvline(lowcut, color='k')
plt.axvline(highcut, color='k')
plt.ylim(-gain*2, gain)
plt.xlim(0, fs/2)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain (dB)')
plt.grid(True)
plt.show()

# 生成正弦波信号
t = np.linspace(0, 1, fs, endpoint=False)
sig = np.sin(2*np.pi*80*t) + 0.5*np.sin(2*np.pi*200*t)

# 应用滤波器
filtered = lfilter(b, a, sig)

# 绘制原始信号和滤波后的信号
fig = plt.figure(figsize=(8, 6))
plt.plot(t, sig, 'b-', label='Signal')
plt.plot(t, filtered, 'g-', linewidth=2, label='Filtered')
plt.xlim(0, 0.01)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.legend(loc='best')
plt.show()