Porting moog ladder filters to audio objects?

Yes, of course. I mistakenly thought we were at the final form or getting really close. Maybe not?

As you said, it would be good to get some more feedback before committing to it. I will keep looking at how to get self-oscilation harmonics that are closer to the analog case, including testing other models, but I have a suspicion that it might take fairly high oversampling rates to achieve that.
 
With max resonance you can clearly hear that analog filter has some "gritty" quality to it (2:52 of your video), while digital is sterile and lacks this effect completely.

I had assumed the digital implementation was not going to be able to emulate the analog filter all that well at full resonance. But I went back for another listen to this clip, and compared the sections at 2:40 (analog) and 3:20 (teensy). At those times they're both sweeping similar low to mid frequencies at full resonance, and they don't actually sound all that different to me. I would expect the difference to increase as you go to high cutoff.
 
Below is today's version, v1.04, which I am hoping will be close to final apart from perhaps some optimisation. I decided that what it may need is simply more resonance, so like a real moog this one allows resonance to be pushed all the way to 1.8, well beyond first sustained oscillation at about 1.0. To achieve that without significant aliasing required oversampling to be increased to 4x, but the good news is that also extends the max_frequency up to 18.7kHz, and CPU usage is still just under 5%.
 
Code:
/* Audio Library for Teensy, Ladder Filter
 * Copyright (c) 2021, Richard van Hoesel
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to deal
 * in the Software without restriction, including without limitation the rights
 * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
 * copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice, development funding notice, and this permission
 * notice shall be included in all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 * THE SOFTWARE.
 */

//--------------------------------------------------------------
// Based on Huovilainen New Moog (HNM) model as per CMJ jun 2006
// Implemented as Teensy37 Audio Library compatible object
// Richard van Hoesel, v. 1.04, Feb. 15 2021
// v1.4 FC extended to 18.7kHz, max res to 1.8, 4x oversampling, 
//      and a minor Q-tuning adjustment 
// v.1.03 adds 2x oversampling, extended resonance,  
//        and exposes parameters input_drive and passband_gain   
// v.1.02 now includes both cutoff and resonance "CV" modulation inputs
// please retain this header if you use this code.  
//--------------------------------------------------------------

// [URL]https://forum.pjrc.com/threads/60488?p=269755&viewfull=1#post269755[/URL]
// [URL]https://forum.pjrc.com/threads/60488?p=269609&viewfull=1#post269609[/URL]

#include <Arduino.h>
#include "filter_ladder.h"
#include <math.h>
#include <stdint.h>
#define MOOG_PI ((float)3.14159265358979323846264338327950288)

//#define  osTimes 1
//#define MAX_RESONANCE ((float)1.1)                                       
//#define MAX_FREQUENCY ((float)(AUDIO_SAMPLE_RATE_EXACT * 0.249f))        

//#define osTimes 2
//#define MAX_RESONANCE ((float)1.2)                                       
//#define MAX_FREQUENCY ((float)(AUDIO_SAMPLE_RATE_EXACT * 0.38f))  

#define osTimes 4
#define MAX_RESONANCE ((float)1.8)                                       
#define MAX_FREQUENCY ((float)(AUDIO_SAMPLE_RATE_EXACT * 0.425f)) 
//#define lfq 0.25

float AudioFilterLadder::LPF(float s, int i)
{
    float ft = s * (1.0f/1.3f) + (0.3f/1.3f) * z0[i] - z1[i];
    ft = ft * alpha + z1[i];
    z1[i] = ft;
    z0[i] = s;
    return ft;
}

void AudioFilterLadder::resonance(float res)
{
    // maps resonance = 0->1 to K = 0 -> 4
    if (res > MAX_RESONANCE) {
        res = MAX_RESONANCE;
    } else if (res < 0.0f) {
        res = 0.0f;
    }
    K = 4.0f * res;
}

void AudioFilterLadder::frequency(float c)
{
    Fbase = c;
    compute_coeffs(c);
}

void AudioFilterLadder:: octaveControl(float octaves)
{
    if (octaves > 7.0f) {
        octaves = 7.0f;
    } else if (octaves < 0.0f) {
        octaves = 0.0f;
    }
    octaveScale = octaves / 32768.0f;
}

void AudioFilterLadder:: passband_gain(float passbandgain) 
{
  pbg = passbandgain;
  if (pbg>0.5) pbg = 0.5f;
  if (pbg<0) pbg = 0.0f;
  input_drive(host_overdrive);
}

void AudioFilterLadder::input_drive(float odrv) 
{
  host_overdrive = odrv;  
  if(host_overdrive >1)
  {
    if(host_overdrive>4) host_overdrive = 4.0f;
    overdrive = 1.0f + (host_overdrive - 1.0f) * (1.0 - pbg);   // max is 4 when pbg = 0, and 2.5 when pbg is 0.5    
  }
  else
  {
    overdrive = host_overdrive;
    if(overdrive < 0) overdrive = 0.0f;
  }
}

void AudioFilterLadder::compute_coeffs(float c)
{
    if (c > MAX_FREQUENCY) {
        c = MAX_FREQUENCY;
    } else if (c < 1.0f) {
        c = 1.0f;
    }
      #ifdef lfq
      if(c<500) lfkmod = 1.0f + (500.0f-c)*(1.0f/500.0f) * lfq; 
      #endif  
     float wc = c * (float)(2.0f * MOOG_PI / (osTimes * AUDIO_SAMPLE_RATE_EXACT));
    float wc2 = wc * wc;
    alpha = 0.9892f * wc - 0.4324f * wc2 + 0.1381f * wc * wc2 - 0.0202f * wc2 * wc2;
     //Qadjust = 1.0029f + 0.0526f * wc - 0.0926 * wc2 + 0.0218* wc * wc2;
    //Qadjust = 1.006f + 0.0536f * wc - 0.095 * wc2 ;     
    Qadjust = 1.006f + 0.0536f * wc - 0.095f * wc2 - 0.05f * wc2 * wc2;   
}

bool AudioFilterLadder::resonating()
{
    for (int i=0; i < 4; i++) {
        if (fabsf(z0[i]) > 0.0001f) return true;
        if (fabsf(z1[i]) > 0.0001f) return true;
    }
    return false;
}

static inline float fast_exp2f(float x)
{
    float i;
    float f = modff(x, &i);
    f *= 0.693147f / 256.0f;
    f += 1.0f;
    f *= f;
    f *= f;
    f *= f;
    f *= f;
    f *= f;
    f *= f;
    f *= f;
    f *= f;
    f = ldexpf(f, i);
    return f;
}

static inline float fast_tanh(float x)
{
    float x2 = x * x;
    return x * (27.0f + x2) / (27.0f + 9.0f * x2);
}

void AudioFilterLadder::update(void)
{
  audio_block_t *blocka, *blockb, *blockc;
  float Ktot;
  bool FCmodActive = true;
  bool QmodActive = true;

  blocka = receiveWritable(0);
  blockb = receiveReadOnly(1);
  blockc = receiveReadOnly(2);
  if (!blocka) {
    if (resonating()) {
      // When no data arrives but the filter is still
      // resonating, we must continue computing the filter
      // with zero input to sustain the resonance
      blocka = allocate();
    }
    if (!blocka) {
      if (blockb) release(blockb);
      if (blockc) release(blockc);
      return;
    }
    for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
      blocka->data[i] = 0;
    }
  }
  if (!blockb) {
    FCmodActive = false;
  }
  if (!blockc) {
    QmodActive = false;
    Ktot = K;
  }
  for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
    float input = blocka->data[i] * (1.0f/32768.0f) * overdrive;
    if (FCmodActive) {
      float FCmod = blockb->data[i] * octaveScale;
      float ftot = Fbase * fast_exp2f(FCmod);
      if (ftot > MAX_FREQUENCY) ftot = MAX_FREQUENCY;
      compute_coeffs(ftot);
    }
    if (QmodActive) {
      float Qmod = blockc->data[i] * (1.0f/32768.0f);
      Ktot =   K + 4.0f * Qmod;     
    }
    #ifdef lfq
    if(Ktot > MAX_RESONANCE * 4 * lfkmod) Ktot = MAX_RESONANCE * 4 * lfkmod;
    #else
    if(Ktot > MAX_RESONANCE * 4 ) Ktot = MAX_RESONANCE * 4;
    #endif
    if (Ktot < 0.0f) Ktot = 0.0f;
    float total = 0;
    for(int os = 0; os < osTimes; os++)
    {
      float u = input - (z1[3] - pbg * input) * Ktot; 
      u = fast_tanh(u);
      float stage1 = LPF(u, 0);
      float stage2 = LPF(stage1, 1);   
      float stage3 = LPF(stage2, 2);  
      float stage4 = LPF(stage3, 3);        
      total += stage4 * 1/osTimes;
    }            
    blocka->data[i] = total * 32767.0f ;    
  }
  transmit(blocka);
  release(blocka);
  if (blockb) release(blockb);
  if (blockc) release(blockc);
}
 
Last edited by a moderator:
I updated github with the latest.

https://github.com/PaulStoffregen/Audio/commit/bd60b95a55652af7e07528d3f3125bce3c80b921

Also moved the initialization of Ktot inside the loop, just in case it's needed for more experiments in dynamically updated Kmax.

I still don't really understand the point of Qadjust. Is this planned to be used in a future version? As nearly as I can tell this variable is never actually used anywhere to affect the sound in any way.
 
A question :
At the beginning of the loop, you do a normalization and multiply the samples whith 1/32768.
Is that needed?
There could be a chance for a little optimization...
At the end of the loop, you scale it back by multiplying it with 32767 (<- why not 32768?).

Just a question.
 
I updated github with the latest.

https://github.com/PaulStoffregen/Audio/commit/bd60b95a55652af7e07528d3f3125bce3c80b921

Also moved the initialization of Ktot inside the loop, just in case it's needed for more experiments in dynamically updated Kmax.

I still don't really understand the point of Qadjust. Is this planned to be used in a future version? As nearly as I can tell this variable is never actually used anywhere to affect the sound in any way.

Good catch Paul. I have overlooked adding Qadjust back into the update loop in the file I posted above. The first line in the filter loop should read:

float u = input - (z1[3] - pbg * input) * Ktot * Qadjust;

Sorry about that. I had removed it in an earlier version thinking it was not necessary, but it allows more refined control of the resonance.
 
@PaulStoffregen, spent a little time with the new code on git, it sounded different, after a bit of testing it appears the code you've loaded on git for v1.4 is different to the code at #129 in a material way.

I know you're optimising it but the git copy has an effect of the filter - I can switch between the git version and the #129 version and get a different result.

If you want to test it dial in max res (1.79) and 20Hz Cutoff. On the #129 version I get a 20Hz tone and on the git version I get a 80Hz tone.
 
A question :
At the beginning of the loop, you do a normalization and multiply the samples whith 1/32768.
Is that needed?
There could be a chance for a little optimization...
At the end of the loop, you scale it back by multiplying it with 32767 (<- why not 32768?).

Just a question.

(edited after I had my morning coffee :)) The tanh won't work if the input isn't scaled down appropriately. And the reason I used 32767 as output scaler was to to stay within the 16bit limit for the DAC, which can't code -32768. Is that not needed?
 
Last edited:
@PaulStoffregen, spent a little time with the new code on git, it sounded different, after a bit of testing it appears the code you've loaded on git for v1.4 is different to the code at #129 in a material way.

I know you're optimising it but the git copy has an effect of the filter - I can switch between the git version and the #129 version and get a different result.

If you want to test it dial in max res (1.79) and 20Hz Cutoff. On the #129 version I get a 20Hz tone and on the git version I get a 80Hz tone.


When I look at the Github code, the compute_coeffs routine has errors and redundancies. The 4x frequency error is because the line:
float wc = c * (float)(2.0f * MOOG_PI / AUDIO_SAMPLE_RATE_EXACT);
should read:
float wc = c * (float)(2.0f * MOOG_PI / (osTimes * AUDIO_SAMPLE_RATE_EXACT));

And Qadjust is calculated twice in the Github code, but we only want the second one.
 
I've fixed the frequency error (a silly oversight on my part while manually merging) and included the resonance adjustment (msg #132) on github.

https://github.com/PaulStoffregen/Audio/commit/c5621587396a507a4b3fb3347bb8aa855c151790

Hopefully we're getting close to the final form? I'm holding off considering optimizations. CPU usage is currently at 4.77% on Teensy 4.0 and 30% on Teensy 3.5. Teensy 3.2 and other boards without FPU have no hope of running any version of this filter.
 
Thanks Paul. I was trying various settings on the non-modulated inputs to decide on final default values for the .h file.
I'm inclined to set the default passband gain, determined by the variable pbg, to 0.05 (nb not 0.5) rather than 0.
The rest seems fine as is.

One thought I had is that we don't do a check to see if the output of the filter exceeds 16 bit limits. Is there a saturating multiply that can be used there rather than introducing more 'if' conditions?

Other than that, from my point of view this version seems more or less ready, unless someone hears something that still needs attention.

I'm still looking at other models too, but this one retains great stability no matter what I throw at it.
 
One other thought. The filter structure apparently also allows high-pass and band-pass outputs by combining different weighting of the four LPF stages. However, I haven't tested that yet. Is there any interest in that?

-- edit --
Never mind this suggestion. I just tried it out briefly and didn't care much for the sound of the resulting hpf and bpf versions. Fine at low resonance, but not at higher resonance.
 
Last edited:
One last minor adjustment. In compute_coeffs(float c), can we please increase the lower cutoff limit to 5Hz, rather than 1Hz?
Thanks,
Richard.
 
Last edited:
@rvh: I have not looked at this in detail. It just caught my eye when I briefly looked over it. If it doesn't work without scaling because of tanh, ok. It also doesn't matter if you multiply by 32767 or 32768.
 
@rvh: I have not looked at this in detail. It just caught my eye when I briefly looked over it. If it doesn't work without scaling because of tanh, ok. It also doesn't matter if you multiply by 32767 or 32768.

Thanks. How is it ensured that -32768 doesn't overflow the 16 bit integer range of the dac?
 
I coded up a simple 'midi synth' with bandlimited sawtooth oscillator running into this filter, and set the filter in self oscillation at Q of about 1.2, and cutoff tracking the keyboard notes. When I got to a few octaves above middle C I heard quite a lot of aliasing. But the problem went away when I turned down the oscillator (filter still resonating just as loudly), so I don't think the problem lies with the filter. Perhaps there's scope to improve the bandlimted oscillators further?
 
I'm not sure what you're asking or if I understand you question.
As long your sample is >=-1.0 and < 1.0 , a multiplication 32768 will do because it fits into 16 bits. 32767 does not make big difference. Don't know if it's worth to write long posts about it ;) Nobody will hear a diffenece, since it's only the volume.
If you're _not_ sure that it fits, there is a simple dsp "asm" saturation-instruction "ssat" that you can use. The cpu can do this in one cycle.
 
When playing the tuned resonant filter with a musical keyboard, you can produce fairly large sudden Fc changes that produce abrupt 'attacks'. Should we include a way of smoothing those changes using some kind or filter-cutoff portamento parameter? Presumably a simple first order LPF running at block-rate (rather than sample rate) would be all that's needed.
 
Last edited:
I put the code on github and make one small fix and adapted the names to be similar to the other filters.

https://github.com/PaulStoffregen/Audio

If anyone wants to give this a try, please grab the latest from github. If you're not familiar with using git, the simplest way is to just download the ZIP file and extract it in {Documents}/Arduino/libraries. Just keep in mind anything to put there overrides all the other locations, so remember to delete this test copy when you later want to use the audio library the Teensyduino installer puts into your copy of Arduino.

Since there aren't any test programs yet, I wrote this very simple program which drives the filter with bursts of a sawtooth waveform.

Code:
#include <Audio.h>

AudioSynthWaveform       waveform1;
AudioFilterLadder        filter1;
AudioOutputMQS           mqs1;
AudioConnection          patchCord1(waveform1, 0, filter1, 0);
AudioConnection          patchCord2(filter1, 0, mqs1, 0);
AudioConnection          patchCord3(filter1, 0, mqs1, 1);

void setup() {
  AudioMemory(40);
  filter1.resonance(1.0);
  filter1.frequency(440);
  waveform1.begin(WAVEFORM_SAWTOOTH);
  waveform1.frequency(100);
}

void loop() {
  waveform1.amplitude(0.9);
  delay(500);
  waveform1.amplitude(0);
  delay(1500);
}

I'm not a synth guy, so I really don't know what this filter is supposed to sound like. But here's the output I see with my oscilloscope at the end of the sawtooth burst.

View attachment 23655

I'm not 100% confident it's really correct. As the test repeats every 2 seconds, the amount of 440 Hz superimposed on the sawtooth seems to vary quite a lot. Maybe later I'll try to capture a time-lapse video of my scope screen....

Really, I'm depending on everyone here with an ear for what the Moog sound is supposed to be to say whether they're happy with this, or if more work is needed before releasing it.

Also, still need a clear confirmation on the MIT license before this can be released.

Probably is a good idea to check on that to cover yourself but the Moog design is over 50 years old now and should be well and truly out of copywrite. The circuit diagrams are out there for free download and there are a lot of copies of the design being made, some good, some bad.

As for some further hints on the sound this clip has a bit of direction to it I think especially demonstrating how the panel controls effect the filter operation.
https://www.youtube.com/watch?v=dVgIf71uWB4
 
I noticed that you are using copy of the same input sample when oversampling and averaging output. You may consider using linear interpolation on input and dropping on output.

@Tomas, thanks for your ongoing comments in this thread. I am reviewing some of them again today, including this one. Interestingly, when I tried dropping the output averaging when down-sampling on the latest 4x oversampling version, aliasing was much worse. I may actually re-instate input interpolation too, just in case it's contributing to the aliasing I hear when I play bandlimited sawtooth waveforms with the filter self-resonating at the matched frequency. In my comment above I suggested the sawtooth itself may be the cause of the aliasing because wirh reduced input level (that doesn't affect the self resonance) the aliasing disappears.

Also, noting your comment about tanh not being bound beyond 4.5, I will review your suggestion to (re)instate +/-3 input checks. My reason is that after more listening, I find I don't really much care for the sound of the filter when it's over-driven. However, I do find that the output level (and hence feedback) of this particular model is remarkably constant and never seems to 'clip', even when over-driven by up to a factor of 4. So maybe my limit of 4x already ensures tanh remains within it's safe range.

On a somewhat related note, I had a short communication with Stefano D'Angelo about these digital models and he recommends using small signal levels to better approximate what goes on in a real moog because apparently the analog signals usually are too. I may experiment with scaling down the input further, an adding the gain back at the output for that purpose.
--- edit --- I realized there is a problem with this because the self resonance is independent of input scaling, so you can't scale the output up. We can already just scale down the input using the 'input_drive' parameter, and it may be worth experimentng a bit more with that to determine the default. Setting it to less than 1 also makes the self resonance relatively more prominent for those who are looking for that.

He also kindly pointed me to some matlab code implementing his 2014 IEEE published model that uses a delay-free loop method. Unfortunately, I have minimal matlab experience, so it may take me a while to figure out how to translate it to C-code. If anyone here is a matlab wiz and would be willing to help out, please let me know!
 
Last edited:
If you did oversampling "by the book" you would need polyphase FIR interpolator at the input and polyphase FIR decimator filter at the output. These are costly. Not doing things by the book results in aliasing. FWIW D'Angelo filter (Arturia) is running at 176kHz rate https://www.kvraudio.com/forum/viewtopic.php?p=6601090&sid=f6b9cd37e87e927ff9025cfe4b24bb49#p6601090

Do you know of any c++ code I could use to implement a polyphase FIR filter (ideally ready to go with coefficients for 4x oversampling) to see whether it offers any significant advantage in this application?
Although costly, perhaps on a teensy 4 they're feasible.
 
The CMSIS has interpolation:
https://arm-software.github.io/CMSI...olate.html#ga62d0757246bbbe5bfaf67b9a95c16c51
It's part of the Teensy core, so nothing to install.

As said, it needs a low pass filter. You can calculate the coefficients with a filter-calculation tool, matlab or whatever.
You can use a const array to store them, then.
Edit: Iowa Hills filter designer: http://iowahills.com/

I have a project where I use it, if you want to see an example. (8x oversampling)
It's meant to be a library some day, but for now it's just a program which includes the Audio-Library part and tabs in Arduino.

In the code the calulation of the coeffs happens in begin() - I used code by WMXZ and DD4WH that generates them.
Note, that not the interpolation is costly, it's just the fact that you have to process 4x (8x whatever) more samples. So it will need a bit more than i.e. 4x the time. Don't forget to scale the samples after that (or at the end of thge loop, when you store the new calculated sample) by the interpolation factor.
As my code uses 2 channels, for stereo, it even does interpolation by 8 two times - still not a big problem for the Teensy 4.

The CMISIS has functions for decimation, too (not used in my project, but you will need it).
 
Last edited:
@Frank, thanks very much. I'll take a look at your example.

-edit - Not really being a programmer, and after having had a quick look, it might take me a while to sort through your example. So if by any chance you'd be willing to show me just the lines I need to upsample a single 128pt audio-object block (by a factor of 4) that would be great.
 
Last edited:
Back
Top