scipy Matlab-style IIR 滤波器设计下(Elliptic \Bessel \IIR-notch\IIR-peak\IIR-comb)

1、Elliptic

椭圆滤波器也被称为Cauer或Zo lotta rev滤波器，它最大限度地提高了频率响应的通带和阻带之间的过渡速率，但却牺牲了两者的纹波，并增加了阶跃响应中的振铃。
当rp接近0时，椭圆滤波器变成切比雪夫II型滤波器（切比2）。当rs接近0时，它变成切比雪夫I型滤波器（切比1）。当两个都接近０时，它变成巴特沃思滤波器（butter）。
等波纹通带具有N个最大值或最小值（例如，一个5阶滤波器有3个最大值和2个最小值）。因此，对于奇数阶滤波器，DC增益是单位，或者-rp dB 为偶数阶滤波器。

接口：

scipy.signal.ellip(N, rp, rs, Wn, btype='low', analog=False, output='ba', fs=None)
"""
rp:float
在通带的单位增益下最大允许纹波（dB）

rsfloat
在阻带的最小衰减（dB）
"""
1
2
3
4
5
6
7
8

示例：

"""
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.ellip(8, 1, 100, 17, 'hp', fs=1000, output='sos')
filtered = signal.sosfilt(sos, sig)
ax2.plot(t, filtered)
ax2.set_title('After 17 Hz high-pass filter')
ax2.axis([0, 1, -2, 2])
ax2.set_xlabel('Time [seconds]')
plt.tight_layout()
plt.show()
1
2
3
4
5
6
7
8
9
10
11

使用ellipord

"""
Design an analog highpass filter such that the passband is within 3 dB above 30 rad/s, while rejecting -60 dB at 10 rad/s. Plot its frequency response, showing the passband and stopband constraints in gray.
"""
N, Wn = signal.ellipord(30, 10, 3, 60, True)
b, a = signal.ellip(N, 3, 60, Wn, 'high', True)
1
2
3
4
5

2、Bessel/Thomson

也被称为汤姆逊滤波器，模拟贝塞尔滤波器具有最大平坦群延迟和最大线性相位响应，在阶跃响应中只有很少的振铃。¹
贝塞尔本质上是一个模拟滤波器。该函数使用双线性变换生成数字贝塞尔滤波器，该变换不保留模拟滤波器的相位响应。因此，它仅在低于大约fs/４频率下近似正确。为了在较高频率下获得最大平坦群延迟，必须使用保相技术对模拟贝塞尔滤波器进行变换。

查看群延迟：

b, a = signal.bessel(5, 1/0.1, 'low', analog=True, norm='delay')
w, h = signal.freqs(b, a)
plt.figure()
plt.semilogx(w[1:], -np.diff(np.unwrap(np.angle(h)))/np.diff(w))
plt.axhline(0.1, color='red')  # 0.1 seconds group delay
plt.title('Bessel filter group delay')
plt.xlabel('Frequency [radians / second]')
plt.ylabel('Group delay [seconds]')
plt.margins(0, 0.1)
plt.grid(which='both', axis='both')
plt.show()
1
2
3
4
5
6
7
8
9
10
11

3、iirnotch²

陷波滤波器是一种带宽窄（高品质因数）的带阻滤波器。它排斥一个狭窄的频带，其余的频谱几乎没有变化。

接口：

scipy.signal.iirnotch(w0, Q, fs=2.0)
"""
Parameters
w0:float
Frequency to remove from a signal. If fs is specified, this is in the same units as fs. By default, it is a normalized scalar that must satisfy 0 < w0 < 1, with w0 = 1 corresponding to half of the sampling frequency.从信号中去除的频率。如果指定了fs，则其单位与fs相同。默认情况下，它是必须满足0＜w0＜1的归一化标量，w0=1对应于采样频率的一半。

Q:float
Quality factor. Dimensionless parameter that characterizes notch filter -3 dB bandwidth bw relative to its center frequency, Q = w0/bw.品质因数。表征陷波器的无量纲参数-3dB相对于其中心频率的带宽bw，Q＝w0/bw

fs:float, optional
The sampling frequency of the digital system.数字系统的采样频率。


