Как получить b, (числитель / знаменатель) для цепных БИХ-фильтров?

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 кГц.

Вот результаты с частотой дискретизации, уменьшенной до 4410 Гц.

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

Затем пропустите выход БИХ-фильтра через 10-кратный интерполяционный фильтр, чтобы вернуться к частоте выборки 44,1 кГц.

ИЛИ, уменьшите порядок фильтра:

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 кГц