How to implement cascaded biquad transfer functions in Python/SciPy

Question

I am trying to implement in python/Scipy the biquad filter with the following S-domain transfer function: $$ H(s)=\frac{w_0^{2N-1}(s+w_z)}{[s^2+(w_0/Q)s+w_0^2]^N} $$

My problem here is that because of the N power denominator, I do not know how to write this transfer function into the M-order numerator and N-order denominator array-form that is required by the freqs or the filtfilt function in scipy.signal. As a work around, I decomposed the filter into a cascade of N-1 2nd-order lowpass and one 2nd-order bandpass section, such as (for N=4):

$$ H(s)=\frac{w_0^2}{s^2+(w_0/Q)s+w_0^2} \times \frac{w_0^2}{s^2+(w_0/Q)s+w_0^2} \times \frac{w_0^2}{s^2+(w_0/Q)s+w_0^2} \times \frac{w_0(s+w_z)}{s^2+(w_0/Q)s+w_0^2} $$ Here is my code:

from scipy import signal
import numpy as np

N=4
Q=10
fc=1000
worN=np.linspace(0, np.pi, int(fs/2))
wc =2*np.pi*fc/fs
w0=wc/(np.sqrt(((N-1)/(2*N-1))*(1-1/(2*Q[ii]**2))) * np.sqrt(1+np.sqrt(1+(1/((N-1)**2/(2*N-1)*(1-1/(2*Q[ii]**2))**2)))))

wz =1/10*w0

#LP transfers function
num_LP=[0, 0, w0**2]
den_LP=[1, w0/Q, w0**2]
w_LP, h_LP=signal.freqs(num_LP,den_LP,worN=worN)


#BP transfers function
num_BP=[0, w0, wz*w0]
den_BP=[1, w0/Q, w0**2]

w_BP, h_BP =signal.freqs(num_BP,den_BP,worN=worN)


# Cascad 3 LP biquads with 1 BP biquad filter
h_casc=h_LP**(N-1)*h_BP

This work around works well but it might not be the most efficient in terms of computation. In addition, my end goal is to optimize the filter's parameters using a Least-square algorithm and therefore the cascaded solution might not be the best one.

My question here is, how can I formulate this filter into a single transfer function so I can implement the actual filtering functions scipy.signal.filtfilt or scipy.signal.lfilter?

filtered_signal=scipy.signal.filtfilt(numerator,denumerator,input_signal)

Should I convert the transfer function to the Z-domain first?

Thank you so much for helping!

Answer 1

I do not know how to write this transfer function into the M-order numerator and N-order denominator form that is required by the filts function in scipy.signal.

What is the "filts" function?

#LP transfers function
num_LP=[0, 0, w0**2]
den_LP=[1, w0/Q, w0**2]

That looks right to me.

#BP transfers function
num_BP=[0, w0, wz**2]
den_BP=[1, w0/Q, w0**2]

Shouldn't that be w0*wz instead of wz**2?

# Cascad 3 LP biquads with 1 BP biquad filter
h_casc=h_LP**(N-1)*h_BP

This is fine for viewing the frequency response, but doesn't implement the actual filter.

My question here is, how can I formulate this filter into a single transfer function?

You shouldn't. You should leave it as four filter stages and apply them consecutively. Cascaded second-order sections are better than high-order transfer functions, which suffer more from numerical error.

Should I convert the transfer function to the Z-domain first?

Yes, your filter is analog, so you can view the analog frequency response with scipy.signal.freqs, but if you want to actually filter a signal, you need to convert it to a digital filter first. What transformation to use and what sample rate to use depends on your ultimate goal.

Answer 2

As endolith pointed out that the cascaded second-order sections are more numerical robust than direct form filters, if you insist to convert a series of cascaded filters into a direct form transfer function, just multiply the numerator and denominator polynomials respectively. $$ H(s) = \frac{B_1(s)}{A_1(s)} \cdot \frac{B_2(s)}{A_2(s)} = \frac{B_1(s)B_2(s)}{A_1(s)A_2(s)} $$ Multiplication of polynomials is equivalent to discrete convolution. See here. Thus the total transfer function is given by

num = conv(conv(conv(num_LP, num_LP), num_LP), num_BP)
den = conv(conv(conv(den_LP, den_LP), den_LP), den_BP)