Nyquist frequency divider feedback not working

Hi all,

I’m new to Nyquist, though I have done some DSP programming before. I have the following code working very well in Faust DSP, but am having a lot of trouble getting it to work in Nyquist. It’s basically a type of mixer/oscillator that precisely divides an audio frequency by two.

The code divides a frequency by two by mixing it with the downmixed output. In this case a 16kHz sinewave input is mixed with the 8Khz downmixed output by feeding the filtered 8kHz output signal back and multiplying it with the 16kHz input to produce the 8kHz difference. This requires the previous output sample of the routine to be sent back to the first stage for multiplication. So far I haven’t found a way to do this in Nyquist that actually works.

The dither ensures the process starts reliably. The working routine needs a limiter in the feedback path to prevent overload. This has been omitted from the example below for simplicity.

It appears as if there is no feedback happening. If the main expression is changed to:
set d = bandpass2(a * (0 + (noise() * dither-level)),8000,100) * 10
… it produces the same result.

I suspect this is due to Nyquist’s “Lazy evaluation”… but how can it be made to work?

The code is shown below.

set a = hzosc(16000,*table*,0)		;generate a 16kHz sinewave input signal for testing

set ditherlevel = db-to-linear(-50)	;set dither level

if ! fboundp(quote(d)) then set d = 0	;set d to zero only if d is unbound

set d = bandpass2(a * (d + (noise() * ditherlevel)),8000,100) * 10	;d = BPF(a * (d + dither)) * 10

return d

Try multiplying by 1000 instead of by 10. in the 4th line (or use a much higher level of noise).
Is that what you want it to do?

No, not at all.

Boosting the gain by increasing the multiplier just gives you higher level filtered noise (dither).

The noise is only there to get things started (like in an analog oscillator). If there was no dither used, there would initially be no output from multiplying a and d. In the Faust DSP language on a 64 bit PC the dither level can be set to -300dB and the code still produces a full scale sinewave output, requiring AGC/Limiting to avoid distortion.

If the routine is working correctly, the Q multiplier effect of the feedback produces a very pure 8kHz sinewave that rapidly grows in amplitude until filter overload occurs. This is not happening.

Swapping d for 0 in the main expression makes no difference. There appears to be no feedback happening. Am I using the variable “d” correctly?

One small correction… in the text above the code “dither-level” should be “ditherlevel”.

In the BANDPASS2 implementation, the biquad parameters are defined to keep the filter stable even at very high Q values.
If you want the filter to become unstable and go into self oscillation, then I think you will need to calculate the biquad parameters yourself and use BIQUAD or BIQUAD-M functions rather than the BANDPASS2. See: http://www.cs.cmu.edu/~rbd/doc/nyquist/part8.html#index460

On reflection, that’s not going to work. The biquad function is protected against becoming unstable and will throw an error.

The error checking for the biquad filter is in the function NYQ:BIQUAD.
You can redefine NYQ:BIQUAD in your own code without the error checking, and then you will able to use BIQUAD or BIQUAD-M with unstable parameters.
This is the definition of NYQ:BIQUAD that you need (it’s in LISP format as I’m not very familiar with SAL)

(defun nyq:biquad (x b0 b1 b2 a0 a1 a2)
  (if (<= a0 0.0)
      (error (format nil "a0 < 0 (unstable parameter a0 = ~A) in biquad~%" a0)))
  (let ((a0r (/ 1.0 a0)))
    (setf a1 (* a0r a1)
          a2 (* a0r a2))
    (snd-biquad x (* a0r b0) (* a0r b1) (* a0r b2)
                  a1 a2 0 0)))

