Another bite of the BLT (Bilinear Transform)

I was a little bit premature in hoping to get into code the other day. One last dabble with some equations and we can march steadily along into some interesting code. Promise!

We've now taken a skim look at the general procedure of the BLT and an overview of the Butterworth low-pass filter design. In doing so we've used some good analytical maths and seen how the workings fall out. Very nice for blogging with the MathJax equations, but not quite so easy to map into a general purpose algorithm for code. So, let's have a look at that. It'll also explain one of the reasons for bashing away at the pole-zero representation.

Taking our analog representation in the s-plane, we can write it in terms of a product of complex roots of the polynomials of the numerator (top) and denominator (bottom), which we call the poles \( \rho\) on the bottom and zeros \( \sigma \) on the top.

$$ H_a(s) = K_a \frac {(s-\sigma_1)(s-\sigma_2)....(s-\sigma_m)}{(s-\rho_1)(s-\rho_2)....(s-\rho_m)}= K_a \frac{\prod^M_{m=1} {(s-\sigma_m)}}{\prod^N_{n=1} {(s-\rho_n)}} $$

Now, we know how to construct a filter design with poles. What if we could do a general purpose transformation of this into a similar representation in the z-plane for the digital implementation without needing to do lots of algebra. Let's have a go!

$$ H_d(z) = H_a(s) |_{s \rightarrow \frac {2}{T}\frac {(1-z^{-1})}{(1+z^{-1})}} = K_a \frac{\prod^M_{m=1} {\left(\frac {2}{T}\frac {(1-z^{-1})}{(1+z^{-1})} -\sigma_m\right)}}{\prod^N_{n=1} {\left(\frac {2}{T}\frac {(1-z^{-1})}{(1+z^{-1})} -\rho_n\right)}} $$

If we now multiply out the \( \frac {1}{T} \frac {1}{(1+z^{-1})} \) terms:

$$ H_d(z) = K_a \frac {T^N}{T^M} \frac {(1+z^{-1})^M}{(1+z^{-1})^N} \frac{\prod^M_{m=1} {\left( 2(1-z^{-1})-\sigma_m{(1+z^{-1})} \right)}}{\prod^N_{n=1} {\left( 2(1-z^{-1}) -\rho_n{(1+z^{-1})}\right)}} $$

Simplifying (a bit)

$$ H_d(z) = K_a T^{N-M}(1+z^{-1})^{M-N} \frac{\prod^M_{m=1} {\left( 2(1-z^{-1})-\sigma_m{(1+z^{-1})} \right)}}{\prod^N_{n=1} {\left( 2(1-z^{-1}) -\rho_n{(1+z^{-1})}\right)}} $$

Now, just taking the zeros (numerator) for the moment:

$$ \begin{align}
\left( 2(1-z^{-1})-\sigma_m{(1+z^{-1})} \right) & = (2-2z^{-1}-\sigma_mT-Tz^{-1}) \\
& = (-z^{-1}(2+\sigma_mT)+(2-\sigma_mT)) \\
& = (2+\sigma_mT)\left(z^{-1}-\frac {(2-\sigma_mT)}{(2+\sigma_mT)}\right)
\end{align}$$


Ad following the same procedure for the bottom:

$$ H_d(z) = K_a T^{N-M} \frac {\prod^M_{m=1}{(2+\sigma_mT)^M}}{\prod^N_{n=1}{(2+\rho_nT)^N}}  (1+z^{-1})^{M-N} \frac{\prod^M_{m=1} {\left(z^{-1}-\frac {(2-\sigma_mT)}{(2+\sigma_mT)}\right)}}{\prod^N_{n=1} {\left(z^{-1}-\frac {(2-\rho_nT)}{(2+\rho_nT)}\right)}} $$

What this huge monstrosity really says is that there is a large constant at the start, followed by some additional zeros matching the difference between the number of analogue complex poles and zeros followed by a set of poles and zeros transformed into the digital domain. What is really helpful is that this can then be made into an algorithm which performs some operations on the complex roots but does not need any algebraic manipulation. Ripe for code!

We'll need one additional function to go along with this to expand out the resulting factorised digital equation into a set of coefficients, but again, that should be a routine algorithm to progressively multiply out the factors.

On next to code...