python实现IIR⾼通低通，带通，带阻滤波器详解及应⽤案例Python实现数字滤波器 ⽂章⽬录
由语⾳的产⽣和感知可知，基⾳频率的范围是60到450赫兹之间，故语⾳信号采集需要提取基⾳时，需要采⽤低通滤波器来获取低频基⾳信号，在采⽤计算机采集语⾳信号时，语⾳常混有50赫兹交流混⾳，也需要采⽤⾼通滤波器将其去除，此篇设计数字滤波器，以及实现他们在语⾔中的简单应⽤。滤波器传递函数如下： 给定数字滤波器的M阶分⼦b和N阶分母a：
jw -jw -jwM
jw B(e ) b[0] + b[1]e + ... + b[M]e
H(e ) = ------ = -----------------------------------
jw -jw -jwN
A(e ) a[0] + a[1]e + ... + a[N]e
b[0]*(jw)**M + b[1]*(jw)**(M-1) + ... + b[M]
H(w) = ----------------------------------------------
a[0]*(jw)**N + a[1]*(jw)**(N-1) + ... + a[N]
1、IIR低通、⾼通、带通和带阻滤波器的设计
1.1、设计滤波器的函数
（1）巴特沃斯滤波器
巴特沃斯数字和模拟滤波器设计函数：
(N, Wn[, btype, analog, output, fs])
功能：设计⼀个N阶数字或模拟巴特沃斯滤波器，然后返回滤波器系数。
调⽤格式：
b,a=scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘ba’, fs=None)
z,p,k=scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘zpk’, fs=None)
sos=scipy.signal.butter(N, Wn, btype=‘low’, analog=False, output=‘sos’, fs=None)
参数：
N（int）（滤波器阶次）
Wn（array、int）（临界频率对于低通和⾼通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于巴特沃斯滤波器，这是增益降⾄通带的1/sqrt(2)的点“-3dB点”，对于数字滤波器，Wn与fs的单位相同，对于模拟滤波器，Wn是⾓频率rad / s） btype：（‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’optional）（滤波类型，默认值是低通其中低通、⾼通、带通、带阻分别对应‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采⽤频率）
返回值：
b（ndarray）（传递函数分⼦）
a（ndarray）（传递函数分母）
z（ndarray）（零点）
p（ndarray）（极点）
k（float）（系统增益）
sos（ndarray）（IIR滤波器⼆阶部分表⽰）
巴特沃斯滤波器在通带中具有最⼤平坦的频率响应，如果要求传递函数形式，则由于根与多项式系数之间的转换是数值敏感的运算，即使对于N> = 4，也会出现数值问题。建议使⽤SOS表⽰形式。 绘制模拟滤波器频率响应：
计算频率响应的语句：
w, h = signal.freqs(b, a)#计算模拟滤波器的频率响应，返回h为频率响应，w为h的⾓频率
w, h = signal.freqz(b, a, worN=512, whole=False, plot=None, fs=6.283185307179586, include_nyquist=False)#计算数字滤波器的频率响应，返回h为频率响应（复数表⽰）w为h的⾓频率单位与fs相同默认情况下，w归⼀化为范围[0，pi）（弧度/样本），include_nyquist包含最后⼀个频率（奈奎斯特频率）filtered = signal.sosfilt(sos, sig)#使⽤级联的⼆阶节沿⼀维过滤数据
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
b, a = signal.butter(4,100,'low', analog=True)# 4阶低通临界频率为100Hz
w, h = signal.freqs(b, a)# h为频率响应,w为频率
plt.figure(1)
plt.semilogx(w,20* np.log10(abs(h)))# ⽤于绘制折线图，两个函数的 x 轴、y 轴分别是指数型的，并且转化为分贝
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0,0.1)# 去除画布四周⽩边
plt.axvline(100, color='green')# 绘制竖线
plt.show()
t = np.linspace(0,1,1000,False)# 1秒，1000赫兹刻度
sig = np.sin(2*np.pi*10*t)+ np.sin(2*np.pi*20*t)# 合成信号
fig,(ax1, ax2)= plt.subplots(2,1, sharex=True)# 2⾏1列的图
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0,1,-2,2])# 坐标范围
⽣成由10 Hz和20 Hz组成的信号，以1 kHz采样
t = np.linspace(0,1,1000,False)# 1秒，1000赫兹刻度
sig = np.sin(2*np.pi*10*t)+ np.sin(2*np.pi*20*t)# 合成信号
fig,(ax1, ax2)= plt.subplots(2,1, sharex=True)# 2⾏1列的图
ax1.plot(t, sig)
电子除垢器ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0,1,-2,2])# 坐标范围
设计⼀个15 Hz的数字⾼通滤波器以消除10 Hz的⾳调，并将其应⽤于信号。（建议在过滤时使⽤⼆阶节格式，以避免传递函数（ba）格式出现数值错误）：
#sos = signal.butter(10, 15, 'hp', fs=1000, output='sos') #10阶，15赫兹
sos = signal.butter(10,15, btype='low', analog=False, output='sos', fs=1000)
filtered = signal.sosfilt(sos, sig)# 滤波
ax2.plot(t, filtered)
ax2.set_title('After 15 Hz high-pass filter')
ax2.axis([0,1,-2,2])
活顶尖ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()
巴特沃斯滤波器阶次选择函数：
(wp, ws, gpass, gstop[, analog, fs])
功能：返回最低阶数字或模拟巴特沃思滤波器的阶次，该滤波器在通带中的损失不超过gpass dB，在阻带中具有⾄少 gstop dB的衰减。调⽤格式：
N, Wn = signal.buttord(wp, ws, gpass, gstop, analog=False, fs=None)
参数：
Wp（float）（通带边缘频率）
Ws（float）（阻带边缘频率）
对于数字滤波器，它们的单位与fs相同。默认情况下， fs是2个半周期/样本，因此将它们从0标准化为1，其中1是奈奎斯特频率。
（因此，wp和ws在半周期/样本中。）例如：
玻璃夹胶机
低通：wp = 0.2，ws = 0.3
⾼通：wp = 0.3，ws = 0.2
带通：wp = [0.2，0.5]，ws = [0.1，0.6]
带阻：wp = [0.1，0.6]，ws = [0.2，0.5]
对于模拟滤波器，wp和ws是⾓频率（例如rad / s）。
gpass（float）（通带中的最⼤损耗dB）
gstop（float）（阻带中的最⼩衰减dB）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采⽤频率）
返回值：
N（int）（滤波器阶次）
Wn（array、int）（临界频率对于低通和⾼通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于
巴特沃斯滤波器，这是增益降⾄通带的1/sqrt(2)的点“-3dB点”，对于数字滤波器，Wn与fs的单位相同，对于模拟滤波器，Wn是⾓频率rad / s）
设计⼀个模拟带通滤波器，其通带范围在20到50 rad / s之间的3 dB之内，⽽在14 rad以下和60 rad / s以上的噪声⾄少抑制-40 dB。绘制其频率响应，以灰⾊显⽰通带和阻带约束。
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
N, Wn = signal.buttord([20,50],[14,60],3,40,True)
b, a = signal.butter(N, Wn,'band',True)
w, h = signal.freqs(b, a, np.logspace(1,2,500))
plt.semilogx(w,20* np.log10(abs(h)))# ⽤于绘制折线图，两个函数的 x 轴、y 轴分别是指数型的，并且转化为分贝
plt.title('Butterworth bandpass filter fit to constraints')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.fill([1,14,14,1],[-40,-40,99,99],'0.9', lw=0)# 区域颜⾊填充阻带
plt.fill([20,20,50,50],[-99,-3,-3,-99], facecolor='gray', alpha=0.2)# 通带，[x_左下, x_左上, x_右上, x_右下(左下⾓顺时针)], [y_左下, y_左上, y_右上, y_右下]
plt.fill([60,60,1e9,1e9],[99,-40,-40,99],'0.9', lw=0)# 阻带
指纹读取器plt.axis([10,100,-60,3])# 画布⼤⼩
plt.show()
（2）切⽐雪夫I型滤波器
切⽐雪夫I型数字和模拟滤波器设计函数：
(N, rp, Wn[, btype, analog, output, fs])
功能：设计⼀个N阶数字或模拟切⽐雪夫I型滤波器，然后返回滤波器系数。
调⽤格式：
b,a=scipy.signal.cheby1(N, rp, Wn, btype=‘low’, analog=False, output=‘ba’, fs=None)
z,p,k=scipy.signal.cheby1(N, rp, Wn, btype=‘low’, analog=False, output=‘zpa’, fs=None)
sos=scipy.signal.cheby1(N, rp, Wn, btype=‘low’, analog=False, output=‘sos’, fs=None)
参数：
发动机调速器N（int）（滤波器阶次）
rp（float）（通带内允许的最⼤纹波低于单位增益。以分贝指定，为正数。）
Wn（array、int）（临界频率对于低通和⾼通Wn是标量，对于带通和带阻Wn是长度为2的序列，对于
I型滤波器，这是过渡带中增益⾸先下降到-rp以下的点，对于数字滤波器，Wn与fs的单位相同，对于模拟滤波器，Wn是⾓频率rad / s）
btype：（‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’optional）（滤波类型，默认值是低通其中低通、⾼通、带通、带阻分别对应‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’默认值为低通）
analog：（bool, optional）（为True时返回模拟滤波器，False时返回数字滤波器）
output：（‘ba’, ‘zpk’, ‘sos’optional）（向后兼容的输出类型）
fs：（float）（数字系统采⽤频率）
返回值：
b（ndarray）（传递函数分⼦）
a（ndarray）（传递函数分母）
z（ndarray）（零点）
p（ndarray）（极点）
k（float）（系统增益）
sos（ndarray）（IIR滤波器⼆阶部分表⽰）
Chebyshev I型滤波器最⼤程度地提⾼了频率响应的通带和阻带之间的截⽌速率，但代价是通带中的纹波和阶跃响应中的振铃增加；I型滤波器⽐II型（cheby2）滚降的速度更快，但II型滤波器的通带中没有任何波动；等纹波通带具有N个最⼤值或最⼩值。
设计⼀个模拟滤波器并绘制其频率响应，显⽰关键点：
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
b, a = signal.cheby1(4,5,100,'low', analog=True)
w, h = signal.freqs(b, a)
玻璃瓶网
plt.semilogx(w,20* np.log10(abs(h)))
plt.title('Chebyshev Type I frequency response (rp=5)')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0,0.1)# 去除画布四周⽩边
plt.axvline(100, color='green')# # 绘制竖线，低通临界频率为100Hz
plt.axhline(-5, color='green')# 绘制横线-rp
plt.show()
⽣成由10 Hz和20 Hz组成的信号，以1 kHz采样
t = np.linspace(0,1,1000,False)
sig = np.sin(2*np.pi*10*t)+ np.sin(2*np.pi*20*t)
fig,(ax1, ax2)= plt.subplots(2,1, sharex=True)
ax1.plot(t, sig)
ax1.set_title('10 Hz and 20 Hz sinusoids')
ax1.axis([0,1,-2,2])
设计⼀个15 Hz的数字⾼通滤波器以消除10 Hz的⾳调，并将其应⽤于信号。（建议在过滤时使⽤⼆阶节格式，以避免传递函数（ba）格式出现数值错误）：
sos = signal.cheby1(10,1,15,'hp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 15 Hz high-pass filter')
ax2.axis([0,1,-2,2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()
切⽐雪夫I型滤波器阶次选择函数：
(wp, ws, gpass, gstop[, analog, fs])
调⽤格式：
N, Wn = scipy.signal.cheb1ord(wp, ws, gpass, gstop, analog=False, fs=None)
参数：
Wp（float）（通带边缘频率）
Ws（float）（阻带边缘频率）
对于数字滤波器，它们的单位与fs相同。默认情况下， fs是2个半周期/样本，因此将它们从0标准化为1，其中1是奈奎斯特频率。
（因此，wp和ws在半周期/样本中。）例如：
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议