2
$\begingroup$

I'm implementing digital Butterworth filter and encounter some numerical problem when filter order is high using direct form, so I wonder how to design the digital Butterworth filter and return its zeros, poles and gain directly, then I can convert zeros and poles to second-order sections. I tried to write a root-finding algorithm to find the zeros and poles from numerator and denominator coefficients of the digital filter, but a robust complex root-finding algorithm is really difficult to achieve. So I turn to design the zeros and poles of the filter.

What I have tried so far is calculate the poles of the lowpass analog Butterworth filter, then calculate the poles in $z$ domain using the bilinear transform: $$ s = c \frac{1-z^{-1}}{1+z^{-1}} = \frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}} \tag{1} $$ where $c = 2/T$ is the bilinear constant.

What confuses me is

  1. In matlab I use [z, p, k] = butter(n, wn) it gives me n zeros of $-1$ for lowpass filter and $1$ for highpass filter. It seems that the bilinear constant $c$ should equal to $1$ and then the analog zeros $s=\pm \infty$ can be mapped to $z=-1$ according to Eq. (1), is that correct?
  2. How to apply frequency transformation (lowpass to highpass, lowpass to bandpass, lowpass to bandstop) directly on the zeros and poles?
$\endgroup$
6
  • $\begingroup$ Matlab do have sos() -function to convert nth order filter to second order sections and also functions for frequency transformations (lp2hp(), lp2bp() and lp2bs()) which converts analog filter prototype or continuous-time state-space prototype to target filter (does plain poles and zeros define a filter type ?) . $\endgroup$
    – Juha P
    Jul 16 at 11:48
  • $\begingroup$ Is this a question about how to use the MATLAB butter() function? Or more fundamentally how to design a butterworth filter in the $z$-plane (digital domain)? $\endgroup$ Jul 16 at 19:59
  • $\begingroup$ @robertbristow-johnson It's about the second one: how to design a Butterworth filter and get its zeros, poles and gain, instead of the filter coefficients, and how to convert the zpk of the lowpass prototype to highpass, bandpass and bandstop variant. $\endgroup$
    – DSP novice
    Jul 17 at 11:43
  • $\begingroup$ @JuhaP Yes I know MATLAB's tf2sos converts transfer function to sos and I actually have tried to implement my own tf2sos in C but suffers numerical problem when finding roots of polynomial. So I want to design the digital filter through zeros, poles and gain directly. $\endgroup$
    – DSP novice
    Jul 17 at 11:52
  • 1
    $\begingroup$ @robertbristow-johnson Thank you Robert! Although my original question is directly getting zpk of digital Butterworth filter from its analog prototype using bilinear transform, your method seems can get each second-order section in $s$ domain and then apply BLT. I'll dive into it! $\endgroup$
    – DSP novice
    Jul 18 at 2:48

1 Answer 1

2
$\begingroup$

This can be designed in the following steps:

  1. Design a $N$-th order lowpass prototype of the analog Butterworth filter and you get an all-pole filter with poles $$ p_{a, n} = e^{j\pi(2n-1+N)/2N}, \ \ \ n=1,\ldots, N \tag{1} $$

  2. Apply bilinear transform to convert the zeros and poles to a digital lowpass filter. $$ s = \frac{2}{T} \frac{1-z^{-1}}{1+z^{-1}} = 2f_s\frac{1-z^{-1}}{1+z^{-1}} \tag{2} $$ Each pole is mapped to $$ p_{d, n} = \frac{2f_s+p_{a, n}}{2f_s-p_{a,n}} \tag{3} $$ and all zeros at $z_{a, n}=\pm\infty$ are mapped to $z_{d, n} = -1$. Then normalize the gain at DC, i.e., let $s=0$ and $z=1$ $$ k_a \frac{1}{\prod_{n=1}^N{(s-p_{a, n})}}\Bigg|_{s=0} = k_d \frac{\prod_{n=1}^N(z-z_{d, n})}{\prod_{n=1}^N(z-p_{d, n})} \Bigg|_{z=1} \tag{4} $$ and calculate $k_d$ from the equation above.

  3. Apply frequency transform to convert the lowpass prototype to lowpass, highpass, bandpass or bandstop. Take lowpass to lowpass for example, which follows the frequency transform: $$ \tilde{z}^{-1} \rightarrow g(z^{-1}) = \frac{z^{-1}-\alpha}{1-\alpha z^{-1}} \tag{5} $$ where $\alpha$ is a real frequency warping coefficient and $|\alpha|<1$. From Eq. (5) the relationship between the transformed zeros and poles and the original zeros and poles can be derived as $$ z = \frac{\tilde{z} + \alpha}{1+\alpha \tilde{z}}\tag{6} $$ Similar to Eq. (4) you need to normalize the magnitude at DC to get the gain of second-order sections. Of course you can combine this step with Eq. (4) and only normalize once with the analog prototype directly.

  4. Finally, convert the zeros, poles and gain to cascaded second-order sections.

$\endgroup$

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.