c++实现巴特沃斯带通滤波器
时间: 2023-05-30 07:07:26 浏览: 67
巴特沃斯带通滤波器是一种常用的数字滤波器,用于去除信号中的高频和低频噪声,保留中心频率附近的信号。其数学模型为:
$$
H(z) = \frac{K}{\prod_{i=1}^n \left(1 - z^{-1}p_i\right)\prod_{j=1}^n \left(1 - z^{-1}q_j\right)}
$$
其中,$n$ 是滤波器的阶数,$p_i$ 和 $q_j$ 是滤波器的极点和零点,$K$ 是滤波器的增益。
在实现巴特沃斯带通滤波器时,可以使用双线性变换将其转换为离散时间滤波器。具体步骤如下:
1. 根据中心频率 $f_c$ 和通带宽度 $B$,计算出模拟滤波器的极点和零点。
2. 使用双线性变换将模拟滤波器转换为离散时间滤波器。
3. 根据所需的滤波器阶数 $n$,计算出滤波器的增益 $K$。
4. 实现离散时间滤波器,即将输入信号按时间序列依次送入滤波器,得到输出信号。
以下是一个 Python 实现巴特沃斯带通滤波器的例子:
```python
import numpy as np
from scipy.signal import butter, filtfilt
# 定义采样频率和通带边界频率
fs = 1000.0 # 采样频率
f_low = 10.0 # 通带低频
f_high = 100.0 # 通带高频
B = f_high - f_low # 通带宽度
# 计算模拟滤波器的极点和零点
omega_low = 2.0 * np.pi * f_low
omega_high = 2.0 * np.pi * f_high
omega_c = np.sqrt(omega_low * omega_high)
Q = omega_c / B
p1 = -np.sinh(np.arcsinh(1.0 / Q) / 2.0) * np.exp(1j * np.pi / 4.0)
p2 = -np.sinh(np.arcsinh(1.0 / Q) / 2.0) * np.exp(-1j * np.pi / 4.0)
z1 = -1.0
z2 = 1.0
# 使用双线性变换将模拟滤波器转换为离散时间滤波器
b, a = butter(2, [2.0 * fs * np.tan(0.5 * omega_low / fs), 2.0 * fs * np.tan(0.5 * omega_high / fs)], btype='bandpass', analog=False)
# 计算滤波器的增益
w, h = np.linspace(0.0, np.pi, 1000, endpoint=True, retstep=True)
H = np.exp(-1j * w * (len(a) - 1)) * np.polyval(b, np.exp(1j * w)) / np.polyval(a, np.exp(1j * w))
K = np.abs(np.sqrt(np.abs(H[0] * H[-1])))
# 实现离散时间滤波器
def butter_bandpass(data, lowcut, highcut, fs, order=2):
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut / nyq
b, a = butter(order, [low, high], btype='band')
y = filtfilt(b, a, data)
return y
# 测试
import matplotlib.pyplot as plt
t = np.linspace(0, 1, 1000)
x = np.sin(2.0 * np.pi * 50.0 * t) + np.sin(2.0 * np.pi * 150.0 * t) + 0.5 * np.random.randn(len(t))
y = butter_bandpass(x, f_low, f_high, fs, order=2)
fig, ax = plt.subplots(2, 1, figsize=(10, 6))
ax[0].plot(t, x)
ax[0].set_title('Input signal')
ax[0].set_xlabel('Time [s]')
ax[0].set_ylabel('Amplitude')
ax[1].plot(t, y)
ax[1].set_title('Output signal')
ax[1].set_xlabel('Time [s]')
ax[1].set_ylabel('Amplitude')
plt.tight_layout()
plt.show()
```
该例子中,我们首先定义了采样频率和通带边界频率,然后计算出模拟滤波器的极点和零点,使用双线性变换将其转换为离散时间滤波器。接着,我们计算滤波器的增益,实现离散时间滤波器,最后用一个包含两个子图的图表显示输入信号和输出信号。
相关推荐
![dll](https://img-home.csdnimg.cn/images/20210720083646.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![-](https://csdnimg.cn/download_wenku/file_type_column_c1.png)
![rar](https://img-home.csdnimg.cn/images/20210720083606.png)