scipy Matlab-style IIR 滤波器设计上(Butterworth\Chebyshev type I \Chebyshev type II )

scipy Matlab-style IIR 滤波器设计上(Butterworth\Chebyshev type I \Chebyshev type II )

各种滤波接口

滤波器接口	含义
butter(N, Wn[, btype, analog, output, fs])	设计Butterworth模拟和数字滤波器
buttord(wp, ws, gpass, gstop[, analog, fs])	自动选择满足条件的最低阶Butterworth
cheby1(N, rp, Wn[, btype, analog, output, fs])	设计Chebyshev type I 模拟和数字滤波器
cheb1ord(wp, ws, gpass, gstop[, analog, fs])	自动选择满足条件的最低阶Chebyshev type I
cheby2(N, rs, Wn[, btype, analog, output, fs])	设计Chebyshev type II 模拟和数字滤波器
cheb2ord(wp, ws, gpass, gstop[, analog, fs])	自动选择满足条件的最低阶Chebyshev type II
ellip(N, rp, rs, Wn[, btype, analog, output, fs])	设计Elliptic (Cauer) 模拟和数字滤波器
ellipord(wp, ws, gpass, gstop[, analog, fs])	自动选择满足条件的最低阶Elliptic (Cauer)
bessel(N, Wn[, btype, analog, output, norm, fs])	设计Bessel/Thomson 模拟和数字滤波器
iirnotch(w0, Q[, fs])	Design second-order IIR notch digital filter.设计二阶IIR陷波数字滤波器
iirpeak(w0, Q[, fs])	Design second-order IIR peak (resonant) digital filter.设计二阶IIR峰值（谐振）数字滤波器。
iircomb(w0, Q[, ftype, fs, pass_zero])	Design IIR notching or peaking digital comb filter.设计IIR陷波或调峰数字梳状滤波器

1、Butterworth

巴特沃斯滤波器在通带中具有最大平坦（ flat frequency）的频率响应。

Butterworth接口参数说明：

'''
scipy.signal.butter(N, Wn, btype='low', analog=False, output='ba', fs=None)[source]

N:int
滤波器阶数，对于“带通”和“带阻”滤波器，最后二阶部分的结果（“sos”）矩阵是2*N，N是所需系统的二次型部分的数目。

Wn:array_like
截止频率，对于低通和高通滤波器，Wn是标量；对于带通和带阻滤波器来说，Wn是长度为2的序列。
对于Butterworth滤波器，这是增益下降到通带的1/sqrt（2）的点（“-3 dB点”）。
对于数字滤波器，如果不指定fs，Wn单位从0归一化为1，其中1是奈奎斯特频率（Wn定义为2*截止频率/fs）。如果指定了fs，则Wn的单位与fs相同HZ。（fs是采样频率）：例如低通滤波器，数据的采样频率fs=25HZ，截止频率10HZ.有两种方式，第一种是Wn=10, fs=25.另一种是Wn=2*10/25,不填fs参数。
对于模拟滤波器，Wn是角频率（如rad/s）。

btype{‘lowpass’, ‘highpass’, ‘bandpass’, ‘bandstop’}
可选的滤波器 Default is ‘lowpass’.

analog:bool, optional
True, 得到模拟滤波器, 否则数字滤波器。

output{‘ba’, ‘zpk’, ‘sos’}, optional
Type of output: 分子/分母 (‘ba’), 极点-零点 (‘zpk’), 或者 二阶型 (‘sos’). 默认  ‘ba’ 用于后向兼容，建议使用'sos'.因为如果要求传递函数形式[b,a]，可能会发生数值问题，因为根和多项式系数之间的转换是一个数值敏感的操作。当N>=4，建议使用'sos'。

fs:float, optional
数字系统的采样频率

Returns（函数返回，和选择的output类型有关）
b, a：ndarray, ndarray
Numerator (b) and denominator (a) polynomials of the IIR filter. Only returned if output='ba'.

z, p, kndarray, ndarray, float
Zeros, poles, and system gain of the IIR filter transfer function. Only returned if output='zpk'.

sosndarray
Second-order sections representation of the IIR filter. Only returned if output=='sos'.
'''

Butterworth示例：

