Given two low-pass digital IIR filters, find bandpass coefficients

Question

In the case of FIR filters it is easy to get a band-pass filter by subtracting the coefficients of two low-pass filter filters or by convolving a high-pass and a low-pass filter:

import matplotlib.pyplot as plt
from scipy import signal
from math import pi

FS = 100.0          # Sampling rate
FC1 = 10/FS
N = 101             # Number of filter taps
a = 1               # Filter denominator
b1 = signal.firwin(N, cutoff=FC1, window='hamming') # Filter numerator
FC2 = 20/FS
b2 = signal.firwin(N, cutoff=FC2, window='hamming') # Filter numerator

h = b2-b1
w, H = signal.freqz(h)
_, ax = plt.subplots(1)
ax.plot(w/(2*pi), abs(H))

However, given two digital IIR filters with possibly different orders, I am not sure if it's possible to get the band-pass coefficients. Since the phase response of IIR filters is nonlinear in general, I can try to subtract the filtered outputs of the two low-pass filters, but this is rather crude. Since I haven't found a similar question on this site, so I am sorry if it doesn't make too much sense.

Answer 1

It's true that only impulse responses of linear phase (FIR) filters can be added or subtracted (if the delays are aligned) such that the resulting filter's magnitude response equals the sum or difference of the individual magnitude responses. This generally doesn't work for non-linear phase filters.

Obviously, cascading always works, because the magnitudes get multiplied and the phases are added, no matter if the filters have linear phase or not.

What can be done instead of adding or subtracting, regardless of the filter's phase response, is transform a given lowpass filter to a bandpass (or highpass, bandstop) filter using a frequency transformation. If $H_{LP}(z)$ is the transfer function of a lowpass filter, a bandpass filter (or any other filter type) can be obtained in the following way:

$$H_{BP}(z)=H_{LP}\big(A(z)\big)\tag{1}$$

where $A(z)$ is an allpass function:

$$A(z)=\frac{\displaystyle\sum_{n=0}^{N}\alpha_nz^{N-n}}{\displaystyle\sum_{n=0}^N\alpha_nz^{-n}},\qquad \alpha_0=1\tag{2}$$

For transformations from lowpass to highpass we use a first-order allpass function $A(z)$ (i.e., $N=1$ in $(2)$). For transformations to bandpass or bandstop filters, we need a second-order allpass function ($N=2$).

Details about how to obtain the desired allpass transformation functions $A(z)$ can be found in Section 7.4. of Discrete-Time Signal Processing by Oppenheim and Schafer (3rd edition). Mathworks also has a good documentation on digital frequency transformations.

Also take a look at this related answer.

Answer 2

If you just want to cascade a highpass with a lowpass you can simply cascade the poles and the zeros.

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

# create low pass and highpass of different order
fs = 48000     # sample rate
fhigh = 100    # low cutoff
flow = 2000    # high cutoff
# crate low pass and high pass
sos_high = signal.butter(3, fhigh * 2 / fs, btype='high', output='sos')
sos_low = signal.butter(2, flow * 2 / fs, btype='low', output='sos')
# concatenate the second order sectsions which concatenates the polse and zeros
sos = np.concatenate((sos_high, sos_low))
# plot it
fr0 = np.logspace(np.log10(20), np.log10(20000), 1000)
[w, h] = signal.sosfreqz(sos, fr0, fs=fs)
plt.clf()
plt.plot(fr0, 20 * np.log10(abs(h)))
plt.xscale("log")
plt.ylim(-60, 3)
plt.xlim(fr0[0], fr0[-1])
plt.grid("on")