2
$\begingroup$

I am currently experimenting with a digital implementation of a Chebyshev Type 1 filter of mine, for audio purposes. I have written the code myself so I was wondering if the behavior I've noticed is due to bad code or something else. My implementation works perfectly fine up to an order of 30. I know that this is ridiculous but that's not the point, it's just for educational purposes. Anyway, when using a filter of order 30 or higher part of the frequency response that should be flat has resonant peaks. For a lowpass filter, it meant boosted low frequencies for a highpass filter it resulted in boosted low-mids. The highpass filter did not behave this way when the cutoff was above the resonant peaks, by the way, and changing the ripple didn't change the trend. Is this normal and a result of the filter getting more and more unstable? Or does this indicate that there is some part of the code not working properly? I am using double precision where I can, but the framework I am working in wants the sample values with 32-bit precision so I can only do the filter coefficients and filter memory with 64-bit precision and have to reduce precision when writing to the audio buffer. Maybe that's an issue?

Edit: To add to my OP, I am using cascaded biquads to realize the filter.

$\endgroup$
  • $\begingroup$ How do you design the filters? In transfer function form or in zero/pole/gain form? Matlab's standard [b,a] = cheby1(10,.1,.2); sos=tf2sos(b,a) doesn't work for higher orders. Do [z,p,k] = cheby1(10,.1,.2); sos = zpk2sos(z,p,k); instead. $\endgroup$ – Hilmar Oct 4 '19 at 12:52
  • $\begingroup$ I am not using Matlab but a framework for audio plugins. I've used the formulas for the pole locations to get my transfer function and then used the bilinear transform to get the digital filter. $\endgroup$ – Agiltohr Oct 5 '19 at 18:38
2
$\begingroup$

With infinite-precision arithmetic they will be stable. However, since you don't have infinite-precision arithmetic, you will have quantization issues even if you use 64-bit precision. These quantization issues can make your filter unstable. Even if your filter is stable, perhaps you will not get the frequency response wanted because of these quantization issues.

First solution : split your order-30 IIR filter into 15 order-2 IIR filters and cascade them. This should fix your problem.

https://www.dsprelated.com/freebooks/filters/Series_Second_Order_Sections.html https://www.mathworks.com/help/signal/ref/tf2sos.html

Alternatively : Instead of cascading, you can parallelize the order-2 IIR filters.

https://www.dsprelated.com/freebooks/filters/Parallel_First_and_or_Second_Order.html

Finally, do you really need an order-30 chebyshev filter? I find it hard to believe to you really more than an order-8 or order-10 IIR filter.

Edit : You edited your post to mention that you already use biquads. Have you tried an simpler filter, like an order-4 IIR? Make sure it works for an order-4 IIR filter (i.e. 2 cascaded biquads) and then try increasing the order. Also, are you sure that your factorization algorithm is correct? Maybe, there's a simple bug hidden somewhere and it will be easier to spot it with a reduced-order filter.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for your input. I forgot to add, that I already am using cascaded biquads, so limited precision is out of the question then? And you're right, I don't need a 30 order filter but I am experimenting because it's fun and I come across issue like these that maybe teach me a thing or two. $\endgroup$ – Agiltohr Oct 3 '19 at 22:46
2
$\begingroup$

For the biquad section that is cascaded, the quantization issues regarding the pole locations are well understood. For a biquad transfer function:

$$\begin{align} H(z) &= \frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}} \\ \\ &= \frac{b_0z^2+b_1z+b_2}{z^2+a_1z+a_2} \\ \\ &= b_0\frac{z^2+\frac{b_1}{b_0}z+\frac{b_2}{b_0}}{z^2+a_1z+a_2} \\ \\ &= b_0\frac{(z-q_1)(z-q_2)}{(z-p_1)(z-p_2)} \\ \\ \end{align}$$

The zeros, $q_1,q_2$ and the poles $p_1,p_2$ are roots to the quadratic equations of $z$ in the numerator and denominator, respectively. The stability of the biquad is determines solely as a function of the poles:

$$ |p_k|<1 $$

The closer the pole is to 1, the less stable the biquad section is.

These to pole locations are:

$$ p_1,p_2 = -\frac{a_1}{2} \pm \sqrt{\frac{a_1^2}{4} - a_2} $$

and if $a_2>\frac{a_1^2}{4}$ the poles are complex conjugate and we factor a $\sqrt{-1}$ out of the radical:

$$ p_1,p_2 = -\frac{a_1}{2} \pm j\sqrt{a_2 - \frac{a_1^2}{4}} $$

For the latter case, it's easy:

$$ |p_k| = \sqrt{a_2} $$

$$ \operatorname{Re} \{ p_k \} = -\frac{a_1}{2} $$

You can right away see that the range of the two denominator coefficients are

$$ 0 < a_2 < 1$$ and $$ -2 < a_1 < +2 $$

to keep these two poles inside the unit circle.

... more to come ...

| improve this answer | |
$\endgroup$
  • $\begingroup$ I see. So you're saying I should just check the coefficients? Because the filter is still stable it just does not exhibit the exact expected frequency response. $\endgroup$ – Agiltohr Oct 4 '19 at 8:06
  • 1
    $\begingroup$ i am not done with this yet, but we need to discuss now what the effects on the quantization of coefficients have on pole location. you can see that the distance the two poles are from the origin (which is 1 minus the distance from the unit circle) depends solely on $a_2$ and the left-right placement of the two poles (the real part) depends solely on $a_1$. it's $a_1$ and $a_2$ that are quantized to fixed-point or floating-point values. $\endgroup$ – robert bristow-johnson Oct 4 '19 at 11:43
  • 1
    $\begingroup$ if you have 32-bit floats, your pole and zero placement should be just fine. but if you have 16-bit fixed, there are serious quantization issues. i might not be able to find time to finish this answer for a while. i'm sorry but several weekend planning issues are consuming my attention now. $\endgroup$ – robert bristow-johnson Oct 4 '19 at 11:48
  • 1
    $\begingroup$ whether it's a Tchebyshev (Type 1 or 2) or Butterworth or Elliptical or whatever prototype, the higher the order, the closer to the unit circle the most dominant poles will be. but as long as you round $a_2$ down when you quantize that coefficient, there is no reason for any biquad section to go unstable. $\endgroup$ – robert bristow-johnson Oct 4 '19 at 11:54
  • $\begingroup$ Thanks Robert! Looks like I'll have to go on a bug hunt. $\endgroup$ – Agiltohr Oct 4 '19 at 11:55

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.