# @Version       :scipy 1.9.3 numpy 1.23.4 matplotlib3.6.2

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np


# 设计模拟滤波器，查看幅频响应
b, a = signal.butter(N=4, Wn=100 ,analog=True)
w, h = signal.freqs(b, a)
plt.semilogx(w, 20 * np.log10(abs(h)))  # x轴转换为对数格式的方式显示数据
plt.title('Butterworth filter frequency response')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.axvline(100, color='green') # cutoff frequency
plt.show()

# 生成 10 Hz and 20 Hz信号, 采样率1 kHz
t = np.linspace(0, 1, 1000, False)  # 1 second
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以下信号,使用'sos'避免'ba'的数值问题。
sos = signal.butter(10, 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()

使用自动最小阶buttord：

接口说明：

scipy.signal.buttord(wp, ws, gpass, gstop, analog=False, fs=None)[source]
"""
Return the order of the lowest order digital or analog Butterworth filter that loses no more than gpass dB in the passband and has at least gstop dB attenuation in the stopband.
得到在通带中损失不超过gpass dB,并且在阻带中具有至少gstop dB衰减损失的最低阶数字或模拟巴特沃斯滤波器的阶数;

Parameters
wp, ws:float
Passband and stopband edge frequencies.通带和阻带边缘频率。

For digital filters, these are in the same units as fs. By default, fs is 2 half-cycles/sample, so these are normalized from 0 to 1, where 1 is the Nyquist frequency. (wp and ws are thus in half-cycles / sample.) For example:(参考前面对Butterworth的注解)

Lowpass: wp = 0.2, ws = 0.3
Highpass: wp = 0.3, ws = 0.2
Bandpass: wp = [0.2, 0.5], ws = [0.1, 0.6]
Bandstop: wp = [0.1, 0.6], ws = [0.2, 0.5]

For analog filters, wp and ws are angular frequencies (e.g., rad/s).

gpass:float
The maximum loss in the passband (dB).

gstop:float
The minimum attenuation衰减 in the stopband (dB).

analog:bool, optional
When True, return an analog filter, otherwise a digital filter is returned.

fs:float, optional（决定wp,ws的单位HZ）
The sampling frequency of the digital system.

Returns
ord:int
The lowest order for a Butterworth filter which meets specs.

wn:ndarray or float
The Butterworth natural frequency (i.e. the “3dB frequency”). Should be used with butter to give filter results. If fs is specified, this is in the same units, and fs must also be passed to butter.巴特沃斯自然频率（即"3dB频率"）。应与butter一起使用得到滤波器结果。如果指定了fs，就是HZ单位，fs也必须传递到butter。（相当于是指定fs单位就是HZ）
"""

示例：

"""
Design an analog bandpass filter with passband within 3 dB from 20 to 50 rad/s, while rejecting at least -40 dB below 14 and above 60 rad/s. Plot its frequency response, showing the passband and stopband constraints in gray.
"""

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np


N, Wn = signal.buttord([20, 50], [14, 60], 3, 40, True)# 没有指定fs，因此Wn指的是2 half-cycles/sample
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)))
plt.title('Butterworth bandpass filter fit to constraints')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')
plt.fill([1,  14,  14,   1], [-40, -40, 99, 99], '0.9', lw=0) # stop
plt.fill([20, 20,  50,  50], [-99, -3, -3, -99], '0.9', lw=0) # pass
plt.fill([60, 60, 1e9, 1e9], [99, -40, -40, 99], '0.9', lw=0) # stop
plt.axis([10, 100, -60, 3])
plt.show()

2、Chebyshev type I

The Chebyshev type I filter maximizes the rate of cutoff between the frequency response’s passband and stopband, at the expense of ripple in the passband and increased ringing in the step response.
Type I filters roll off faster than Type II (cheby2), but Type II filters do not have any ripple in the passband.
The equiripple passband has N maxima or minima (for example, a 5th-order filter has 3 maxima and 2 minima). Consequently, the DC gain is unity for odd-order filters, or -rp dB for even-order filters.
切比雪夫I型滤波器的频率响应的通带和阻带之间的截止速率最大化，在通带的纹波ripple 和阶跃响应的振铃ringing增加。
类型I型滤波器比II型（cheby2）的滤波器滚转速度 roll off 更快，但II类滤波器在通带中没有任何波纹。
等纹波通带具有N个最大值或最小值（例如，一个5阶滤波器有3个最大值和2个最小值）。因此，对于奇数 odd-order阶滤波器，直流DC增益是1 unity，或者偶数even-order 阶滤波器-rp dB 。

