No wonder that it looks alien 
From Audacity bug 152:
I can confirm that with unprotected biquad and lowpass2 functions, garbage is returned with lowpass2 frequencies > samplerate/2 but I do get a random garbage noise, not -1 samples. I hope that the -1 samples problem was caused by any Audacity weirdness that was fixed as Martyn wrote.
With the new “dspprims.lsp” from post 180014 on page 2 of this thread the highpass2, lowpass2 as well as eq-lowshelf functions work as expected, with highpass2 and lowpass2 frequencies > samplerate/2 I get an “frequency out of range” error, what is absolutely OK.
In case that -1 samples appear again, I wrote a Lisp replacement for the C code in “biquadfilt.c”, where the biquad filter coefficients and the sample values of each filter iteration can be printed in the Debug window:
(defun nyq:biquad (snd b0 b1 b2 a0 a1 a2)
;; this is what the Lisp code in nyq:biquad does
(let ((a0r (/ 1.0 a0)))
(setf b0 (* a0r b0) b1 (* a0r b1) b2 (* a0r b2)
a1 (* a0r a1) a2 (* a0r a2)))
;; print the results of nyq:biquad
(format t "b0 ~a~%b1 ~a~%b2 ~a~%" b0 b1 b2)
(format t "a0 ~a~%a1 ~a~%a2 ~a~%" a0 a1 a2)
;; copy all samples from the sound into an array
(let* ((size (round len))
(snd-array (snd-fetch-array snd size size)))
;; this is what the C code in biquadfilt.c does
(let ((z0 0.0) ; current input sample
(z1 0.0) ; previous input sample
(z2 0.0) ; second-previous input sample
(out 0.0)) ; current output sample
(dotimes (index size)
;; fetch a sample from the array
(setf z0 (aref snd-array index))
;; compute the a coefficients
(setf z0 (+ z0 (* a1 z1) (* a2 z2)))
;; compute the b coefficients
(setf out (+ (* z0 b0) (* z1 b1) (* z2 b2)))
;; print the values of the current iteration
(format t "index ~a~%" index)
(format t " z2 = ~a~%" z2)
(format t " z1 = ~a~%" z1)
(format t " z0 = ~a~%" z0)
(format t " out = ~a~%" out)
;; write the output sample back to the array
(setf (aref snd-array index) out)
;; shift the values of the sample variables
(setf z2 z1)
(setf z1 z0))
;; return the sound from the array
(snd-from-array 0.0 *sound-srate* snd-array))))
Just copy the function into the Audacity Nyquist prompt and the normal Nyquist filter function call below. The function will overwrite the built-in nyq:biquad function and use Lisp code instead of C code to compute the DSP signal. But be aware that this is terribly slow, so it’s better not to use this function with signals of more than approx. 1000 (one thousand) samples length!
Yes, but the problem is that the trap would need to be executed with every sample anew, what would make the filter much slower. In case the -1 samples appear again it would be better to find a way how to test the a and b coefficients before they are passed to the C function.
The problem with C code is that it always looks alien, even if written by humans.
- edgar
(highpass2 s 24000) applied to a track with a 44.1 kHz sample rate.
The sample values are rapidly escalating.
By the 500th sample we have
index 500
z2 = 1.0835e+46
z1 = -1.08191e+46
z0 = 9.74385e+45
out = 1.02527e+45
by the 1000th sample:
index 1000
z2 = 5.2483e+92
z1 = -6.30421e+92
z0 = 7.31081e+92
out = 6.11209e+91
and by 1999 (as far a I’ve taken it)
index 1999
z2 = -2.06755e+184
z1 = 8.61847e+184
z0 = -1.78256e+185
out = -9.01728e+183
When we get up to 3340 samples we pass infinity.
index 3338
z2 = -5.64808e+307
z1 = 1.31419e+308
z0 = -inf
out = -inf
index 3339
z2 = 1.31419e+308
z1 = -inf
z0 = -nan
out = -nan
index 3340
z2 = -inf
z1 = -nan
z0 = -nan
out = -nan
index 3341
z2 = -nan
z1 = -nan
z0 = -nan
out = -nan
index 3342
z2 = -nan
z1 = -nan
z0 = -nan
out = -nan
by the way, it runs a lot faster if it outputs the text to a file rather than the debug window.
(defun nyq:biquad (snd b0 b1 b2 a0 a1 a2)
(setq outfile "/home/<user>/Desktop/nyqbiquad-data.txt")
(setq f (open outfile :direction :output))
;; this is what the Lisp code in nyq:biquad does
(let ((a0r (/ 1.0 a0)))
(setf b0 (* a0r b0) b1 (* a0r b1) b2 (* a0r b2)
a1 (* a0r a1) a2 (* a0r a2)))
;; print the results of nyq:biquad
(format f "b0 ~a~%b1 ~a~%b2 ~a~%" b0 b1 b2)
(format f "a0 ~a~%a1 ~a~%a2 ~a~%" a0 a1 a2)
;; copy all samples from the sound into an array
(let* ((size (round len))
(snd-array (snd-fetch-array snd size size)))
;; this is what the C code in biquadfilt.c does
(let ((z0 0.0) ; current input sample
(z1 0.0) ; previous input sample
(z2 0.0) ; second-previous input sample
(out 0.0)) ; current output sample
(dotimes (index size)
;; fetch a sample from the array
(setf z0 (aref snd-array index))
;; compute the a coefficients
(setf z0 (+ z0 (* a1 z1) (* a2 z2)))
;; compute the b coefficients
(setf out (+ (* z0 b0) (* z1 b1) (* z2 b2)))
;; print the values of the current iteration
(format f "index ~a~%" index)
(format f " z2 = ~a~%" z2)
(format f " z1 = ~a~%" z1)
(format f " z0 = ~a~%" z0)
(format f " out = ~a~%" out)
;; write the output sample back to the array
(setf (aref snd-array index) out)
;; shift the values of the sample variables
(setf z2 z1)
(setf z1 z0))
(close f)
;; return the sound from the array
(snd-from-array 0.0 *sound-srate* snd-array))))
(highpass2 s 24000)
; test sample values
; (setq outfile "/home/<user>/Desktop/sample-data.txt")
; (setq f (open outfile :direction :output))
; (dotimes (i (snd-length s ny:all))
; (format f "~a~%" (snd-fetch s)))
; (close f)
; (format nil "Done")
Now be prepared for one of Edgar’s famous box-and-arrow pictures. To me the biquad DSP code from above looks like this:
+------- * a2 --------------+
| +----- * a1 -----+ |
v v | |
+---+ +----+ +----+ +----+
in -->| + |-->| z0 |-->| z1 |-->| z2 |
+---+ +----+ +----+ +----+
| | |
* b0 * b1 * b2
v v v
+---------------------+
| + |--> out
+---------------------+
The current z0 sample is computed from the sum of the input sample together with two feedback signals from the previous z1 and the second-previous z2 samples, multiplied by the a1 and a2 coefficients. This means that only a1 and a2 can produce feedback desaster, but not b0, b1, or b2, so IMO feedback overflow can only happen if a1 or a2 are greater than 1, where I’m afraid that this test is not enough, because the overflow probably depends on a combination of both values.
How does the picture compare e.g. to the Wikipedia biquad math (translated to Nyquist variable names):
b0 + b1*z1 + b2*z2
out = --------------------
1 + a1*z1 + a2*z2
And what exactly does the a0r “normalisation” in the original nyq:biquad function do? Is this “normalisation” the key to transform the Wikipedia math to the code in the picture above? I can see similarities but I still do not understand the details.
Don’t ask me why all this produces an audio filter, I first need to find out more about the practical implementation of biquad filters.
- edgar
Intermediate result - nyq:biquad produces infinite feedback if after normalising a1 and a2, (abs (* a1 a2)) is greater than 2, or in code:
(defun nyq:biquad (x b0 b1 b2 a0 a1 a2)
(let ((a0r (/ 1.0 a0)))
(setf a1 (* a0r a1) a2 (* a0r a2))
(if (< 2 (abs (* a1 a2)))
(error "NYQ:BIQUAD - parameters out of range")
(snd-biquad x (* a0r b0) (* a0r b1) (* a0r b2) a1 a2 0 0))))
This seems to work but I do not recommend to implement it before I have undestood why.
- edgar
Bad news: The code from above, as somehow already expected, does only work if the coefficients are preprocessed e.g. by trigonometric functions in the highpass2 and lowpass2 functions, so they already are in a semi-meaninful range before they enter the nyq:biquad function. If random numbers are used as a and b coefficients then it is sufficient that only one of the normalized a1 or a2 coefficients is greater than 1.0 to make snd-biquad run towards positive or negative infinity and it only matters if the number of iterations is big enough to make a floating-point overflow happen.
I have talked to a physics engineer who said that because the biquad theorem is just a mathematical abstraction the problem is the same as with nearly all mathematical functions. The possible input parameter combinations that produce meaningful results are only a tiny fraction compared to the much higher number of random combinations that produce completely useless nonsense results.
And it’s even worse. Because IIR filters are based on recursion it is impossible to predict from the values of the input parameters if the result will be useful at all. The only reliable way to find this out is first to compute the result and then (afterwards) test if the result is useful or not. This would completely contradict the Nyquist “lazy evaluation” principle, where samples are only computed if really needed.
The only other way to fix snd-biquad would be to use floating-point traps in the inner DSP loops, what would make snd-biquad a lot slower, so this is also not a solution that makes the slightest sense.
Conclusion: What would be possible is to add samplerate/2 tests to all high-level Lisp filter functions that call the nyq:biquad or nyq:biquad-m functions. Frequency values above samplerate/2 would lead to aliased (=useless) signals anyway, and the samplerate/2 tests in the highpass2 and lowpass2 functions have already proved to be useful, so I suggest this as the way to go.
Nyquist investigations finished for today - I’m tired now.
- edgar
New version of shelf filter that should work with Audacity 2.0
shelf.ny (2.75 KB)
Could someone test the latest version before I upload it to the wiki. Thanks
Latest version available on the wiki: http://wiki.audacityteam.org/wiki/Nyquist_Effect_Plug-ins#Shelf_Filter
Nice work! appreciate it.
btw, would be good if you posted in the first post, that “You can download latest version from wiki[link]” so, instead of attached old file, link to wiki.
thnx
Looks to me like the first post (Shelf Filter) already does that. 