2
$\begingroup$

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!

$\endgroup$
1
  • 1
    $\begingroup$ So are you saying you want to see the cascaded response most simply, and then after you optimize that, you will then decompose it into 2nd order biquad sections? If that is the case, then just multiply out your denominator polynomial to have an Nth order polynomial and the coefficients of that will be your denominator vector. I am not confident that is really your question, could you clarify further what it is you want to do? $\endgroup$ Feb 22 at 23:22

2 Answers 2

2
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ Thank you very much for your answer. You are right, the filter is not implemented here and this is where I struggle. The functions to implement the filter are 'scipy.signal.filtfilt' or 'scipy.signal.lfilter'. They take as input the filter's numerator, the denumerator and the signal to be filtered. According to your answer I should implement each single second order stage separately, such as if N=4, the filtering function has to be implemented 4 times. Is that correct? You are correct regarding wz**2 (I'll edit that). $\endgroup$
    – papaya
    Feb 23 at 15:41
  • $\begingroup$ @papaya You will most likely want to convert your 2nd-order analog stages into 2nd-order digital stages, and then apply them to a sampled signal. What is the signal you are trying to filter? What is your ultimate goal? I'm not sure what you mean by N=4. You'd apply the LPF three times and then the BPF once. $\endgroup$
    – endolith
    Feb 23 at 16:01
  • $\begingroup$ Ok, thanks for confirming. N is the number of second order stages. Is a bilinear transformation suited for the conversion to the Z-domain (scipy.signal.bilinear)? The signal to filter will be gaussian noise/speech within 0-44kHz. My goal is to be able to automatically fit (though least square optimization) the cascaded filter frequency response to a given frequency response. $\endgroup$
    – papaya
    Feb 23 at 16:14
  • $\begingroup$ @papaya Bilinear should be fine if you're doing LPF and BPF. It always adds a zero at fs/2, so that is the region where it differs from the analog. But if you're applying LPFs it shouldn't make any noticeable difference. You don't mean that you want your upper cutoff to be 44 kHz, right? You mean that you have a signal that stops at 44.1 kHz, and you want to filter below that? $\endgroup$
    – endolith
    Feb 23 at 16:37
  • 1
    $\begingroup$ Sorry my bad. The signal will be within the range 20Hz-20kHz (audio) so probably using a sampling rate of 44 or 48kHz. Yes, I saw that the frequency response would be 0 at fs/2 but that's not a problem. Thanks again for your answers. $\endgroup$
    – papaya
    Feb 23 at 16:54
2
$\begingroup$

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)
$\endgroup$
1
  • 1
    $\begingroup$ Ok, that a very helpful point. I will keep the filter in cascaded second-order sections then. $\endgroup$
    – papaya
    Feb 23 at 15:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.