接口说明：

scipy.signal.cheby1(N, rp, Wn, btype='low', analog=False, output='ba', fs=None)
'''
rp：float
The maximum ripple allowed below unity gain in the passband. Specified in decibels, as a positive number.在通带中允许低于单位增益的最大纹波。用正数,单位分贝dB表示的。其他参数和Butterworth一样
'''

示例：

sos = signal.cheby1(10, 1, 15, 'hp', fs=1000, output='sos')

使用cheb1ord设计auto:

scipy.signal.cheb1ord(wp, ws, gpass, gstop, analog=False, fs=None)
"""
参数和buttord一样
"""

示例：

"""
Design a digital lowpass filter such that the passband is within 3 dB up to 0.2*(fs/2), while rejecting at least -40 dB above 0.3*(fs/2). Plot its frequency response, showing the passband and stopband constraints in gray.相当于截止频率0.3*(fs/2)HZ，转换成0.2rad/s和0.3rad/s，就不用fs参数，不然就需要
"""

from scipy import signal
import matplotlib.pyplot as plt
N, Wn = signal.cheb1ord(0.2, 0.3, 3, 40)
b, a = signal.cheby1(N, 3, Wn, 'low')
w, h = signal.freqz(b, a)
plt.semilogx(w / np.pi, 20 * np.log10(abs(h)))
plt.title('Chebyshev I lowpass filter fit to constraints')
plt.xlabel('Normalized frequency')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')
plt.fill([.01, 0.2, 0.2, .01], [-3, -3, -99, -99], '0.9', lw=0) # stop
plt.fill([0.3, 0.3,   2,   2], [ 9, -40, -40,  9], '0.9', lw=0) # pass
plt.axis([0.08, 1, -60, 3])
plt.show()

3、Chebyshev type II

切比雪夫II型滤波器的频率响应的通带和阻带之间的截止速率最大化，在通带的纹波ripple 和阶跃响应的振铃ringing增加。
II型过滤器不会像I型（cheby1）那样快速地滚落 roll of。

接口：

scipy.signal.cheby2(N, rs, Wn, btype='low', analog=False, output='ba', fs=None)

"""
rs:float
The minimum attenuation required in the stop band. Specified in decibels, as a positive number.在阻带中所需的最小衰减。正数dB。
"""

示例：

"""
Design a digital high-pass filter at 17 Hz to remove the 10 Hz tone, and apply it to the signal. (It’s recommended to use second-order sections format when filtering, to avoid numerical error with transfer function (ba) format):
"""
sos = signal.cheby2(12, 20, 17, 'hp', fs=1000, output='sos')

auto设计cheby2:

"""
Design a digital bandstop filter which rejects -60 dB from 0.2*(fs/2) to 0.5*(fs/2), while staying within 3 dB below 0.1*(fs/2) or above 0.6*(fs/2). Plot its frequency response, showing the passband and stopband constraints in gray.
"""

from scipy import signal
import matplotlib.pyplot as plt
N, Wn = signal.cheb2ord([0.1, 0.6], [0.2, 0.5], 3, 60)
b, a = signal.cheby2(N, 60, Wn, 'stop')
w, h = signal.freqz(b, a)
plt.semilogx(w / np.pi, 20 * np.log10(abs(h)))
plt.title('Chebyshev II bandstop filter fit to constraints')
plt.xlabel('Normalized frequency')
plt.ylabel('Amplitude [dB]')
plt.grid(which='both', axis='both')
plt.fill([.01, .1, .1, .01], [-3,  -3, -99, -99], '0.9', lw=0) # stop
plt.fill([.2,  .2, .5,  .5], [ 9, -60, -60,   9], '0.9', lw=0) # pass
plt.fill([.6,  .6,  2,   2], [-99, -3,  -3, -99], '0.9', lw=0) # stop
plt.axis([0.06, 1, -80, 3])
plt.show()