Making the BBD implementation correct

[jpcima]

@b3ll, I'll describe in some detail what I'm doing here, and I let you experiment what may work or not.
Considering the paper is not really clear at all times, I've taken a few guesses.

If the works reaches a point of apparent stability, I'll merge it.
For now, I've added 2 commits in the unit-delay branch.

• 971a7a8 for fixing the output filter

Filter couldn't operate according to the structure shown, because the unit delay variable was not kept, meaning either algorithm or structure has been incorrect.
Most probably algorithm, so I reformulated the filter with use of unit delays (as Xout_mem).

• e7f197e BBD's summed output, going to output filter

You have the processing made of 2 loops: (as denoted in the paper)

• the k loop is the outer loop, runs at the sample rate (being n in code, don't mistake it)
• the n loop the inner one, operating at BBD clock rate

n-loop implements as:

bbd-delay-experimental/bbd_line.cc

Line 72 in 543cd7e

for (unsigned tick = 0; tick < tick_count; ++tick) {

Further, this n-loop does two different jobs on even/odd ticks:

• even: pushes in a sample into 1 BBD stage
• odd: sums BBD-rate outputs for entry into output filter, which operates at audio-rate

You will observe odd-case will perform this add:

bbd-delay-experimental/bbd_line.cc

Line 89 in 543cd7e

Xout[m] += Gout[m] * delta;

I have found suspicious that BBD's output Xout can be allowed to grow unlimited, as this sum is never reset.

From the paper, I was able to read this information.

If I correctly interpret "over one interval of the k side", it suggest to reset the Xout sum at each frame. Hence I did.

bbd-delay-experimental/bbd_line.cc

Lines 70 to 71 in e7f197e

	for (unsigned m = 0; m < Mout; ++m)
	Xout[m] = 0;

Thinking of it, it means at some frames, accumulation may remain 0 if BBD doesn't tick (the case of a slow clock).
But also intuitively, it makes sense: in analog world, when BBD doesn't tick, the bucket will empty itself (capacitor discharge), and it's the similar idea when zero are fed into the lowpass following.

---

Hope this little comment makes sense.
I didn't try this action, you can tell me if this has worked or not.

[b3ll]

Ok, so I tried both of these, both together and individually and I didn't notice anything different, unfortunately. These things do seem like they're correct regardless (especially zeroing out the ever increasing output filter, and the emptying bucket when there's no clock ticks).

However, I did some more analysis and it looks like the issue is even worse when the delay time is being modulated.

Assuming that what I'm doing is correct (this should be a single delay line, with step size of 256, modulating a delay value from 0.00166 to 0.00535 over a period of 1hz with a simple sine wave:

  const float modFrequency = 1.0f;
  const float minDelay = 0.00166f;
  const float maxDelay = 0.00535f;

  float modDelayBase = (maxDelay + minDelay) / 2.0f;
  float modDelayDepth = (maxDelay - minDelay) / 2.0f;
  bucketSize = 256;
  auto *leftChannel = buffer.getWritePointer(0);
  float leftClock[numSamples];
  for (int i = 0; i < numSamples; i++) {
    float modDelayValue = sin(2 * M_PI * _leftTime);

    float delayTime = modDelayBase + (modDelayDepth * modDelayValue);

    leftClock[i] = ((bucketSize * 2.0) / delayTime) / sampleRate;

    _leftTime += (modFrequency / sampleRate);
    _leftTime -= (int)_leftTime;
  }
  _leftDelayLine->process(numSamples, leftChannel, leftClock);

It ends up still shelving a lot of artifacts in the top end of the spectrum (regardless of what note is being played). I attached a sound file with before / after, when you hear the delay being enabled you'll hear a buzzing noise (alongside the delayed signal). There's also a video showing the filter output after playing without the delay enabled, playing the same note again without the delay enabled, and then the delay being enabled. You'll see the top end is oscillating.

Part of me wants to think that somewhere along (perhaps the start of) the BBD implementation, some noise gets introduced into the delay and then that noise is persisted indefinitely. When I initially read this I thought that the Xout growing unlimited could've been the issue, but I don't think that's the cause. The buzzing is 100% related to the modulation since if I set a slower or faster modulation rate, the buzzing changes to match.

filter-output.zip

[jpcima]

For now I don't know, I'll just note a couple things.

• there exists a reference available for comparing. It has a description in the paper's results section, and sound samples on the webpage. https://www.hsu-hh.de/ant/en/team/martin-holters/dafx2018-bbd
  It will be useful to replicate the experimental protocol for verifying.

• the analog BBD will also have inherent aliasing, being also a sampled system.
  in digital, the model translates to a sampled system inside of another. (getting somewhat meta)

• oversampling by factor 2 is suggested to lower the aliasing. The website also has oversampling samples.
  

[b3ll]

Yeah, I was about to try oversampling to see what happens. I imagine changing the sample rate to 88200hz isn't oversampling right? To oversample, I should be taking the input, and interpolating the input to 88200hz (by doubling and averaging, right?)

i.e. [sample1, sample2, sample3, …] becomes [sample1, (sample1+sample2) / 2, sample2, (sample2+sample3) / 2, …]

Also, in the technical writeup they talk about having the ƒBBD at 50000hz (for their test sources), how would I go about doing that in your implementation? Would I be supplying the delay as 0.000020s without modulation?

[jpcima]

Yeah, I was about to try oversampling to see what happens. I imagine changing the sample rate to 88200hz isn't oversampling right? To oversample, I should be taking the input, and interpolating the input to 88200hz (by doubling and averaging, right?)

No, it's a little more elaborate. More like:

1. upsample input: fill gaps with zero samples, and smooth it with lowpass
2. setup the bbd element with N× times scaled frequencies (filters rate divided by N, BBD clocks divided by N also)
3. run then lowpass output
4. decimate by N (keep every Nth output)

Anyway, better to fix the ordinary case first, before doing the oversampling.

Also, in the technical writeup they talk about having the ƒBBD at 50000hz (for their test sources), how would I go about doing that in your implementation?

It's a simple matter of passing (50000/Fs) as clock value.
It handles the frequencies over Nyquist, it's not a problem to clock a line to be this high.

[b3ll]

Ok, now we're getting somewhere!

I setup the bbd as such:

  // setup
  _leftDelayLine = new BBD_Line(44100.0, 256, bbd_fin_j60, bbd_fout_j60);

  // inside process block
  auto *leftChannel = buffer.getWritePointer(0);
  float leftClock[numSamples];
  for (int i = 0; i < numSamples; i++) {
    leftClock[i] = 50000.0 / 44100.0;
  }
  _leftDelayLine->process(numSamples, leftChannel, leftClock);

  buffer.copyFrom(1, 0, buffer.getReadPointer(0), buffer.getNumSamples());

It turns out that when running this with the sample wav files supplied (they're attached), there's no sort of modulating / delaying effect heard at all. It just sounds super noisy and has lots of static added (which is probably what that static noise is when the BBD is modulated.

I looked at the sample files, and it looks like we're not making the same shape at all

Here's a photo of the original C chord input file:

Then this is what the research paper's output file looked like:

But this is what ours looks like (using what I wrote above):

I'm currently trying to understand the whole algorithm to see if there's something missing at the beginning or if there's an error in how the signal is being delayed.

Here's a zip of the input chord, their output model chord, and what we're generating (as a before / after).
chords.zip

[jpcima]

My next step in the plan is trying to implement "3.3 Real-valued systems".

The goal is to rewrite the complex-valued filter in a real-valued one.
It takes a modification of the filter structure.
(eg. the case M=5, it converts to 1 structure of first order and 2 structures of second order)

When it's real-valued, it's a return to a more conventional filter structure, and it makes the analysis easier.

Also now, the branch contains a test program which works on audio files.
The buzzing distorsion is reproducible, even on the simple case, a sine wave at 1kHz.

[b3ll]

Okay, that sounds great. I tried of altering a lot of things last night, but I couldn’t seem to get the buzzing to go away… oversampling did help but only with absurdly high sample rates (i.e. 192000)

I was trying to see how the filter was done, but I genuinely don’t understand how the complex values map to a filter cutoff point. Filters I’ve normally written are usually something like (anything greater than some frequency is dropped, and others are dropped by a certain percent, depending on the curve).

Is there anything I can read to understand how the complex values / filter works? I’d very much like to help make the real value one

[jpcima]

You haven't picked the simplest of exercises for a start. it's for sure :D

This filter discretization is very unusual, but I think it exploits the fact of BBD being a sampling system.
In this logic, it can treat signal input as a series of instantaneous pulses which occurs at the BBD ticks.

As I understand, it's because of this special case, that the filter can be simplified and expressed as such.

So yeah, considering it's so application-specific, you will have surely a hard time to get in the detail of this.

But if you search a generic explanation of filter theory, probably Miller Puckette's book has the best lesson I have found. (The Theory and Technique of Electronic Music)

[b3ll]

hah, yeah, it's fun to learn piece by piece though!

Upon checking out the book it has a section on delays which deals with complex numbers and it seems to share a lot of the same stuff as mentioned in the research paper, so I'll read that stuff too!

[jpcima]

Adam, I think the error is finally fixed at unit-delay, in the commit aadf660.
As it seems, the output filter was just generated with a wrong formula.
Please check with the new one.

[b3ll]

IT WORKS! 💯🔥🎹✨

And it sounds great!! Very nice catch! :D :D

[b3ll]

Ok I'm gonna get this working in the chorus effect plugin I'm working on and I'll post the results here :D

[jpcima]

Also, regarding the two fixes made regarding the algorithm.

(1) the reset of the sum must absolutely be made, or signal will go excentered around the DC, which is not the case in any of the website's output recordings

(2) about the filter restructuring and its unit delay node: it seems signal is outputting with more edge before the fix, and as more dull after applying. as it woulds seem, the website recordings appear to exhibit this tiny high frequency content.

So, who knows if it's algorithm which is right, or it's the structure pictured...
It's testable with or without patch, and see what will happen...

[b3ll]

Is it possible to mix the filter output? Perhaps its too harsh and clipping more than it should…

Here's ours (with the filter using 1024 as an interp_size)

vs. theirs (theirs has additional modulation that I'm not seeing in our delay though so that's kinda weird)

These look rather similar and I don't think it's clipping too much, it still sounds great and significantly better than the buzzing beforehand

[b3ll]

In addition, I ran the test program against the input_chord.wav and compared it to the model_chord.wav and they look nearly identical. The model_chord.wav does have some oscillating frequencies in the very top end (around 20k), but that's the only difference I'm seeing… the article only mentions a clock of 50kHz and nothing about modulation, so I'm not sure about that

input_chord.wav ➡️ their file
input_model.wav ➡️ input_chord.wav ran through bbd with 256 steps at 50kHz
model_chord.wav ➡️ their file

chords2.zip

[jpcima]

It's quite possible that "theirs" has just misimplemented the output filter, and it's explaining the high frequency content. It's as I said above.

Did you compare with 971a7a8 modifications reverted ?

[b3ll]

Yeah, reverting 971a7a8 introduces the buzzing issue and isn't even close

[jpcima]

Yep, It's the same observation here. By the algorithm's way of implementing, output filter is inoperational.

Is it possible to mix the filter output?

Do you mean, to have a wet/dry control? you can't make this on chorus side?

[b3ll]

Oh sorry, I meant asking that in relation to the "dullness" of the filter. I think it sounds good but I was proposing the idea of having the output filter having a wet/dry control to add some of that harshness / edge back if desired

[jpcima]

Well it was just the implementation error, which makes the filter disabled. 🤡

The thing about these filters, it's they're designed from a fixed specification, it's not directly controllable.
As you can see by BBD_Filter_Spec, and values are taken from article's "Table 1".

It's description of the "decomposed" filter form in M sections that the document's I/O filter sections are about. They must relate to filter's analog zeros and poles somewhat, but it's not told exactly how, or how are these values obtained from the Juno.

[b3ll]

LOL fair.

Yeah that part was definitely confusing to me, I had no idea how those values were computed (nor was it mentioned in the research paper).

Regardless, I'm happy with how this sounds! :D

[b3ll]

Here's some more good news, this is my plugin which aims to repro setting number II on the chorus unit vs. a Roland System 8 using its chorus II setting which uses hardware modeling to be nearly identical to the real unit:

Both are playing the note C3 endlessly, that's pretty close

[b3ll]

Also here's a C Major Chord

Mine:

Theirs:

[jpcima]

I've tried plotting these I/O filters in numpy according to their R and P values.
These describe analog transfer functions, R stands for residue, and P for pole.

But the bode plots came out as inverted, I don't know what's up with that.
Here it is: filter.py.gz

[b3ll]

uh wait what, that's really weird… if this is inverted, that would explain the tail of the frequency spectrum getting magnified a lot, no?

[jpcima]

No, it means I don't have right interpretation of these values, or I don't know to use scipy software correctly 🤡

Anyway now's time for bed, and I don't have brainpower to check this

[b3ll]

I don't understand yet what the values mean, but from the docs it looks like you're using scipy correctly :P

Thank you for all your help though! I really appreciate it

[jpcima]

So, I've returned to this matter and figured out the filter.
It's formulation known as Partial Fraction Expansion expressed as:

H(s) = sum for i=[1:M] : (Ri / (s - Pi))
having s=jω=j2πf

and added in the last commit. By plotting the response, I found a minor glitch which occurred on the input filter in the neighbourhood of 100kHz, and I fixed the coefficient.

Transfer function program: https://gist.github.com/jpcima/38c464509213620841385a3a09f96547

And a logarithmic plot:

[b3ll]

oh that's awesome! That's super cool, I'll def read up on that… you've gotta see the plugin I've built with this, it sounds really great!

[jpcima]

I post here the Octave script creating custom filters.
It's a 5th-order Butterworth 8kHz for this case.

Probably, the Juno one didn't use Butterworth, it used a steeper model for its cutoff.
The TCA-350-Y chip from the Solina is known to have Butterworth output from its datasheet.

For other filter design kinds:
https://fr.mathworks.com/help/signal/ug/comparison-of-analog-iir-lowpass-filters.html

The script output is R and P arrays, with M being filter's order.

format longG;

cutoff=8000;

[b,a]=butter(5,2*pi*cutoff,'low','s');
%[h,w]=freqs(b,a);
w=logspace(1,5);
h=freqs(b,a,w);
dB = mag2db(abs(h));

plot(w/(2*pi),dB);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
grid on;
grid minor;

[r,p,k]=residue(b,a)

[b3ll]

Oh this is awesome, yeah the Juno had a -24/db lowpass filter, which seems right. Happy to finally understand where those magic numbers in the table came from hah