Or you could use RESON (http://www.cs.cmu.edu/~rbd/doc/nyquist/part8.html#index447)

Filter stability refers to the internal filter operation (whether there is output with no input). The filter used here is required to be stable, or else the output frequency would be uncontrolled, and also there would be no way to apply AGC to prevent distortion.

So back to the original question… why isn’t the variable “d” providing feedback? Would “d” in the expression always be equal to the previous output “d” sample value? (it needs to be for the routine to work) Is Nyquist always evaluating the “d” variable in the expression, or is it evaluating it once only?

Is there a way to make sure Nyquist is re-evaluating “d” in the expression for each sample value passing through the routine?

This code:

set d = bandpass2(a * (d + (noise() * ditherlevel)),8000,100) * 10   ;d = BPF(a * (d + dither)) * 10

could be written as:

set input = a * (d + (noise() * ditherlevel))
set d = bandpass2(input, 8000, 100) * 10

In fact, your code is identical to this:

set input = hzosc(16000) * (0 + (noise() * db-to-linear(-50)))
return bandpass2(input, 8000, 100) * 10

or as one line:

return bandpass2( hzosc(1600) * noise() * db-to-linear(-50), 8000, 100) * 10

The Nyquist filter:

bandpass2(<sound>, 8000, 100)

is equivalent to the Matlab biquad filter with parameters:

b0: 0.00454278
b1: 0
b2: -0.00454278
a0: 1.00454
a1: -0.83554
a2: 0.995457

so your code could be written as (in SAL):

set b0 = 0.00454278
set b1 = 0
set b2 = -0.00454278
set a0 = 1.00454
set a1 = -0.83554
set a2 = 0.995457
set input = hzosc(16000) * (noise() * db-to-linear(-50))
return biquad-m(input, b0, b1, b2, a0, a1, a2) * 10

Those later examples have no feedback path, so they can’t perform the required function.

To generate output at 8kHz from a 16kHz input, the 16kHz signal has to be multiplied with an 8kHz signal to perform frequency mixing.

To ensure the output is exactly half the input frequency, the 8kHz mixing signal needs to be from the 8kHz output signal. This requires a path from the signal output back to the signal input.

In Faust DSP, you can’t do this by just using the variable “d”, as it flags an infinite loop error. Instead, the recursion operator “~” and the wire connection operators “_” must be used. How would this be done in Nyquist?

Here is a Faust DSP routine that actually works. (It has some aliasing distortion due to the simplified limiter code.) It can be run using the Faust online editor. Be warned, it produces a loud 8kHz audio tone when run.

import("stdfaust.lib");
process = d;
			
dither = no.noise * ba.db2linear(-300); //dither at -300dB is adequate for 64 bit systems
a = os.oscp(16000,0); //16kHz sinewave input test signal
d = (a * (_ + dither):fi.bandpass(2,7950,8050):limit) ~_; //frequency divider by mixing
limit(x) = x*(0.2/max(abs(x),0.09)); //simple limiter to avoid massive filter overload

OK, I see what you are wanting to do, and it’s quite tricky in Nyquist. See this section of the Nyquist manual: http://www.cs.cmu.edu/~rbd/doc/nyquist/part6.html#49
Of course, if all you wanted was the 8kHz sine wave, then in practical code you would just write hzosc(8000) [ or in LISP syntax, (hzosc 8000) ]

Here’s a block diagram of the Faust code in my previous post:

In the application this code is used in, a Hilbert clipper is normally used to hard limit the signal without causing harmonic distortion or aliasing. That wasn’t done here because it would complicate the code.

In this basic example, taking the output from the filter directly works better, as it reduces aliasing on the output caused by the harmonic distortion of the simplified limiter (a peak clipper) that has been used.

import("stdfaust.lib");
process = d;
         
dither = no.noise * ba.db2linear(-300); //dither at -300dB is adequate for 64 bit systems
a = os.oscp(16000,0); //16kHz sinewave input test signal
d = (a * (_ + dither):fi.bandpass(2,7950,8050)) ~ limit:_; //frequency divider by mixing
limit(x) = x*(0.5/max(abs(x),0.1)); //simple limiter to avoid massive filter overload

Here’s the block diagram of this code:

Can Nyquist do this? What function would allow this type of recursion?

Thanks for the link Steve. I have read those pages before, but I will study them more closely.

I am amazed how difficult it is in Nyquist to get one previous output sample to feed back into a function input! Maybe this needs to be addressed in a future release?

This code is to be part of a high performance AM synchronous demodulator. It uses a squaring loop to allow carrier recovery from the modulation sidebands even during carrier fades. But a squaring loop doubles the carrier frequency, so a divide by two is needed. This routine seems to be the simplest and cleanest way to precisely divide a frequency by two in DSP.

Faust code for all this is working quite well, but it would be useful to be able to process recorded IQ files from SDR receivers faster than real-time. Faust code can be compiled into VST plugins, but there seems to be a problem getting them to run in Audacity under Windows. Hence my interest in Nyquist…

The issue is that you are needing to create a new DSP block that incorporates a mixer and limiter into the inner loop of a a bandpass filter. Because Nyquist in an interpreted language, looping through samples is expensive, which is why Nyquist’s DSP blocks are written in C.

In Faust, DSP blocks are created by defining the algorithm in an algorithmic language, and then compiling that algorithm into low level code (such as C or C++). This is the main purpose of Faust.


In Nyquist you can do the same. The algorithm can be defined as an ALG definition, then compiled into C code Appendix 1: Extending Nyquist
This is not the main purpose of Nyquist, and is much less convenient to do than in Faust, but it can be done.

It is perfectly possible to create low distortion hard limiting, using Nyquist’s built-in DSP (no need to create new DSP blocks). https://github.com/audacity/audacity/blob/master/plug-ins/limiter.ny
However, you don’t want to do that. You specifically want to create a DSP in the way that you do in Faust, which while possible, is not easy in Nyquist.
Why do you want to use Nyquist for this, when it is so much easier in Faust?


In Nyquist, the easiest way to create new DSP blocks is as Lisp objects. See DSP in Lisp. This is considered to be an “advanced topic” for Nyquist programming.

My previous post answers that question:

It seems the way Nyquist processes audio isn’t compatible with writing routines that require a feedback path, unless you can compile your own routines as extensions.

If only Audacity offered Faust support! That would enable very powerful processing routines to be easily written.

Perhaps I should abandon Nyquist and focus on getting Faust VSTs to work in Audacity…

Faust is very good, but it’s execution options are a bit limited. In Windows the Faust VSTs compiled using the online tools won’t load in Audacity, so the Faust code can only be executed in real time - using the online Faust tools, in FaustLive, or compiled into Webassembly and running in a browser. These methods can’t operate directly on files, which means additional software might be needed for routing audio to and from the applications, and there can be audio dropouts due to other resource demands on the PC. It also means the code always runs in real-time, so a 1 hour file takes 1 hour to process.

or as a DSP “class” (as described in the link that I gave previously)

I rarely use Windows, but on Linux, LV2 plug-ins generated by the Faust on-line compiler work in Audacity.
It may be worth contacting the Faust people to see if they know why their VST effects don’t work in Audacity, or if they know a way to compile them in a way that is supported by Audacity.