It is sometimes necessary to use very steep low pass filters, for example for preventing aliasing when resampling or pitch shifting.
Nyquist provides a fairly steep 8th order low pass filter (48 dB/Oct) with:
(lowpass8 sound frequency)
If this is not steep enough for a particular application, the filter may be “nested” so that it is applied multiple times, for example:
(lowpass8
(lowpass8 sound frequency)
frequency)
Nesting up to about 10 lowpass8 filters can be used to good effect by looping through the filter:
;;; steep low pass filter
(defun lp-steep (sound hz)
(setq hz (* 0.94 hz))
(dotimes (count 10)
(setf sound (lowpass8 sound hz)))
sound)
; apply a steep 1000 Hz filter to the selection
(lp-steep s 1000)
This provides a pretty steep cut-off, though there is a fair amount of attenuation before the set filter frequency.
The lp-steep function lowers the cut-off frequency to 94% of its input value so as to provide a high degree of attenuation at the set frequency.
For an even sharper cut off, here is a 100th order Butterworth filter.
The characteristics are that it has an extremely steep cut-off with very high attenuation at the filter frequency (around -90 dB at 1 kHz).
Many thanks to kai.fisher for doing most of the leg work on this https://forum.audacityteam.org/t/chebyshev-type-i-and-butterworth-hp-lp-filters/24241/1 and to S. Smith, author of “The Scientist and Engineer’s Guide to Digital Signal Processing” for publishing the biquad equations on which it is based.
(defun lp100 (sig frequency)
(defun butterworth (s1 fc n)
(defun Coeffs (ffc p n)
(let* (;Calculate pole location on unit circle in z-plane
(ReP (- (cos (+ (/ pi (* n 2.0)) (* (- p 1.0) (/ pi n))))))
(ImP (sin (+ (/ pi (* n 2.0)) (* (- p 1.0) (/ pi n)))))
;s-domain to z-domain conversion
(t (* 2.0 (tan 0.5)))
(t2 (* t t))
(wfc (* 2.0 pi ffc))
(m2 (+ (* ReP ReP) (* ImP ImP)))
(d (+ 4.0 (- (* 4.0 ReP t)) (* m2 t2)))
;Lowpass prototype coefficients
(x0 (/ t2 d))
(x1 (/ (* 2.0 t2) d))
(x2 (/ t2 d))
(y1 (/ (- 8.0 (* 2.0 m2 t2)) d))
(y2 (/ (- -4.0 (* 4.0 ReP t) (* m2 t2)) d))
;Transform Lp prototype filter coefficients to actual LP or to HP
(k (/ (sin (- 0.5 (/ wfc 2.0))) (sin (+ (/ wfc 2.0) 0.5))))
(k2 (* k k))
(d (+ 1.0 (* y1 k) (- (* y2 k2))))
;derive coefficients
(a0 1.0)
(a1 (/ (+ (* 2.0 k) y1 (* y1 k2) (- (* 2.0 y2 k))) d))
(a2 (/ (+ (- k2) (- (* y1 k)) y2) d))
(b0 (/ (+ x0 (- (* x1 k)) (* x2 k2)) d))
(b1 (/ (+ (- (* 2.0 x0 k)) x1 (* x1 k2) (- (* 2.0 x2 k))) d))
(b2 (/ (+ (* x0 k2) (- (* x1 k)) x2) d))
;Determine gain of filter
(ga (+ a0 (- a1) (- a2)))
(gb (+ b0 b1 b2 ))
(g (/ ga gb))
;Apply gain so that passband is unity.
(b0 (* b0 g))
(b1 (* b1 g))
(b2 (* b2 g)))
;Output coefficients
(vector a0 a1 a2 b0 b1 b2)))
;convert cut-off frequency to a fraction of sampling rate
(let ((ffc (/ (* 0.94 fc) *sound-srate*)))
;loop once for each biquad
(dotimes (p (/ n 2))
(setf coef (Coeffs ffc (1+ p) n))
(setf a0 (aref coef 0))
(setf a1 (aref coef 1))
(setf a2 (aref coef 2))
(setf b0 (aref coef 3))
(setf b1 (aref coef 4))
(setf b2 (aref coef 5))
(setf s1 (biquad s1 b0 b1 b2 a0 a1 a2)))
s1))
; 100th order butterworth
(butterworth sig frequency 100))
; apply a 'brick wall' 1000 Hz filter to the selection
(lp100 s 1000)