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

关注

原创最新推荐文章于 2026-05-31 10:40:32 发布 1.9k 阅读 AI 写同款· GEO 优化 ›

在matlab里设计使用Bessel Filter（贝塞尔滤波器）进行数据降噪的简单例子作为一个土木工程系的研究生，现在大家做数据处理也用得越来越多了，自己接触的东西也多。作为土木工程师，很大部分时间其实并不想弄清楚滤波器的原理，只是希望拿来用，拿来把自己采集到的数据降噪一下，拿来做个傅立叶分析看看结构的固有频率变没变。但是matlab自带的文档很多时候都默认大家知道滤波器的原理，很多专有名词搞得不太友好，所以在这里以问题解答的形式直接给出我在科研中遇到的一个简单的用贝塞尔滤波器的例... 阅读详情

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）
"""

示例：

"""
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()

使用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)

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()

在这里插入图片描述

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.
"""

示例：

"""
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()

在这里插入图片描述

4、iirpeak²

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

接口：

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

示例：

#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()

在这里插入图片描述

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’
"""

示例：

#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()

在这里插入图片描述

引用：

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. ↩︎
Sophocles J. Orfanidis, “Introduction To Signal Processing”, Prentice-Hall, 1996 ↩︎ ↩︎ ↩︎ ↩︎