python - 如何获得链式 IIR 滤波器的 b、a(分子/分母)？

Question

假设我们依次应用了 3 个过滤器:

b, a = iirfilter(...)  # or bilinear(...) or anything else producing b, a
y = lfilter(b, a, x)
b, a = iirfilter(...) 
y = lfilter(b, a, y)
b, a = iirfilter(...) 
y = lfilter(b, a, y)

如何获取系数b2 , a2相当于3个过滤器，这样我们只需通过一次lfilter就可以找到结果:

y = lfilter(b2, a2, x)

<强>？

<小时/>

编辑:卷积似乎不起作用:

fs = 44100
b2, a2 = iirfilter(2, 2.0/fs * np.asarray([40, 60]), btype='bandstop')  # 50 hz reject
b3, a3 = iirfilter(2, 2.0/fs * np.asarray([85, 115]), btype='bandstop')  # 100 hz reject
b = np.convolve(b2, b3)
a = np.convolve(a2, a3)
w, h = signal.freqz(b, a, worN=10000)

给出:

我尝试过 same , full , valid np.convolve 的参数，但都没有解决问题。

Answer 1

参见https://dsp.stackexchange.com/questions/38675/how-to-plot-magnitude-and-phase-response-of-2-cascaded-filters-in-matlab

您可以分别对分子和分母进行卷积

import scipy as sp
import scipy.signal as sig

# Individual filters
b1, a1 = sig.iirfilter(...)
b2, a2 = sig.iirfilter(...)

# Cascaded filter
a = sp.convolve(a1, a2)
b = sp.convolve(b1, b2)
y = sig.lfilter(b, a, x)

例如，您的采样率太高，并且复合滤波器的阶数不够长，无法对靠近的空值提供那么多的拒绝。降低采样率，然后内插至 44.1 kHz。

以下是采样率降低至 4410 Hz 时的结果。

fs = 4410.0
b2, a2 = sig.iirfilter(2, 2.0/fs * sp.asarray([40, 60]), btype='bandstop')  # 50 hz reject
w2, h2 = sig.freqz(b2, a2, worN=4096)

b3, a3 = sig.iirfilter(2, 2.0/fs * sp.asarray([85, 115]), btype='bandstop')  # 100 hz reject
w3, h3 = sig.freqz(b3, a3, worN=4096)

b = sp.convolve(b2, b3)
a = sp.convolve(a2, a3)
w, h = sig.freqz(b, a, worN=4096)

f = w/2.0*fs

然后将 IIR 滤波器的输出通过 10x 插值滤波器，以恢复到 44.1 kHz 采样率。

或者，减少过滤器阶数:

fs = 44100.0
b2, a2 = sig.iirfilter(1, 2.0/fs * sp.asarray([40, 60]), btype='bandstop')  # 50 hz reject
w2, h2 = sig.freqz(b2, a2, worN=4096)

b3, a3 = sig.iirfilter(1, 2.0/fs * sp.asarray([85, 115]), btype='bandstop')  # 100 hz reject
w3, h3 = sig.freqz(b3, a3, worN=4096)

b = sp.convolve(b2, b3, 'full')
a = sp.convolve(a2, a3, 'full')
w, h = sig.freqz(b, a, worN=4096)

其原始采样率为 44.1 kHz