Returns
b, a:ndarray, ndarray
Numerator (b) and denominator (a) polynomials of the IIR filter.
"""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

示例：

"""
Design and plot filter to remove the 60 Hz component from a signal sampled at 200 Hz, using a quality factor Q = 30
"""
from scipy import signal
import matplotlib.pyplot as plt
fs = 200.0  # Sample frequency (Hz)
f0 = 60.0  # Frequency to be removed from signal (Hz)
Q = 30.0  # Quality factor
# Design notch filter
b, a = signal.iirnotch(f0, Q, fs)
# Frequency response
freq, h = signal.freqz(b, a, fs=fs)
# Plot
fig, ax = plt.subplots(2, 1, figsize=(8, 6))
ax[0].plot(freq, 20*np.log10(abs(h)), color='blue')
ax[0].set_title("Frequency Response")
ax[0].set_ylabel("Amplitude (dB)", color='blue')
ax[0].set_xlim([0, 100])
ax[0].set_ylim([-25, 10])
ax[0].grid(True)
ax[1].plot(freq, np.unwrap(np.angle(h))*180/np.pi, color='green')
ax[1].set_ylabel("Angle (degrees)", color='green')
ax[1].set_xlabel("Frequency (Hz)")
ax[1].set_xlim([0, 100])
ax[1].set_yticks([-90, -60, -30, 0, 30, 60, 90])
ax[1].set_ylim([-90, 90])
ax[1].grid(True)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

4、iirpeak²

设计二阶IIR峰值(谐振)数字滤波器。峰值滤波器是带宽窄(高品质因子)的带通滤波器.它在窄带之外拒绝组件。

接口：

scipy.signal.iirpeak(w0, Q, fs=2.0)# 含义和iirnotch一样
1

示例：

#Design and plot filter to remove the frequencies other than the 300 Hz component from a signal #sampled at 1000 Hz, using a quality factor Q = 30

from scipy import signal
import matplotlib.pyplot as plt
fs = 1000.0  # Sample frequency (Hz)
f0 = 300.0  # Frequency to be retained (Hz)
Q = 30.0  # Quality factor
# Design peak filter
b, a = signal.iirpeak(f0, Q, fs)
# Frequency response
freq, h = signal.freqz(b, a, fs=fs)
# Plot
fig, ax = plt.subplots(2, 1, figsize=(8, 6))
ax[0].plot(freq, 20*np.log10(np.maximum(abs(h), 1e-5)), color='blue')
ax[0].set_title("Frequency Response")
ax[0].set_ylabel("Amplitude (dB)", color='blue')
ax[0].set_xlim([0, 500])
ax[0].set_ylim([-50, 10])
ax[0].grid(True)
ax[1].plot(freq, np.unwrap(np.angle(h))*180/np.pi, color='green')
ax[1].set_ylabel("Angle (degrees)", color='green')
ax[1].set_xlabel("Frequency (Hz)")
ax[1].set_xlim([0, 500])
ax[1].set_yticks([-90, -60, -30, 0, 30, 60, 90])
ax[1].set_ylim([-90, 90])
ax[1].grid(True)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

5、iircomb²

陷波梳状滤波器由具有窄带宽（高品质因数）的规则间隔带阻滤波器组成。每一个都排斥一个狭窄的频带，使其余的频谱几乎没有变化。
峰值梳状滤波器由具有窄带宽（高品质因数）的规则间隔带通滤波器组成。每个都拒绝一个狭窄的频带之外的组件。

实现的详细信息，请参阅². 梳状滤波器的TF实现是数值稳定的，即使在高阶由于使用一个单一的重复极点，这将不会遭受精度损失。

接口：

scipy.signal.iircomb(w0, Q, ftype='notch', fs=2.0, *, pass_zero=False)
"""
Parameters

f:type{‘notch’, ‘peak’}
The type of comb filter generated by the function. If ‘notch’, then the Q factor applies to the notches. If ‘peak’, then the Q factor applies to the peaks. Default is ‘notch’.
函数生成的梳状滤波器的类型。如果是“凹口”，则Q因子应用于凹口。如果“峰值”，则Q因子应用于峰值。默认值为“槽口”。

pass_zero:bool, optional
If False (default), the notches (nulls) of the filter are centered on frequencies [0, w0, 2*w0, …], and the peaks are centered on the midpoints [w0/2, 3*w0/2, 5*w0/2, …]. If True, the peaks are centered on [0, w0, 2*w0, …] (passing zero frequency) and vice versa.
如果为False（默认值），则滤波器的凹口（零值）以频率【0,w0,2*w0，...】为中点，而峰值则以中点【w0/2,3*w2,5*wO/2，…】为中点。如果为真，则峰值以【0，w0，2*w0,...】为中点（通过零频率），反之亦然。

Returns
b, andarray, ndarray
Numerator (b) and denominator (a) polynomials of the IIR filter.

Raises
ValueError
If w0 is less than or equal to 0 or greater than or equal to fs/2, if fs is not divisible by w0, if ftype is not ‘notch’ or ‘peak’
"""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

示例：

#Design and plot notching comb filter at 20 Hz for a signal sampled at 200 Hz, using quality factor #Q = 30

from scipy import signal
import matplotlib.pyplot as plt
fs = 200.0  # Sample frequency (Hz)
f0 = 20.0  # Frequency to be removed from signal (Hz)
Q = 30.0  # Quality factor
# Design notching comb filter
b, a = signal.iircomb(f0, Q, ftype='notch', fs=fs)
# Frequency response
freq, h = signal.freqz(b, a, fs=fs)
response = abs(h)
# To avoid divide by zero when graphing
response[response == 0] = 1e-20
# Plot
fig, ax = plt.subplots(2, 1, figsize=(8, 6), sharex=True)
ax[0].plot(freq, 20*np.log10(abs(response)), color='blue')
ax[0].set_title("Frequency Response")
ax[0].set_ylabel("Amplitude (dB)", color='blue')
ax[0].set_xlim([0, 100])
ax[0].set_ylim([-30, 10])
ax[0].grid(True)
ax[1].plot(freq, (np.angle(h)*180/np.pi+180)%360 - 180, color='green')
ax[1].set_ylabel("Angle (degrees)", color='green')
ax[1].set_xlabel("Frequency (Hz)")
ax[1].set_xlim([0, 100])
ax[1].set_yticks([-90, -60, -30, 0, 30, 60, 90])
ax[1].set_ylim([-90, 90])
ax[1].grid(True)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

引用：

---

1. Thomson, W.E., “Delay Networks having Maximally Flat Frequency Characteristics”, Proceedings of the Institution of Electrical Engineers, Part III, November 1949, Vol. 96, No. 44, pp. 487-490. ↩︎

2. Sophocles J. Orfanidis, “Introduction To Signal Processing”, Prentice-Hall, 1996 ↩︎ ↩︎ ↩︎ ↩︎