Goertzel DSP implementation #1

theloni-monk · 2022-07-12T22:39:02Z

I think this mostly does what I want it to. Could probably use more testing.

WeirdConstructor

I've made some comments in the code for whoever picks it up. I just want to make sure the overall code quality of the DSP code in HexoDSP stays consistent.

WeirdConstructor · 2022-07-13T06:00:59Z

src/dsp/goertzel.rs

+
+    pub fn reset(&mut self){
+        self.buff.clear();
+        self.ideal_buffsize = (2.0 * (1.0/self.srate) / self.reference_tuning).floor() as usize;


Store sample rate not as 1.0/sample_rate, to save a division and be consistent with the other DSP implementations.

WeirdConstructor · 2022-07-13T06:04:29Z

src/dsp/goertzel.rs

+        let mut Q1 = 0.0;
+        let mut Q2 = 0.0;
+
+        for i in 0..self.ideal_buffsize{


tick() currently iterates all samples for each sample in the current block. That is quite expensive. It should be explored, if there is a significant downside to executing the Goertzel algorithm only each M samples.
Maybe with a user configurable M as "latency" parameter?

See also my comment below about the Gz3Filt.

WeirdConstructor · 2022-07-13T06:05:11Z

src/dsp/goertzel.rs

+            Q2 = Q1;
+            Q1 = Q0;
+        }
+        let mag_squared = (Q1.powf(2.0) + (Q2.powf(2.0)) - (Q1*Q2*coeff)) as f32;


WeirdConstructor · 2022-07-13T06:07:04Z

src/dsp/goertzel.rs

+        let mut Q1 = 0.0;
+        let mut Q2 = 0.0;
+
+        for i in 0..self.ideal_buffsize{


Change this loop to iterate over a slice of self.buff, that saves the overhead of checking unwrap() and reflects the intent better:

for s in &self.buff[0..self.ideal_buffsize] { Q0 = coeff * Q1 - Q2 + *s }

WeirdConstructor · 2022-07-13T06:10:26Z

src/dsp/goertzel.rs

+    pub srate: f32,
+    reference_tuning: f32, // assumed that target freqs will be int multiples of ref tuning
+    ideal_buffsize: usize, // calculated with respect to target freq and ref tuning
+    buff: VecDeque<f32>


VecDeque is unlikely in the current implementation to reallocate, which is fine.
But maybe should be replaced by a simpler implementation that writes into a fixed array similar to helpers::DelayBuffer - maybe even use the DelayBuffer in the first place. That makes accidental reallocations more unlikely.

WeirdConstructor · 2022-07-13T06:16:48Z

src/dsp/node_goertzel.rs

+           || (cgain - self.ogain).abs() > 0.0001
+        {
+            // recalculate coeffs of all in the cascade
+            self.computer.target_freq = cfreq;


Set frequency with a set method on the Goertzel structure.

WeirdConstructor · 2022-07-13T06:18:20Z

tests/node_goertzel.rs

+    let (out_l, _) = run_for_ms(&mut node_exec, 25.0);
+    let rms_minmax = calc_rms_mimax_each_ms(&out_l[..], 10.0);
+    eprintln!("RMS: {:?}", rms_minmax);
+    assert!(rms_minmax[1].2 - rms_minmax[1].1 < 0.01); // the output should be const for const freq input


Improve this test case and check the expected output with an absolute value additionally to this min/max difference check.

WeirdConstructor · 2022-07-13T06:20:21Z

tests/node_goertzel.rs

+
+    (matrix, node_exec)
+}
+


Add test cases:

that check the Goertzel output against noise as input.

that change the Goertzel frequency from it's default and checks if the output signal of the Goertzel algorithm changes.

WeirdConstructor · 2022-07-13T06:22:23Z

src/dsp/node_goertzel.rs

+use crate::dsp::goertzel::*;
+
+#[macro_export]
+macro_rules! fa_goertzel_type { ($formatter: expr, $v: expr, $denorm_v: expr) => { {


remove this, fa_goertzel_type is not used anywhere.

WeirdConstructor · 2022-07-13T06:23:37Z

tests/node_goertzel.rs

@@ -0,0 +1,51 @@
+// Copyright (c) 2021 Weird Constructor <[email protected]>


fix Copyright to that of goertzel.rs

WeirdConstructor · 2022-07-13T06:32:48Z

I have some concerns regarding the performance of this filter. Currently it recalculates the output over the complete 100 or 200 samples in it's buffer for every input sample. This should be checked in HexoSynth and the DSP CPU should be observed for this filter (eg. by testing as plugin in a DAW or via qjackctl). And maybe a note about the performance should be made in the help text for this node.

This is part of my DSP node checklist for adding this to HexoDSP:

Implement boilerplate (done)
Document boilerplate (done)
DSP implementation (done, except necessary improvements)
Parameter fine tuning (test in HexoSynth, if the parameters feel right)
DSP tests for all params
Ensure Documentation is properly formatted for the GUI
Add CHANGELOG.md entry to HexoSynth
Add table entry in README.md in HexoSynth

WeirdConstructor · 2022-07-13T07:01:33Z

Regarding the performance, I could see a different kind of node implementation as valuable for using it in a sound synthesis patch:

Have a Gz3Filt, which provides 3 outputs for 3 diffrerent frequencies
With 3 input parameters to configure the frequencies that should be detected.
Have an input setting parameter that defines the static resolution (the implied output delay) of the filter. Have the user decide if they want the result after X=1, 2, 4, 10, ... milliseconds. And mention in the documentation the accuracy implications this has (the longer they can wait, the more accurate the detection).
The outputs should output the most recently detected value and change each X milliseconds
Apply a user configurable slew limiter to each output as convenience

This exploits the efficiency of the Goertzel algorithm better in context of a modular synthesizer in my opinion. Because you get a precise (depending on the output delay settings) very fine multi band bandpass filter.

theloni-monk · 2022-07-13T11:40:23Z

Have an input setting parameter that defines the static resolution (the implied output delay) of the filter. Have the user decide if they want the result after X=1, 2, 4, 10, ... milliseconds. And mention in the documentation the accuracy implications this has (the longer they can wait, the more accurate the detection).
The outputs should output the most recently detected value and change each X milliseconds

I'm interested in this implementation, so I have some questions. Feel free to answer as many as would not be annoying. For this would it be too expensive to have that latency but window it such that the value still changes every frame, i.e. take overlapping windows and compute. Also for the input frequencies, would these be parameters or inputs? could they be both? I agree the GzFilt3 would be more interesting. Are CV signals restricted to (-1,1) and we would have to map that to a frequency range? or can they be whatever?

WeirdConstructor · 2022-07-13T12:07:45Z

I'm interested in this implementation, so I have some questions. Feel free to answer as many as would not be annoying. For this would it be too expensive to have that latency but window it such that the value still changes every frame, i.e. take overlapping windows and compute.

How would the window work? The performance hit comes from calculating Q0, Q1 and Q2 for the whole buffer again and again per input sample. You can of course amortize this by recalculating only each 5, 10 or 50 samples. But that is not as cheap as just accepting the longer delays and run goertzel Q0/Q1/Q2 calculations only once per new input sample.

I believe, correct me if I am wrong, that is also the whole point for Goertzel vs. FFT, as detecting only a few known frequencies does not require a full FFT anymore, but only running Q0/Q1/Q2 calculation once per input sample.

Also for the input frequencies, would these be parameters or inputs?

In HexoDSP there are two kinds of inputs for a node: Input parameters inp and atoms at. Input parameters get the iconic octagon knobs and accept sample accurate inputs from other nodes. Atoms are settings that are not automateable, they are settings that are manually changed by a user typically.

could they be both? I agree the GzFilt3 would be more interesting. Are CV signals restricted to (-1,1) and we would have to map that to a frequency range? or can they be whatever?

All signals in the DSP graph are (kind of) restricted to [-1,1]. "Kind of" because that is more a convention than a technical restriction. Even though most nodes clamp their inputs to that range.
The mapping between the DSP graph signal and more meaningful values that are used by the actual mathematic calculations is done by "denormalizing" the values. This is specified in the big parameter/input matrix in mod.rs and hidden in statements like these:

   let freq = inp::BOsc::freq(inputs);

theloni-monk · 2022-07-13T12:12:59Z

But that is not as cheap as just accepting the longer delays and run goertzel Q0/Q1/Q2 calculations only once per new input sample. I believe, correct me if I am wrong, that is also the whole point for Goertzel vs. FFT, as detecting only a few known frequencies does not require a full FFT anymore, but only running Q0/Q1/Q2 calculation once per input sample.

Yes and no... so the more samples you take the thinner the bins of the equivalent fft become, but the more you are averaging over the sample. So computing a new Q0/Q1/Q2 for each new sample instead of each new buffer is eqivalent to finding the amplitude that a given frequency has had averaged over the entire life of the node. This is why you want to compute it for a window over the input stream, instead of the whole input stream. Currently, I am making a new window every sample, this is needlessly expensive however. So as you say, it would make sense to have a param for how often to calculate with respect to a new window.

WeirdConstructor · 2022-07-13T12:44:35Z

Yes and no... so the more samples you take the thinner the bins of the equivalent fft become, but the more you are averaging over the sample. So computing a new Q0/Q1/Q2 for each new sample instead of each new buffer is eqivalent to finding the amplitude that a given frequency has had averaged over the entire life of the node. This is why you want to compute it for a window over the input stream, instead of the whole input stream. Currently, I am making a new window every sample, this is needlessly expensive however. So as you say, it would make sense to have a param for how often to calculate with respect to a new window.

Yes, I would love to give the user that control of how big the window is, which means how accurate the frequency they pass in via parameter is detected. And I would see Gz3Filt not using overlapping windows (which requires saving input samples in a buffer, which is often used for visualizing a spectrum of frequencies), but by resetting Q0/Q1/Q2 to 0 (like described in the website you linked) and restarting for the next window. That would also get rid of any internal buffering.

The output would then be the most recently computed amplitude of the frequency.

theloni-monk · 2022-07-13T12:58:12Z

Sounds good, I'll work on these over the next week or so when I have free time

theloni-monk and others added 7 commits July 9, 2022 15:30

started

856b35a

Merge remote-tracking branch 'upstream/master'

dbbfecf

wrote basic node... untested

986d0d4

it builds now

5144ed4

started test and added with_capacity by default

50aaf77

wrote and passed tests

6e669a8

added dynamic sample rate

d6b9ca2

WeirdConstructor reviewed Jul 13, 2022

View reviewed changes

theloni-monk added 9 commits July 13, 2022 21:35

eliminated private buffer

5f82c8b

implemented 3 gfilters

501662f

added latency param

9a1c934

started debugging tests

7023313

cache calculation on previous window

1735adf

merge formatting

c5abf7f

fix merge error

3c20d17

multi-output -bugged

178f20d

Merge remote-tracking branch 'upstream/master'

a280c57

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Goertzel DSP implementation #1

Goertzel DSP implementation #1

theloni-monk commented Jul 12, 2022

WeirdConstructor left a comment

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor Jul 13, 2022

WeirdConstructor commented Jul 13, 2022 •

edited

Loading

WeirdConstructor commented Jul 13, 2022 •

edited

Loading

theloni-monk commented Jul 13, 2022

WeirdConstructor commented Jul 13, 2022

theloni-monk commented Jul 13, 2022 •

edited

Loading

WeirdConstructor commented Jul 13, 2022

theloni-monk commented Jul 13, 2022

		@@ -0,0 +1,51 @@
		// Copyright (c) 2021 Weird Constructor <[email protected]>

Goertzel DSP implementation #1

Are you sure you want to change the base?

Goertzel DSP implementation #1

Conversation

theloni-monk commented Jul 12, 2022

WeirdConstructor left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

WeirdConstructor commented Jul 13, 2022 • edited Loading

WeirdConstructor commented Jul 13, 2022 • edited Loading

theloni-monk commented Jul 13, 2022

WeirdConstructor commented Jul 13, 2022

theloni-monk commented Jul 13, 2022 • edited Loading

WeirdConstructor commented Jul 13, 2022

theloni-monk commented Jul 13, 2022

WeirdConstructor commented Jul 13, 2022 •

edited

Loading

WeirdConstructor commented Jul 13, 2022 •

edited

Loading

theloni-monk commented Jul 13, 2022 •

edited

Loading