# Cascade Filter (recursive function) - Help

Hi All

I am trying to create a cascaded BiQuad filter using nyquist. I have created this code. I have also provided (as an example) the coefficients required for a 10th order low-pass filter with Fc = 1412.54 Hz at 44100 Hz sampling rate… (FYI this is the upper cut-off frequency of the 1kHz octave band). But my question is, why doesn’t my code work…

``````(defun cascade (s1 i n coeffs)
;s1 = sound to be filtered
;i  = current filter (0 to (n/2 - 1))
;n  = order of filter (integer multiple of 2)
;coeffs = (a0_0 a1_0 a2_0 b0_0 b1_0 b2_0 a0_1 ...)
;this filter will call itself so that the filters are
;cascaded, each step will implement the next set of
;biquad coefficients until i = 0 when this function will
;stop recursion and return only the last biquaded sound.

(setf a0 (aref coeffs (+ (* 6 i) 0)))
(setf a1 (aref coeffs (+ (* 6 i) 1)))
(setf a2 (aref coeffs (+ (* 6 i) 2)))
(setf b0 (aref coeffs (+ (* 6 i) 3)))
(setf b1 (aref coeffs (+ (* 6 i) 4)))
(setf b2 (aref coeffs (+ (* 6 i) 5)))

(if (<= i 0)
(biquad s1 b0 b1 b2 a0 a1 a2)
(biquad (cascade s1 (- i 1) n coeffs) b0 b1 b2 a0 a1 a2)
)
)
``````

If I run this like this

``````(setf coeffs (vector 1.000000 1.636525 -0.670236 0.008428 0.016855 0.008428 1.000000 1.663372 -0.697635 0.008566 0.017132 0.008566 1.000000 1.716946 -0.752313 0.008842 0.017683 0.008842 1.000000 1.796591 -0.833599 0.009252 0.018504 0.009252 1.000000 1.900213 -0.939355 0.009786 0.019571 0.009786))

``````

I should get a 10th order low-pass filter with fc=1412.54 Hz. But I don’t I only get something that looks about 4th order and the cutoff frequency is all wrong (presumably because the whole set of filters have not worked).

Compare the “plot spectrum…” frequency response of some white noise from the above to the result from the following

``````(biquad (biquad (biquad (biquad (biquad s 0.008428 0.016855 0.008428 1.000000 1.636525 -0.670236) 0.008566 0.017132 0.008566 1.000000 1.663372 -0.697635) 0.008842 0.017683 0.008842 1.000000 1.716946 -0.752313) 0.009252 0.018504 0.009252 1.000000 1.796591 -0.833599) 0.009786 0.019571 0.009786 1.000000 1.900213 -0.939355)
``````

This second bit of code is what I am trying to achieve, and as far as I can tell my cascade filter should output exactly this. The coefficients are identical (and I know these coefficients are correct). But my cascade function only seems to spit out about half of the filtering it should… is this a syntax problem? a problem with recursion or am I doing something strange?

Any help would be appreciated.

There’s two problems.
The first is a minor coding error - you’ve forgotten to set s1 to the evaluated value on each iteration.
After running:

``````  (if (<= i 0)
(biquad s1 b0 b1 b2 a0 a1 a2)
(biquad (cascade s1 (- i 1) n coeffs) b0 b1 b2 a0 a1 a2)
``````

“s1” is unchanged, whereas:

`````` (setf s1 (biquad s1 b0 b1 b2 a0 a1 a2))
``````

would appear to evaluate the biquad function and set “s1” to the result.

The second is a bit less obvious, but an important feature of Nyquist, which is lazy evaluation.
Even with the above change it will not evaluate as expected but will only evaluate the biquad function on the last iteration.
Try iterating with a “do” loop instead, for example:

``````(defun cascade (s1 times coeffs)
(do ((i times (setq i (1- i))))
((< i 0) s1)
(let ((a0 (aref coeffs (+ (* 6 i) 0)))
(a1 (aref coeffs (+ (* 6 i) 1)))
(a2 (aref coeffs (+ (* 6 i) 2)))
(b0 (aref coeffs (+ (* 6 i) 3)))
(b1 (aref coeffs (+ (* 6 i) 4)))
(b2 (aref coeffs (+ (* 6 i) 5))))
(format t "Coefficients for i = ~a:~%a0 ~a   a1 ~a   a2 ~a~%b0 ~a   b1 ~a   b2 ~a~%~%"
i a0 a1 a2 b0 b1 b2); print to debug window
(setf s1 (biquad s1 b0 b1 b2 a0 a1 a2)))))

;; Set coefficients
(setf coeffs
(vector 1.000000 1.636525 -0.670236 0.008428 0.016855 0.008428
1.000000 1.663372 -0.697635 0.008566 0.017132 0.008566
1.000000 1.716946 -0.752313 0.008842 0.017683 0.008842
1.000000 1.796591 -0.833599 0.009252 0.018504 0.009252
1.000000 1.900213 -0.939355 0.009786 0.019571 0.009786))
(setq *float-format* "%f") ; decimal notation
``````

Thanks very much - worked a treat.

Don’t know why I didn’t just think to use a loop in the first place… recursive must have been on the brain - all those recursive filter coefficients.

I’ve run into the same sort of problem myself and spent hours and hours trying to work out why it didn’t work. Out of interest, you said that you knew the coefficients were correct - I was wondering how you worked them out.

I have coded up an nth order butterworth filter which outputs the biquad coefficients needed to go into a cascade to get the desired cut-off frequency. Then, together with a LP and HP version of this you can achieve a flat-top band pass filter. I am also working on getting chebyshev type I filters working and that is almost done. when its finished I’ll post the code.

I am an acoustic consultant and more and more frequently have to analyse recorded wav files (recorded on our level metre). Its very handy to be able to do this in audacity directly. So I have been trying to get ISO standard compliant octave band and 1/3 octave band filters coded. 99% done now.

FYI its based on the code in chapter 20 of the free online book from here http://www.dspguide.com/.

Excellent. Do you also do this sort of thing “for fun”? If so I may have questions for you in the future, so I hope you stick around even when your current project is complete If you’d like / be happy for your filter to be made generally available, you may be interested in having a look at this topic about suggested guidelines for Audacity Nyquist plug-ins: https://forum.audacityteam.org/t/conventions-for-nyquist-plug-ins/16324/1
If you’d prefer to just post the ‘raw’ filter code that’s fine (though it would be helpful if you could explicitly state if you intend it to be GPL v2 compatible)

What an excellent resource - thanks.

Happy to GPL v2 the code.

I’ll check out the guide and try to comply. As you might have noticed my lisp is letting me down! Just learned it over the last couple of weeks so might needs some help “optimising” the code.

I’m only self taught with Lisp (with some help from Edgar Franke and a lot of help from the documentation) but I’m happy to help where I can.
In case you’ve not got these: