Porting moog ladder filters to audio objects?

Below is version 1.5 for people to try out. By default it uses a 32-tap polyphase FIR interpolator/decimator now. That adds around 1.5% to the CPU load (per filter) compared to simple linear interpolation. You can switch between linear and polyphase FIR interpolation by setting the filter parameter "interpMethod" to LINEAR or FIR_POLY respectively.
Let me know what you think.
Cheers,
Richard.
 
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.
 */

//-----------------------------------------------------------
// Huovilainen New Moog (HNM) model as per CMJ jun 2006
// Implemented as Teensy Audio Library compatible object
// Richard van Hoesel, Feb. 21 2021
// v1.5 adds polyphase FIR or Linear interpolation 
// v1.4 FC extended to 18.7kHz, max res to 1.8, 4x oversampling,
//      and a minor Q-tuning adjustment
// v.1.03 adds 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 "arm_math.h"
#include <stdint.h>
#define MOOG_PI ((float)3.14159265358979323846264338327950288)    

#define MAX_RESONANCE ((float)1.8)                                       
#define MAX_FREQUENCY ((float)(AUDIO_SAMPLE_RATE_EXACT * 0.425f))     


void AudioFilterLadder::initpoly()
{
  if (arm_fir_interpolate_init_f32(&interpolation, INTERPOLATION, interpolation_taps, interpolation_coeffs, interpolation_state, NUM_SAMPLES))
  {
    Serial.println("Init of interpolation failed");
    while (1);
  }
  if (arm_fir_decimate_init_f32(&decimation, interpolation_taps, INTERPOLATION, interpolation_coeffs, decimation_state, I_NUM_SAMPLES))
  {
    Serial.println("Init of decimation failed");
    while (1);
  } 
}  

void AudioFilterLadder::interpMethod(int imethod)
{
  if(imethod==FIR_POLY)      
    polyOn=true;    
  else
    polyOn=false;
}

float AudioFilterLadder::LPF(float s, int i)
{
  float ft = s * float(1.0f/1.3f) + float(0.3f/1.3f) * z0[i] - z1[i];
  ft = ft * (1 * 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 < 5.0f) {
    c = 5.0f;
  }
  float wc = c * (float)(2.0f * MOOG_PI / (INTERPOLATION * 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.095f * wc2 - 0.05f * wc2 * wc2;     // revised hfQ (rvh - feb 14 2021)
}

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)
{
  if(x>3) return(1.0f);
  if(x<-3) return(-1.0f);
  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;
  }
  if(firstpoly) {
      initpoly();
      firstpoly=false; 
  }
  
  if(polyOn ==true)
  {
    /*----------------------- upsample -------------------------*/  
      float blockOS[I_NUM_SAMPLES], blockIn[AUDIO_BLOCK_SAMPLES], blockOutOS[I_NUM_SAMPLES], blockOut[AUDIO_BLOCK_SAMPLES];;
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++)
        blockIn[i] = blocka->data[i]  * overdrive *  float(INTERPOLATION /32768.0f);
      arm_fir_interpolate_f32(&interpolation, blockIn, blockOS, NUM_SAMPLES);     
    
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {      
        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;     
        }
        if(Ktot > MAX_RESONANCE * 4 ) Ktot = MAX_RESONANCE * 4;
          if (Ktot < 0.0f) Ktot = 0.0f;  
            for(int os = 0; os < INTERPOLATION; os++)  
            {      
              float input = blockOS[i*4 + os];   
              float u = input - (z1[3] - pbg * input) * Ktot * Qadjust; 
              u = fast_tanh(u);
              float stage1 = LPF(u, 0);
              float stage2 = LPF(stage1, 1);   
              float stage3 = LPF(stage2, 2);  
              float stage4 = LPF(stage3, 3);            
              blockOutOS[i*4 + os] = stage4 ;
            }      
      }
      arm_fir_decimate_f32(&decimation, blockOutOS, blockOut, I_NUM_SAMPLES);  
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++)
        blocka->data[i] = blockOut[i] * float(32768.0); 
  }
  else  // linear interpolation
  {
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
        float input = blocka->data[i]  * overdrive * (1.0f/32768.0f);
        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;
        float interp = 0;
        for(int os = 0; os < INTERPOLATION; os++)
        {      
          float u = (interp * oldinput + (1.0f-interp)*input)- (z1[3] - pbg * input) * Ktot * Qadjust; 
          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 * float(1.0/INTERPOLATION);         
          interp += float(1.0f/INTERPOLATION);    
        }
        oldinput=input;
        blocka->data[i] = total * 32768.0f ;    
      }
  }     
  transmit(blocka);
  release(blocka);
  if (blockb) release(blockb);
  if (blockc) release(blockc);
}
 
Last edited:
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.
 */

//-----------------------------------------------------------
// Huovilainen New Moog (HNM) model as per CMJ jun 2006
// Implemented as Teensy Audio Library compatible object
// Richard van Hoesel, Feb. 21 2021
// v1.5 adds polyphase FIR or Linear interpolation
// v1.4 FC extended to 18.7kHz, max res to 1.8, 4x oversampling,
//      and a minor Q-tuning adjustment
// v.1.03 adds 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=269609&viewfull=1#post269609[/URL]

#ifndef filter_ladder_h_
#define filter_ladder_h_

#define LINEAR 0
#define FIR_POLY 1
#include "Arduino.h"
#include "AudioStream.h"
#include "arm_math.h"
#define LINEAR 0
#define FIR_POLY 1
#define INTERPOLATION  4
#define INTERPOLATIONscaler  4
#define I_SAMPLERATE  AUDIO_SAMPLE_RATE_EXACT * INTERPOLATION
#define NUM_SAMPLES  128
#define I_NUM_SAMPLES  NUM_SAMPLES * INTERPOLATION
#define interpolation_taps 32
static float interpolation_coeffs[interpolation_taps] = {
444.8708144065367380E-6, 0.002132195798620151, 0.004238318820060649, 0.005517336264360800, 0.004277228394758379,-780.0288538232418890E-6,-0.009468419595216885,-0.019306296307004593,
-0.025555801376824804,-0.022415559173490765,-0.005152700802329096, 0.027575475722578378, 0.072230543044204509, 0.120566259375849666, 0.161676773663769008, 0.185322994251944373,
 0.185322994251944373, 0.161676773663769008, 0.120566259375849666, 0.072230543044204509, 0.027575475722578378,-0.005152700802329096,-0.022415559173490765,-0.025555801376824804,
-0.019306296307004593,-0.009468419595216885,-780.0288538232418890E-6, 0.004277228394758379, 0.005517336264360800, 0.004238318820060649, 0.002132195798620151, 444.8708144065367380E-6
};

class AudioFilterLadder: public AudioStream
{
public:
  AudioFilterLadder() : AudioStream(3, inputQueueArray) {};
  void frequency(float FC);
  void resonance(float reson);
  void octaveControl(float octaves);
  void passband_gain(float passbandgain); 
  void input_drive(float drv);    
  void interpMethod(int im);   
  void  initpoly(); 
  virtual void update(void);
  
private:
  float interpolation_state[(NUM_SAMPLES-1) + interpolation_taps / INTERPOLATION];
  arm_fir_interpolate_instance_f32 interpolation;
  float decimation_state[(I_NUM_SAMPLES-1) + interpolation_taps];
  arm_fir_decimate_instance_f32 decimation;
  float LPF(float s, int i);
  void  compute_coeffs(float fc);
  bool  resonating();  
  bool  firstpoly=true;
  bool  polyOn = true;            // linear default  
  float alpha = 1.0;
  float beta[4] = {0.0, 0.0, 0.0, 0.0};
  float z0[4] = {0.0, 0.0, 0.0, 0.0};
  float z1[4] = {0.0, 0.0, 0.0, 0.0};
  float K = 1.0f;
  float Fbase = 1000;
  float Qadjust = 1.0f;
  float octaveScale = 1.0f/32768.0f;
  float pbg = 0.5f;
  float overdrive = 0.5f;
  float host_overdrive = 1.0f;
  float oldinput = 0; 
  audio_block_t *inputQueueArray[3];
};
#endif
 
Last edited:
Thanks Mark, the trade-off here off course is stopband suppression versus potentially audible HF rolloff in the passband within the constraint of no more than 32 points.
I have tried a few more 32-point designs and have settled on this for now:

View attachment 23794

You might get a better response with the Remez algorithm rather than the window method?
 
You might get a better response with the Remez algorithm rather than the window method?

I tried a few designs using Parks McClellan, and decided in the end to let the filter be just marginally longer. This one is 36 taps, derived from a 40-tap specification that had very small end coefficients.

FIR36.jpg

I think I'll leave it at this for now. It sounds pretty good to me. Here is the updated .h file. If you prefer the previous 32-tap filter, you can uncomment those coefficients instead.

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.
 */

//-----------------------------------------------------------
// Huovilainen New Moog (HNM) model as per CMJ jun 2006
// Implemented as Teensy Audio Library compatible object
// Richard van Hoesel, Feb. 21 2021
// v1.5 adds polyphase FIR or Linear interpolation
// v1.4 FC extended to 18.7kHz, max res to 1.8, 4x oversampling,
//      and a minor Q-tuning adjustment
// v.1.03 adds 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.
//-----------------------------------------------------------

// https://forum.pjrc.com/threads/60488?p=269609&viewfull=1#post269609

#ifndef filter_ladder_h_
#define filter_ladder_h_

#define LINEAR 0
#define FIR_POLY 1
#include "Arduino.h"
#include "AudioStream.h"
#include "arm_math.h"
#define LINEAR 0
#define FIR_POLY 1
#define INTERPOLATION  4
#define INTERPOLATIONscaler  4
#define I_SAMPLERATE  AUDIO_SAMPLE_RATE_EXACT * INTERPOLATION
#define NUM_SAMPLES  128
#define I_NUM_SAMPLES  NUM_SAMPLES * INTERPOLATION

/*
#define interpolation_taps 32
static float interpolation_coeffs[interpolation_taps] = {
444.8708144065367380E-6, 0.002132195798620151, 0.004238318820060649, 0.005517336264360800, 0.004277228394758379,-780.0288538232418890E-6,-0.009468419595216885,-0.019306296307004593,
-0.025555801376824804,-0.022415559173490765,-0.005152700802329096, 0.027575475722578378, 0.072230543044204509, 0.120566259375849666, 0.161676773663769008, 0.185322994251944373,
 0.185322994251944373, 0.161676773663769008, 0.120566259375849666, 0.072230543044204509, 0.027575475722578378,-0.005152700802329096,-0.022415559173490765,-0.025555801376824804,
-0.019306296307004593,-0.009468419595216885,-780.0288538232418890E-6, 0.004277228394758379, 0.005517336264360800, 0.004238318820060649, 0.002132195798620151, 444.8708144065367380E-6
};
*/

#define interpolation_taps 36
static float interpolation_coeffs[interpolation_taps] = {
-14.30851541590154240E-6,  0.001348560352009071, 0.004029285548698377, 0.007644563345368599, 0.010936856250494802, 0.011982063548666887, 0.008882946305001046, 826.6598116471556070E-6,
-0.011008071930708746,-0.023014151355548934,-0.029736402750934567,-0.025405787911977455,-0.006012006772274640, 0.028729626071574525, 0.074466890595619062, 0.122757573409695370,
 0.163145421379242955, 0.186152844567746417, 0.186152844567746417, 0.163145421379242955, 0.122757573409695370, 0.074466890595619062, 0.028729626071574525,-0.006012006772274640,
-0.025405787911977455,-0.029736402750934567,-0.023014151355548934,-0.011008071930708746, 826.6598116471556070E-6, 0.008882946305001046, 0.011982063548666887, 0.010936856250494802,
 0.007644563345368599, 0.004029285548698377, 0.001348560352009071,-14.30851541590154240E-6
};


class AudioFilterLadder: public AudioStream
{
public:
  AudioFilterLadder() : AudioStream(3, inputQueueArray) {};
  void frequency(float FC);
  void resonance(float reson);
  void octaveControl(float octaves);
  void passband_gain(float passbandgain); 
  void input_drive(float drv);    
  void interpMethod(int im);   
  void  initpoly(); 
  virtual void update(void);
  
private:
  float interpolation_state[(NUM_SAMPLES-1) + interpolation_taps / INTERPOLATION];
  arm_fir_interpolate_instance_f32 interpolation;
  float decimation_state[(I_NUM_SAMPLES-1) + interpolation_taps];
  arm_fir_decimate_instance_f32 decimation;
  float LPF(float s, int i);
  void  compute_coeffs(float fc);
  bool  resonating();  
  bool  firstpoly=true;
  bool  polyOn = true;            // linear default  
  float alpha = 1.0;
  float beta[4] = {0.0, 0.0, 0.0, 0.0};
  float z0[4] = {0.0, 0.0, 0.0, 0.0};
  float z1[4] = {0.0, 0.0, 0.0, 0.0};
  float K = 1.0f;
  float Fbase = 1000;
  float Qadjust = 1.0f;
  float octaveScale = 1.0f/32768.0f;
  float pbg = 0.5f;
  float overdrive = 0.5f;
  float host_overdrive = 1.0f;
  float oldinput = 0; 
  audio_block_t *inputQueueArray[3];
};
#endif
 
Hi Richard

Had a quick go with the FIR version (your 1.5 + 36 taps). I couldn't get to to run initially which I now know is because I was running with a smaller than standard AUDIO_BLOCK_SAMPLES. This seems to impacts both the linear and FIR configs. If I switch to the default AUDIO_BLOCK_SAMPLES = 128 it runs fine.

Sound really good - I was looking at latency as want to use as part of an effect rather than for synth.

With 128 sample blocks and a simple scheme of just the filter and an input mixer I get about 6.6ms latency from line in to line out. This is a bit long for use in effects which is what I'm looking to do.

If I change your code so that NUM_SAMPLES = AUDIO_BLOCK_SAMPLES (#define NUM_SAMPLES AUDIO_BLOCK_SAMPLES) in the header I can reduce the block size.

At 16 sample blocks latency in the same set up reduces to around 1.6ms, which is very useable. Processor usage is about the same at 6%.

Does changing NUM_SAMPLES impact on the calculation you are running - I haven't spent a long time listening but initial view is it sounds the same. If it doesn't I'd suggest changing to match AUDIO_BLOCK_SAMPLES for general use.

Thanks for your work on this, will be great to have a nice sounding LPF. Cheers, Paul
 
bro @rvh you are absolutely killing it! and shout out to OG paul for pulling/merging so fast. i really can't wait to try this one myself and give you guys some more synth person feedback, and while most of my synth building is still based on the www.axoloti.com platform for now, i am eagerly watching the teensy audio lib development and the patcher GUI, to maybe one day make a real project based on this.
 
Before merging this to the audio library, I want to see these new constants moved either inside the class definition (eg, at static const types) or to the cpp file. I can do it, if you would like to go with whatever style I choose. Or you can update in the way you prefer.

But we just can't have a lot of pretty generic names defined with global scope, where they become a long-term liability for conflicts with other stuff in the audio library and possibly other libraries that user might wish to include in their project. These constants need to be confined to the cpp file, or declared as const stuff inside the class, so they can't conflict with other libraries.
 
Hi Richard

Had a quick go with the FIR version (your 1.5 + 36 taps). I couldn't get to to run initially which I now know is because I was running with a smaller than standard AUDIO_BLOCK_SAMPLES. This seems to impacts both the linear and FIR configs. If I switch to the default AUDIO_BLOCK_SAMPLES = 128 it runs fine.

Sound really good - I was looking at latency as want to use as part of an effect rather than for synth.

With 128 sample blocks and a simple scheme of just the filter and an input mixer I get about 6.6ms latency from line in to line out. This is a bit long for use in effects which is what I'm looking to do.

If I change your code so that NUM_SAMPLES = AUDIO_BLOCK_SAMPLES (#define NUM_SAMPLES AUDIO_BLOCK_SAMPLES) in the header I can reduce the block size.

At 16 sample blocks latency in the same set up reduces to around 1.6ms, which is very useable. Processor usage is about the same at 6%.

Does changing NUM_SAMPLES impact on the calculation you are running - I haven't spent a long time listening but initial view is it sounds the same. If it doesn't I'd suggest changing to match AUDIO_BLOCK_SAMPLES for general use.

Thanks for your work on this, will be great to have a nice sounding LPF. Cheers, Paul

Hi Paul,

Yes, your suggested change is a good idea. It should make no difference to the processing as long as the blocks contain enough samples for the interpolation filter. At the moment that would set a lower limit of 9 samples because the 4x OS LPF has 36 taps. If you switch to linear interpolation, you could go smaller.

@PaulStoffregen, if you can't see any reason not to, can you please change to the NUM_SAMPLES define in the header file to #define NUM_SAMPLES AUDIO_BLOCK_SAMPLES as per Paul's is suggestion? And maybe in the documentation include a warning about not using very small block sizes in the POLY_FIR interpolation mode. Thanks.

I'm glad to hear you're enjoying the sound of the filter so far. I plan to keep working on it, as well as possible alternatives, but they will take some time.

Cheers,
Richard
 
Before merging this to the audio library, I want to see these new constants moved either inside the class definition (eg, at static const types) or to the cpp file. I can do it, if you would like to go with whatever style I choose. Or you can update in the way you prefer.

But we just can't have a lot of pretty generic names defined with global scope, where they become a long-term liability for conflicts with other stuff in the audio library and possibly other libraries that user might wish to include in their project. These constants need to be confined to the cpp file, or declared as const stuff inside the class, so they can't conflict with other libraries.

Hi Paul,
I'm totally fine with you making alterations as you see fit to make the code more compliant with the rest of the framework.
Thanks for the offer.
Cheers,
Richard
 
I've merged the latest code, moved the constants into the class private members, and taken a few small liberties with the public API to make it similar to the rest of the audio library.

Documentation still needs descriptions filled in.

https://www.pjrc.com/teensy/gui/index2.html?info=AudioFilterLadder

If anyone wants to recommend words to include in the 3 "TODO: add info here..." places, now's the time. Please try to write from the perspective of advising an Arduino user who's enthusiastic about getting the most from the filter, but really just wants to know how to best use these functions to achieve their synth goals. Try to avoid unnecessary discussion of the internal operation on the filter, or at least keep it simple and focused on an end user's desires. If they really want to learn how the filter works, they can read the code or this forum thread (which is linked from the source code).
 
I've merged the latest code, moved the constants into the class private members, and taken a few small liberties with the public API to make it similar to the rest of the audio library.

Documentation still needs descriptions filled in.

https://www.pjrc.com/teensy/gui/index2.html?info=AudioFilterLadder

If anyone wants to recommend words to include in the 3 "TODO: add info here..." places, now's the time. Please try to write from the perspective of advising an Arduino user who's enthusiastic about getting the most from the filter, but really just wants to know how to best use these functions to achieve their synth goals. Try to avoid unnecessary discussion of the internal operation on the filter, or at least keep it simple and focused on an end user's desires. If they really want to learn how the filter works, they can read the code or this forum thread (which is linked from the source code).

Thanks for doing that Paul. I agree that it needs to be simplified. We probably have more parameters than we need at the moment. Having used the filter for a bit, the input_drive doesn't sound very good to me when you use values larger than 1, and I never use linear interpolation because the poly_fir seems to run fast enough for most applications. So I would suggest leaving those out of the description (but leave them in the code?) and just describe the key controls: cutoff & resonance (and their modulators), octavecontrol, and optionally passband_gain. You presumably already have a description for most of those in the state-variable filter.

The function of the passband_gain (0-0.5) is to partly restore the level of the passband signal, which is increasingly attenuated in ladder filters as resonance goes up. If you want typical analog ladder filter behavior, leave it at 0 (or close to it).
 
For those of you who are interested in this sort of thing, here is a comparison of the current filter based on Huovilainen's simplified model with his full model. Input is a very low level 70Hz saw, and both models are running at strong self-oscillation with the filter tuned at about 550Hz. You can see the simplified model shows higher level harmonic distortion of the 550Hz self oscillation frequency compared to the full model (eg around 3x the fundamental at around 1600Hz), presumably due to the use of a single tanh function rather than one for each LPF stage. In agreement with the statements in the simplified-model paper, it's a reasonably good approximation to the full model at lower resonance, but less so at higher resonance, where it starts to sound different.

edit - I should mention this was with the filter set to an extreme resonance of 1.8 to get the self-oscillation level to better match that of the full model level. You can audibly reduce the sine distortion by winding the resonance back closer to 1.0, but the self-oscillation becomes significantly less prominent. Personally, I think I prefer a resonance limit of around 1.2 with this filter. There is also significantly increased aliasing when you push the filter to very high frequencies at these extreme resonance settings. I'd like to hear from anyone who has tried the filter, whether they agree that 1.8 is too high, and whether we should perhaps reduce the maximum available value of that parameter.
 

Attachments

  • Comparison spectra.jpg
    Comparison spectra.jpg
    199.1 KB · Views: 101
Last edited:
On a related note, what are people's views on CPU demands versus performance with these filters? Currently the simplified model version uses between 6-7% at 4x oversampling. I think I can get a 4x oversampling version of the full model running at about double that, and if I run the full model at 2x oversampling, around 9%. I'm fairly sure most people would find the sound of the full model closer to an analog ladder filter. Is that what we are aiming for? Do we want more than one filter option, and if so, how do we integrate them into the library?
 
On a related note, what are people's views on CPU demands versus performance with these filters?
this might be the obvious answer you weren't looking for but as a user i'd ideally want to have both options. there are valid use cases for both, and for example a lot of commercial VST audio plugins also offer several quality settings to find the right balance between CPU load and sound quality.
there will be even faster teensys in the (not so far) future, and some people might even rather add multiple MCUs (for more CPU power) to a project to not compromise on sound qualityl

not sure how that would work out in terms of efficiency and bloating the code but the most convenient solution would be to have ONE audio library object that includes a "quality" or "model" switch.
 
this might be the obvious answer you weren't looking for but as a user i'd ideally want to have both options. there are valid use cases for both, and for example a lot of commercial VST audio plugins also offer several quality settings to find the right balance between CPU load and sound quality.
there will be even faster teensys in the (not so far) future, and some people might even rather add multiple MCUs (for more CPU power) to a project to not compromise on sound qualityl

not sure how that would work out in terms of efficiency and bloating the code but the most convenient solution would be to have ONE audio library object that includes a "quality" or "model" switch.

Thanks for the suggestions. I've managed to get a 2x over-sampled version of the new model (which I'll describe in the next post) sounding as good as a 4x over-sampled version for all intents and purposes, so I don't plan to release the latter anymore.

I can understand why you might like both the old model and the new one in one object, but there are reasons not to as well. The new model is not a higher quality version of the old one, it's a substantially different method. They sound quite different especially at higher resonance, and they don't quite use the same parameters either.
 
Here is the new full-model filter, which for now I've called AudioFilterLadder2. It's based on Huovilainen's full model, presented at DAFX 2004. Below are the filter_ladder2.cpp and .h files if you would like to try it out.

It's been a very long time since I owned a minimoog, but this filter to me sounds clearly more like what I remember than does the simplified model version that we currently have in the Audio library. However, I have had to make some parameter decisions that aren't described in the paper (as is common). It's set up so that at full resonance the filter can produce those obnoxiously loud resonance sweeps that can be used to good effect when creating electro-percussion a-la Kraftwerk etc. To achieve that from typical 0-1 positive envelopes, I've also extended the octave control range to allow up to 10 octaves, whereas the limit with previous filters in the lib has been 7.

-- edit -- I've now attached a zip file containing an mp3 demo of the new filter generating this kind of sound

If you prefer more subdued self-resonance levels, you can do so by specifying a resonance just slightly shy of 1.0.

Note that in the paper this model has been set up expecting resonance values strictly between 0 and 1, which is perhaps why Ive heard it said elsewhere that this model isn't capable of self-oscillation. In my implementation I expand that to allow a max value of 1.015, and have tuned the filter with that exact range in mind. I don't recommend changing either that value (max res) or the max cutoff freq that's defined in these files. Doing so will disrupt the tuning, potentially lead to instability, and may cause much more aliasing. You have been warned :)

It works the same as the previous filter in terms of signals/modulators, and the parameters are mostly similar (keeping in mind the extended octaveControl range)

However, there's no .passbandGain control, and although the .inputDrive parameter is still in place, it's not really a useful and is limited to 1-2. Even pushing it to 10 hasn't done much in the way of interesting sound changes imo so I may take it out in future.

There is also a new parameter .portamento that is a first order LPF coefficient that allows gradual rather than immediate filter cutoff changes whenever ithe .frequency parameter is changed. A value of 0 means instantaneous changes (i.e. the way filters have behaved to date), whereas a value of 1 means the filter stops changing altogether. Slow glides result from values close to, but not equal to 1.






Here's the filter_ladder2.cpp file:

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.
 */

//-----------------------------------------------------------
// AudioFilterLadder2 for Teensy Audio Library
// Based on Huovilainen's full model, presented at DAFX 2004
// Richard van Hoesel, March. 4, 2021
// v.1.0 Uses 2x oversampling. Host resonance range is 0-1.0. 
// Self oscillation at lower levels is available by specifying 
// values just below 1.0   
//
// please retain this header if you use this code.
//-----------------------------------------------------------

#include <Arduino.h>
#include "filter_ladder2.h"
#include "arm_math.h"
#include <math.h>
#include <stdint.h>

#define MAX_RESONANCE ((float)1.015)
#define MAX_FREQUENCY ((float)(AUDIO_SAMPLE_RATE_EXACT * 0.425f))
#define fI_NUM_SAMPLES2x  AUDIO_BLOCK_SAMPLES * INTERPOLATION2x

//=== 2 xOS, 32taps (88.2kHz)
static float AudioFilterLadder2 ::finterpolation_coeffs2x[AudioFilterLadder2::finterpolation_taps2x] = {    
-129.6800658538650740E-6, 0.001562843504973615, 0.004914813078250063, 0.007186231102125209, 0.002728024844356039,-0.007973931496234599,-0.013255882861166815,-841.1646993462468340E-6,
 0.022309967341058432, 0.027240071522569135,-0.007073857882406733,-0.056523975383091875,-0.055929860712812876, 0.042278471323346570, 0.207056768528768836, 0.335633514003638445,
 0.335633514003638445, 0.207056768528768836, 0.042278471323346570,-0.055929860712812876,-0.056523975383091875,-0.007073857882406733, 0.027240071522569135, 0.022309967341058432,
-841.1646993462468340E-6,-0.013255882861166815,-0.007973931496234599, 0.002728024844356039, 0.007186231102125209, 0.004914813078250063, 0.001562843504973615,-129.6800658538650740E-6
};

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;
}

void AudioFilterLadder2::inputDrive(float drive) 
{
  if(drive<1)
      overdrive = 1.0f;
  else 
  if(drive>2)
     overdrive =  2.0f;  
  else
    overdrive = drive;   
}

void AudioFilterLadder2::portamento(float pc)
{
  FcPorta = pc; 
}

void AudioFilterLadder2::compute_coeffs2x(float freq, float res) // tuned for 2x os
{
    float fc, freq2, freq3;          
    if (res > MAX_RESONANCE) res = MAX_RESONANCE ;
    if (res < 0.0f) res = 0.0f;    if (freq < 5) freq = 5;
    if (freq > MAX_FREQUENCY) freq = MAX_FREQUENCY;       
    fc = freq *  FC_SCALER2x;                             
    freq2 = freq*freq;
    freq3 = freq2 * freq;      
    float qa = 0.998f + 3.5e-5f * freq - 8e-10f * freq2 - 4e-14f * freq3; 
    K = 4.0f * res * qa;       
    float fa = 1.0036f - 2e-5f *freq + 1e-9f*freq2   - 1.75e-14f * freq3;  
    Gx2Vt = float((1.0f - exp(-fc * fa)) * VTx2);                                      
}

void AudioFilterLadder2::frequency(float c)
{
    if (c < 1) c = 1;
    if (c > MAX_FREQUENCY) c = MAX_FREQUENCY;  
    targetFbase = c; 
    //Fbase = c;   
}

void AudioFilterLadder2::resonance(float res)
{  
    res *= MAX_RESONANCE;                                 // normalize so host specifies 0-1 rather than 0-1.015
    if (res > MAX_RESONANCE) res = MAX_RESONANCE ;
    if (res < 0.0f) res = 0.0f;
    Qbase = res;
}

void AudioFilterLadder2::octaveControl(float octaves)
{
  if (octaves > 10.0f) {             // increased range compared to previous filters to span full audio range with 0-1 envelope
    octaves = 10.0f;
  } else if (octaves < 0.0f) {
    octaves = 0.0f;
  }
  octaveScale = octaves / 32768.0f;
}

static inline float fast_tanh(float x)
{
  if(x>3) return 1;
  if(x<-3) return -1;
  float x2 = x * x;
  return x * (27.0f + x2) / (27.0f + 9.0f * x2);
}

bool AudioFilterLadder2::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;
}

void AudioFilterLadder2::initpoly()  
{
    if (arm_fir_interpolate_init_f32(&interpolation2x, INTERPOLATION2x, finterpolation_taps2x, finterpolation_coeffs2x, finterpolation_state2x, AUDIO_BLOCK_SAMPLES))  
    {
      Serial.println("Init of interpolation 2x failed");
      while (1);
    }
    if (arm_fir_decimate_init_f32(&decimation2x, finterpolation_taps2x, INTERPOLATION2x, finterpolation_coeffs2x, fdecimation_state2x, fI_NUM_SAMPLES2x))
    {
      Serial.println("Init of decimation 2x failed");
      while (1);
    }     
} 

void AudioFilterLadder2::update(void)
{
  audio_block_t *blocka, *blockb, *blockc;
  float ftot, qtot;
  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;
        ftot = Fbase;
      }
      if (!blockc) {
        QmodActive = false;
        qtot = Qbase;
      }
      if(QmodActive == false && FCmodActive == false) 
        compute_coeffs2x(Fbase,Qbase) ;    

      float blockIn[AUDIO_BLOCK_SAMPLES], blockOut[AUDIO_BLOCK_SAMPLES];   
      float blockOS[fI_NUM_SAMPLES2x], blockOutOS[fI_NUM_SAMPLES2x];                
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++)
        blockIn[i] = blocka->data[i]  * overdrive *  float(INTERPOLATION2x /32768.0f);
      arm_fir_interpolate_f32(&interpolation2x, blockIn, blockOS, AUDIO_BLOCK_SAMPLES);             
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {                   
        Fbase = FcPorta * Fbase + (1-FcPorta) * targetFbase;             
        if (FCmodActive) {
          float FCmod = blockb->data[i] * octaveScale;
          ftot = Fbase * fast_exp2f(FCmod);                    
        }
        if (QmodActive) {
          float Qmod = blockc->data[i] * (1.0f/32768.0f);
          qtot = Qbase + Qmod;
        }            
        if(QmodActive || FCmodActive) compute_coeffs2x(ftot,qtot) ;                          
        for(int os = 0; os < INTERPOLATION2x; os++)  
        {      
          float input = blockOS[i * INTERPOLATION2x + os];                            
          filter_y1 = filter_y1 + Gx2Vt * (fast_tanh((input - K * filter_out) * inv2VT) - save_tan1);
          save_tan1=  fast_tanh(filter_y1 * inv2VT);                  
          filter_y2 = filter_y2 + Gx2Vt * (save_tan1 - save_tan2);
          save_tan2 = fast_tanh(filter_y2 * inv2VT);                  
          filter_y3 = filter_y3 + Gx2Vt * (save_tan2 - save_tan3);  
          save_tan3 = fast_tanh(filter_y3 * inv2VT);                   
          filter_y4 = filter_y4 + Gx2Vt * (save_tan3 - fast_tanh(filter_y4 * inv2VT));                             
          filter_out = (filter_y4 + filter_y5) * 0.5f;
          filter_y5 = filter_y4;   
          blockOutOS[i*INTERPOLATION2x + os] = filter_out;
        }      
      }
      arm_fir_decimate_f32(&decimation2x, blockOutOS, blockOut, fI_NUM_SAMPLES2x);                  
      float absblk;
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {  
        absblk = 1.2f * blockOut[i];
        if(absblk<0) absblk=-absblk;        
        if(absblk > 1)
           if(absblk > peak)
              peak = absblk;
      }     
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {         
        if(peak>1){ 
          pgain = 0.995f * pgain + 0.005f/peak;    // fast attack    
          peak *=0.99995f;
        }   
        blocka->data[i] = blockOut[i]  * pgain * float(0.9f * float(32767.0));  
      } 
      transmit(blocka);
      release(blocka);
      if (blockb) release(blockb);
      if (blockc) release(blockc);
}

And filter_ladder2.h

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.
 */

//-----------------------------------------------------------
// AudioFilterLadder2 for Teensy Audio Library
// Based on Huovilainen's full model, presented at DAFX 2004
// Richard van Hoesel, March. 4, 2021
// v.1.0 Uses 2x oversampling. Host resonance range is 0-1.0. 
// Self oscillation at lower levels is available by specifying 
// values just below 1.0   
//
// please retain this header if you use this code.
//-----------------------------------------------------------

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

#ifndef filter_ladder2_h_
#define filter_ladder2_h_

#include "Arduino.h"
#include "AudioStream.h"
#include "arm_math.h"
#define MOOG_PI ((float)3.14159265358979323846264338327950288)  

class AudioFilterLadder2: public AudioStream
{
public:
  AudioFilterLadder2() : AudioStream(3, inputQueueArray) { initpoly(); };
  void frequency(float FC);
  void resonance(float reson);
  void octaveControl(float octaves);
  void inputDrive(float drive);
  void portamento(float pcf);   
  void compute_coeffs2x(float freq, float res);  
  void initpoly(); 
  virtual void update(void);
  
private:
  static const int INTERPOLATION2x = 2;
  static const int finterpolation_taps2x = 32;
  static const float FC_SCALER2x = 2.0f * MOOG_PI / (INTERPOLATION2x * AUDIO_SAMPLE_RATE_EXACT);
  float finterpolation_state2x[(AUDIO_BLOCK_SAMPLES-1) + finterpolation_taps2x / INTERPOLATION2x];
  arm_fir_interpolate_instance_f32 interpolation2x;
  float fdecimation_state2x[(AUDIO_BLOCK_SAMPLES*INTERPOLATION2x-1) + finterpolation_taps2x];
  arm_fir_decimate_instance_f32 decimation2x;   
  static float finterpolation_coeffs2x[finterpolation_taps2x];
  
  float Gx2Vt, K, filter_res, filter_y1, filter_y2, filter_y3, filter_y4, filter_y5, filter_out; 
  float VTx2 = 10; 
  float inv2VT = 1.0/VTx2;
  float Fbase=1000, Qbase=0.5;
  float overdrive=1;
  float ftot=Fbase, qtot=Qbase;
  bool  resonating();
  float octaveScale = 1.0f/32768.0f;
  float z0[4] = {0.0, 0.0, 0.0, 0.0};
  float z1[4] = {0.0, 0.0, 0.0, 0.0};
  float pbg=0;
  float save_tan1=0;
  float save_tan2=0;
  float save_tan4=0;
  float save_tan3=0;
  float peak=1.0;
  float pgain=1.0;
  float targetFbase = 1000;
  float FcPorta = 0;     
  audio_block_t *inputQueueArray[3];
};

#endif
Feel free to send me your impressions.
Cheers,
Richard.
 

Attachments

  • AudioFilterLadder2_self_osc_demo.zip
    74.5 KB · Views: 99
Last edited:
Here is v1.01 (already....). Rather than dialing back resonance parameter to try to reduce self oscillation loudness, which has the disadvantage of often cutting out the self oscillation unintentionally, I've added a parameter "resonanceLoudness(float)" that specifically controls the loudness but not the 'sustained oscillation' while leaving the resonance parameter unchanged. It takes values between 1.0 (soft) and 10.0 (loud).

filter_ladder2.cpp:
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.
 */

//-----------------------------------------------------------
// AudioFilterLadder2 for Teensy Audio Library
// Based on Huovilainen's full model, presented at DAFX 2004
// Richard van Hoesel, March. 4, 2021
// v.1.0 Uses 2x oversampling. Host resonance range is 0-1.0. 
// Self oscillation at lower levels is available by specifying 
// values just below 1.0   
//
// please retain this header if you use this code.
//-----------------------------------------------------------

#include <Arduino.h>
#include "filter_ladder2.h"
#include "arm_math.h"
#include <math.h>
#include <stdint.h>

#define MAX_RESONANCE ((float)1.015)
#define MAX_FREQUENCY ((float)(AUDIO_SAMPLE_RATE_EXACT * 0.425f))
#define fI_NUM_SAMPLES2x  AUDIO_BLOCK_SAMPLES * INTERPOLATION2x

//=== 2 xOS, 32taps (88.2kHz)
static float AudioFilterLadder2 ::finterpolation_coeffs2x[AudioFilterLadder2::finterpolation_taps2x] = {    
-129.6800658538650740E-6, 0.001562843504973615, 0.004914813078250063, 0.007186231102125209, 0.002728024844356039,-0.007973931496234599,-0.013255882861166815,-841.1646993462468340E-6,
 0.022309967341058432, 0.027240071522569135,-0.007073857882406733,-0.056523975383091875,-0.055929860712812876, 0.042278471323346570, 0.207056768528768836, 0.335633514003638445,
 0.335633514003638445, 0.207056768528768836, 0.042278471323346570,-0.055929860712812876,-0.056523975383091875,-0.007073857882406733, 0.027240071522569135, 0.022309967341058432,
-841.1646993462468340E-6,-0.013255882861166815,-0.007973931496234599, 0.002728024844356039, 0.007186231102125209, 0.004914813078250063, 0.001562843504973615,-129.6800658538650740E-6
};

void AudioFilterLadder2::resonanceLoudness(float v)  // input 1-10
{
  if(v<1) v=1;
  if(v>10) v=10;
  VTx2 =  v;   
  inv2VT = 1.0/VTx2;
}

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;
}

void AudioFilterLadder2::inputDrive(float drive) 
{
  if(drive<1)
      overdrive = 1.0f;
  else 
  if(drive>2)
     overdrive =  2.0f;  
  else
    overdrive = drive;   
}

void AudioFilterLadder2::portamento(float pc)
{
  FcPorta = pc; 
}



void AudioFilterLadder2::compute_coeffs2x(float freq, float res) // tuned for 2x os
{
    float fc, freq2, freq3;          
    if (res > MAX_RESONANCE) res = MAX_RESONANCE ;
    if (res < 0.0f) res = 0.0f;    if (freq < 5) freq = 5;
    if (freq > MAX_FREQUENCY) freq = MAX_FREQUENCY;       
    fc = freq *  FC_SCALER2x;                             
    freq2 = freq*freq;
    freq3 = freq2 * freq;      
    float qa = 0.998f + 3.5e-5f * freq - 8e-10f * freq2 - 4e-14f * freq3; 
    K = 4.0f * res * qa;       
    float fa = 1.0036f - 2e-5f *freq + 1e-9f*freq2   - 1.75e-14f * freq3;  
    Gx2Vt = float((1.0f - exp(-fc * fa)) * VTx2);                                      
}

void AudioFilterLadder2::frequency(float c)
{
    if (c < 1) c = 1;
    if (c > MAX_FREQUENCY) c = MAX_FREQUENCY;  
    targetFbase = c; 
    //Fbase = c;   
}

void AudioFilterLadder2::resonance(float res)
{  
    res *= MAX_RESONANCE;                                 // normalize so host specifies 0-1 rather than 0-1.015
    if (res > MAX_RESONANCE) res = MAX_RESONANCE ;
    if (res < 0.0f) res = 0.0f;
    Qbase = res;
}

void AudioFilterLadder2::octaveControl(float octaves)
{
  if (octaves > 10.0f) {             // increased range compared to previous filters to span full audio range with 0-1 envelope
    octaves = 10.0f;
  } else if (octaves < 0.0f) {
    octaves = 0.0f;
  }
  octaveScale = octaves / 32768.0f;
}

static inline float fast_tanh(float x)
{
  if(x>3) return 1;
  if(x<-3) return -1;
  float x2 = x * x;
  return x * (27.0f + x2) / (27.0f + 9.0f * x2);
}

bool AudioFilterLadder2::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;
}

void AudioFilterLadder2::initpoly()  
{
    if (arm_fir_interpolate_init_f32(&interpolation2x, INTERPOLATION2x, finterpolation_taps2x, finterpolation_coeffs2x, finterpolation_state2x, AUDIO_BLOCK_SAMPLES))  
    {
      Serial.println("Init of interpolation 2x failed");
      while (1);
    }
    if (arm_fir_decimate_init_f32(&decimation2x, finterpolation_taps2x, INTERPOLATION2x, finterpolation_coeffs2x, fdecimation_state2x, fI_NUM_SAMPLES2x))
    {
      Serial.println("Init of decimation 2x failed");
      while (1);
    }     
} 

void AudioFilterLadder2::update(void)
{
  audio_block_t *blocka, *blockb, *blockc;
  float ftot, qtot;
  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;
        ftot = Fbase;
      }
      if (!blockc) {
        QmodActive = false;
        qtot = Qbase;
      }

      float blockIn[AUDIO_BLOCK_SAMPLES], blockOut[AUDIO_BLOCK_SAMPLES];   
      float blockOS[fI_NUM_SAMPLES2x], blockOutOS[fI_NUM_SAMPLES2x];                
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++)
        blockIn[i] = blocka->data[i]  * overdrive *  float(INTERPOLATION2x /32768.0f);
      arm_fir_interpolate_f32(&interpolation2x, blockIn, blockOS, AUDIO_BLOCK_SAMPLES);             
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {                   
        Fbase = FcPorta * Fbase + (1-FcPorta) * targetFbase;             
        if (FCmodActive) {
          float FCmod = blockb->data[i] * octaveScale;
          ftot = Fbase * fast_exp2f(FCmod);                    
        }
        if (QmodActive) {
          float Qmod = blockc->data[i] * (1.0f/32768.0f);
          qtot = Qbase + Qmod;
        }            
        compute_coeffs2x(ftot,qtot) ;                          
        for(int os = 0; os < INTERPOLATION2x; os++)  
        {      
          float input = blockOS[i * INTERPOLATION2x + os];                            
          filter_y1 = filter_y1 + Gx2Vt * (fast_tanh((input - K * filter_out) * inv2VT) - save_tan1);
          save_tan1=  fast_tanh(filter_y1 * inv2VT);                  
          filter_y2 = filter_y2 + Gx2Vt * (save_tan1 - save_tan2);
          save_tan2 = fast_tanh(filter_y2 * inv2VT);                  
          filter_y3 = filter_y3 + Gx2Vt * (save_tan2 - save_tan3);  
          save_tan3 = fast_tanh(filter_y3 * inv2VT);                   
          filter_y4 = filter_y4 + Gx2Vt * (save_tan3 - fast_tanh(filter_y4 * inv2VT));                             
          filter_out = (filter_y4 + filter_y5) * 0.5f;
          filter_y5 = filter_y4;   
          blockOutOS[i*INTERPOLATION2x + os] = filter_out;
        }      
      }
      arm_fir_decimate_f32(&decimation2x, blockOutOS, blockOut, fI_NUM_SAMPLES2x);                  
      float absblk;
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {  
        absblk = 1.2f * blockOut[i];
        if(absblk<0) absblk=-absblk;        
        if(absblk > 1)
           if(absblk > peak)
              peak = absblk;
      }     
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {         
        if(peak>1){ 
          pgain = 0.995f * pgain + 0.005f/peak;    // fast attack    
          peak *=0.99995f;
        }   
        blocka->data[i] = blockOut[i]  * pgain * float(0.9f * float(32767.0));  
      } 
      transmit(blocka);
      release(blocka);
      if (blockb) release(blockb);
      if (blockc) release(blockc);
}

filter_ladder2.h
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.
 */

//-----------------------------------------------------------
// AudioFilterLadder2 for Teensy Audio Library
// Based on Huovilainen's full model, presented at DAFX 2004
// Richard van Hoesel, March. 4, 2021
// v.1.01 Adds a "resonanceLoudness" parameter (1.0 - 10.0)
//        to control resonance loudness 
// v.1.0 Uses 2x oversampling. Host resonance range is 0-1.0. 
//
// please retain this header if you use this code.
//-----------------------------------------------------------

// https://forum.pjrc.com/threads/60488?p=269609&viewfull=1#post269609

#ifndef filter_ladder2_h_
#define filter_ladder2_h_

#include "Arduino.h"
#include "AudioStream.h"
#include "arm_math.h"
#define MOOG_PI ((float)3.14159265358979323846264338327950288)  

class AudioFilterLadder2: public AudioStream
{
public:
  AudioFilterLadder2() : AudioStream(3, inputQueueArray) { initpoly(); };
  void frequency(float FC);
  void resonance(float reson);
  void octaveControl(float octaves);
  void inputDrive(float drive);
  void portamento(float pcf);   
  void compute_coeffs2x(float freq, float res);  
  void initpoly(); 
  void resonanceLoudness(float v);
  virtual void update(void);
  
private:
  static const int INTERPOLATION2x = 2;
  static const int finterpolation_taps2x = 32;
  static const float FC_SCALER2x = 2.0f * MOOG_PI / (INTERPOLATION2x * AUDIO_SAMPLE_RATE_EXACT);
  float finterpolation_state2x[(AUDIO_BLOCK_SAMPLES-1) + finterpolation_taps2x / INTERPOLATION2x];
  arm_fir_interpolate_instance_f32 interpolation2x;
  float fdecimation_state2x[(AUDIO_BLOCK_SAMPLES*INTERPOLATION2x-1) + finterpolation_taps2x];
  arm_fir_decimate_instance_f32 decimation2x;   
  static float finterpolation_coeffs2x[finterpolation_taps2x];
  
  float Gx2Vt, K, filter_res, filter_y1, filter_y2, filter_y3, filter_y4, filter_y5, filter_out; 
  float VTx2 = 10; 
  float inv2VT = 1.0/VTx2;
  float Fbase=1000, Qbase=0.5;
  float overdrive=1;
  float ftot=Fbase, qtot=Qbase;
  bool  resonating();
  float octaveScale = 1.0f/32768.0f;
  float z0[4] = {0.0, 0.0, 0.0, 0.0};
  float z1[4] = {0.0, 0.0, 0.0, 0.0};
  float pbg=0;
  float save_tan1=0;
  float save_tan2=0;
  float save_tan4=0;
  float save_tan3=0;
  float peak=1.0;
  float pgain=1.0;
  float targetFbase = 1000;
  float FcPorta = 0;     
  audio_block_t *inputQueueArray[3];
};

#endif
 
...and although the .inputDrive parameter is still in place, it's not really a useful and is limited to 1-2.

I have changed my mind about this parameter for this second filter. I'm finding now that allowing it to range between 1-3 seems to be a good way to 'fatten' the sound without altering the character of the filter too much otherwise (which was a problem for me in the first filter model). Note that you'll need to change the inputDrive function in the .cpp file to set it beyond 2.
 
Here is the revised .cpp file that allows inputDrive 1-3, and also makes a small change to a time constant.
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.
 */

//-----------------------------------------------------------
// AudioFilterLadder2 for Teensy Audio Library
// Based on Huovilainen's full model, presented at DAFX 2004
// Richard van Hoesel, March. 4, 2021
// v.1.01 Adds a "resonanceLoudness" parameter (1.0 - 10.0)
//        to control resonance loudness 
// v.1.0 Uses 2x oversampling. Host resonance range is 0-1.0. 
//
// please retain this header if you use this code.
//-----------------------------------------------------------

#include <Arduino.h>
#include "filter_ladder2.h"
#include "arm_math.h"
#include <math.h>
#include <stdint.h>

#define MAX_RESONANCE ((float)1.015)
#define MAX_FREQUENCY ((float)(AUDIO_SAMPLE_RATE_EXACT * 0.425f))
#define fI_NUM_SAMPLES2x  AUDIO_BLOCK_SAMPLES * INTERPOLATION2x

//=== 2 xOS, 32taps (88.2kHz)
static float AudioFilterLadder2 ::finterpolation_coeffs2x[AudioFilterLadder2::finterpolation_taps2x] = {    
-129.6800658538650740E-6, 0.001562843504973615, 0.004914813078250063, 0.007186231102125209, 0.002728024844356039,-0.007973931496234599,-0.013255882861166815,-841.1646993462468340E-6,
 0.022309967341058432, 0.027240071522569135,-0.007073857882406733,-0.056523975383091875,-0.055929860712812876, 0.042278471323346570, 0.207056768528768836, 0.335633514003638445,
 0.335633514003638445, 0.207056768528768836, 0.042278471323346570,-0.055929860712812876,-0.056523975383091875,-0.007073857882406733, 0.027240071522569135, 0.022309967341058432,
-841.1646993462468340E-6,-0.013255882861166815,-0.007973931496234599, 0.002728024844356039, 0.007186231102125209, 0.004914813078250063, 0.001562843504973615,-129.6800658538650740E-6
};

void AudioFilterLadder2::resonanceLoudness(float v)  // input 1-10
{
  if(v<1) v=1;
  if(v>10) v=10;
  VTx2 =  v;   
  inv2VT = 1.0/VTx2;
}

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;
}

void AudioFilterLadder2::inputDrive(float drive) 
{
  if(drive<1)
      overdrive = 1.0f;
  else 
  if(drive>3)
     overdrive =  3.0f;  
  else
    overdrive = drive;   
}

void AudioFilterLadder2::portamento(float pc)
{
  FcPorta = pc; 
}

void AudioFilterLadder2::compute_coeffs2x(float freq, float res) // tuned for 2x os
{
    float fc, freq2, freq3;          
    if (res > MAX_RESONANCE) res = MAX_RESONANCE ;
    if (res < 0.0f) res = 0.0f;    if (freq < 5) freq = 5;
    if (freq > MAX_FREQUENCY) freq = MAX_FREQUENCY;       
    fc = freq *  FC_SCALER2x;                             
    freq2 = freq*freq;
    freq3 = freq2 * freq;      
    float qa = 0.998f + 3.5e-5f * freq - 8e-10f * freq2 - 4e-14f * freq3; 
    K = 4.0f * res * qa;       
    float fa = 1.0036f - 2e-5f *freq + 1e-9f*freq2   - 1.75e-14f * freq3;  
    Gx2Vt = float((1.0f - exp(-fc * fa)) * VTx2);                                      
}

void AudioFilterLadder2::frequency(float c)
{
    if (c < 1) c = 1;
    if (c > MAX_FREQUENCY) c = MAX_FREQUENCY;  
    targetFbase = c; 
    //Fbase = c;   
}

void AudioFilterLadder2::resonance(float res)
{  
    res *= MAX_RESONANCE;                                 // normalize so host specifies 0-1 rather than 0-1.015
    if (res > MAX_RESONANCE) res = MAX_RESONANCE ;
    if (res < 0.0f) res = 0.0f;
    Qbase = res;
}

void AudioFilterLadder2::octaveControl(float octaves)
{
  if (octaves > 10.0f) {             // increased range compared to previous filters to span full audio range with 0-1 envelope
    octaves = 10.0f;
  } else if (octaves < 0.0f) {
    octaves = 0.0f;
  }
  octaveScale = octaves / 32768.0f;
}

static inline float fast_tanh(float x)
{
  if(x>3) return 1;
  if(x<-3) return -1;
  float x2 = x * x;
  return x * (27.0f + x2) / (27.0f + 9.0f * x2);
}

bool AudioFilterLadder2::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;
}

void AudioFilterLadder2::initpoly()  
{
    if (arm_fir_interpolate_init_f32(&interpolation2x, INTERPOLATION2x, finterpolation_taps2x, finterpolation_coeffs2x, finterpolation_state2x, AUDIO_BLOCK_SAMPLES))  
    {
      Serial.println("Init of interpolation 2x failed");
      while (1);
    }
    if (arm_fir_decimate_init_f32(&decimation2x, finterpolation_taps2x, INTERPOLATION2x, finterpolation_coeffs2x, fdecimation_state2x, fI_NUM_SAMPLES2x))
    {
      Serial.println("Init of decimation 2x failed");
      while (1);
    }     
} 

void AudioFilterLadder2::update(void)
{
  audio_block_t *blocka, *blockb, *blockc;
  float ftot, qtot;
  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;
        ftot = Fbase;
      }
      if (!blockc) {
        QmodActive = false;
        qtot = Qbase;
      }

      float blockIn[AUDIO_BLOCK_SAMPLES], blockOut[AUDIO_BLOCK_SAMPLES];   
      float blockOS[fI_NUM_SAMPLES2x], blockOutOS[fI_NUM_SAMPLES2x];                
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++)
        blockIn[i] = blocka->data[i]  * overdrive *  float(INTERPOLATION2x /32768.0f);
      arm_fir_interpolate_f32(&interpolation2x, blockIn, blockOS, AUDIO_BLOCK_SAMPLES);             
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {                   
        Fbase = FcPorta * Fbase + (1-FcPorta) * targetFbase;             
        if (FCmodActive) {
          float FCmod = blockb->data[i] * octaveScale;
          ftot = Fbase * fast_exp2f(FCmod);                    
        }
        if (QmodActive) {
          float Qmod = blockc->data[i] * (1.0f/32768.0f);
          qtot = Qbase + Qmod;
        }            
        compute_coeffs2x(ftot,qtot) ;                          
        for(int os = 0; os < INTERPOLATION2x; os++)  
        {      
          float input = blockOS[i * INTERPOLATION2x + os];                            
          filter_y1 = filter_y1 + Gx2Vt * (fast_tanh((input - K * filter_out) * inv2VT) - save_tan1);
          save_tan1=  fast_tanh(filter_y1 * inv2VT);                  
          filter_y2 = filter_y2 + Gx2Vt * (save_tan1 - save_tan2);
          save_tan2 = fast_tanh(filter_y2 * inv2VT);                  
          filter_y3 = filter_y3 + Gx2Vt * (save_tan2 - save_tan3);  
          save_tan3 = fast_tanh(filter_y3 * inv2VT);                   
          filter_y4 = filter_y4 + Gx2Vt * (save_tan3 - fast_tanh(filter_y4 * inv2VT));                             
          filter_out = (filter_y4 + filter_y5) * 0.5f;
          filter_y5 = filter_y4;   
          blockOutOS[i*INTERPOLATION2x + os] = filter_out;
        }      
      }
      arm_fir_decimate_f32(&decimation2x, blockOutOS, blockOut, fI_NUM_SAMPLES2x);                  
      float absblk;
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {  
        absblk = 1.25f * blockOut[i];
        if(absblk<0) absblk=-absblk;        
        if(absblk > 1)
           if(absblk > peak)
              peak = absblk;
      }     
      for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {         
        if(peak>1){ 
          pgain = 0.99f * pgain + 0.01f/peak;    // fast attack    
          peak *=0.99995f;
        }   
        blocka->data[i] = blockOut[i]  * pgain * float(0.85f * float(32767.0));  
      } 
      transmit(blocka);
      release(blocka);
      if (blockb) release(blockb);
      if (blockc) release(blockc);
}

And also a revised .h that initializes the resonanceLoudness to a mid-level loudness of 5 rather 10:

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.
 */

//-----------------------------------------------------------
// AudioFilterLadder2 for Teensy Audio Library
// Based on Huovilainen's full model, presented at DAFX 2004
// Richard van Hoesel, March. 4, 2021
// v.1.01 Adds a "resonanceLoudness" parameter (1.0 - 10.0)
//        to control resonance loudness 
// v.1.0 Uses 2x oversampling. Host resonance range is 0-1.0. 
//
// please retain this header if you use this code.
//-----------------------------------------------------------

// https://forum.pjrc.com/threads/60488?p=269609&viewfull=1#post269609

#ifndef filter_ladder2_h_
#define filter_ladder2_h_

#include "Arduino.h"
#include "AudioStream.h"
#include "arm_math.h"
#define MOOG_PI ((float)3.14159265358979323846264338327950288)  

class AudioFilterLadder2: public AudioStream
{
public:
  AudioFilterLadder2() : AudioStream(3, inputQueueArray) { initpoly(); };
  void frequency(float FC);
  void resonance(float reson);
  void octaveControl(float octaves);
  void inputDrive(float drive);
  void portamento(float pcf);   
  void compute_coeffs2x(float freq, float res);  
  void initpoly(); 
  void resonanceLoudness(float v);
  virtual void update(void);
  
private:
  static const int INTERPOLATION2x = 2;
  static const int finterpolation_taps2x = 32;
  static const float FC_SCALER2x = 2.0f * MOOG_PI / (INTERPOLATION2x * AUDIO_SAMPLE_RATE_EXACT);
  float finterpolation_state2x[(AUDIO_BLOCK_SAMPLES-1) + finterpolation_taps2x / INTERPOLATION2x];
  arm_fir_interpolate_instance_f32 interpolation2x;
  float fdecimation_state2x[(AUDIO_BLOCK_SAMPLES*INTERPOLATION2x-1) + finterpolation_taps2x];
  arm_fir_decimate_instance_f32 decimation2x;   
  static float finterpolation_coeffs2x[finterpolation_taps2x];
  
  float Gx2Vt, K, filter_res, filter_y1, filter_y2, filter_y3, filter_y4, filter_y5, filter_out; 
  float VTx2 = 5; 
  float inv2VT = 1.0/VTx2;
  float Fbase=1000, Qbase=0.5;
  float overdrive=1;
  float ftot=Fbase, qtot=Qbase;
  bool  resonating();
  float octaveScale = 1.0f/32768.0f;
  float z0[4] = {0.0, 0.0, 0.0, 0.0};
  float z1[4] = {0.0, 0.0, 0.0, 0.0};
  float pbg=0;
  float save_tan1=0;
  float save_tan2=0;
  float save_tan4=0;
  float save_tan3=0;
  float peak=1.0;
  float pgain=1.0;
  float targetFbase = 1000;
  float FcPorta = 0;     
  audio_block_t *inputQueueArray[3];
};

#endif
 
Last edited:
Hello all.
First, many thanks for all the effort, the filter sounds great! I am using the "simplified model", that is currently in the official Audio Library.
But it really is ressource-hungry. On my Teensy 4 running at 600MHz, the filter needs a bit over 7 % CPU time. With linear interpolation, it still needs 6 % CPU time.
For my 8 voice polyphony synth, who still has to handle lots of stuff in the background, this is too much and crashes the synth, evey time I use all voices.

Now would be my question:
Is there a way, I can reduce the CPU load further? Do I need to stick with version 1.02 (I need cutoff CV input) or is there any other way to work with the filter, that is currently in the library (e.g. disable the oversampling?) If I can manage to get under 4% CPU load, that would be great.

Greetings, Florian
 
Hello all.
First, many thanks for all the effort, the filter sounds great! I am using the "simplified model", that is currently in the official Audio Library.
But it really is ressource-hungry. On my Teensy 4 running at 600MHz, the filter needs a bit over 7 % CPU time. With linear interpolation, it still needs 6 % CPU time.
For my 8 voice polyphony synth, who still has to handle lots of stuff in the background, this is too much and crashes the synth, evey time I use all voices.

Now would be my question:
Is there a way, I can reduce the CPU load further? Do I need to stick with version 1.02 (I need cutoff CV input) or is there any other way to work with the filter, that is currently in the library (e.g. disable the oversampling?) If I can manage to get under 4% CPU load, that would be great.

Greetings, Florian

Hi Florian, good to hear you're enjoying the filter. On my system, I get 7% for POLY_FIR, and 5% for LINEAR. To reduce it further I suggest going to 2x oversampling. However, that requires a different set of FIR coefficients. The ones from the newer filter that already runs at 2x should be fine, but they are 32pt rather than 36 so the changes to the .h file include changing both INTERPOLATION from 4 to 2, and interpolation_taps form 36 to 32.

The coefficients needed in the cpp file are:

Code:
static float AudioFilterLadder ::interpolation_coeffs[AudioFilterLadder::interpolation_taps] = {    
-129.6800658538650740E-6, 0.001562843504973615, 0.004914813078250063, 0.007186231102125209, 0.002728024844356039,-0.007973931496234599,-0.013255882861166815,-841.1646993462468340E-6,
 0.022309967341058432, 0.027240071522569135,-0.007073857882406733,-0.056523975383091875,-0.055929860712812876, 0.042278471323346570, 0.207056768528768836, 0.335633514003638445,
 0.335633514003638445, 0.207056768528768836, 0.042278471323346570,-0.055929860712812876,-0.056523975383091875,-0.007073857882406733, 0.027240071522569135, 0.022309967341058432,
-841.1646993462468340E-6,-0.013255882861166815,-0.007973931496234599, 0.002728024844356039, 0.007186231102125209, 0.004914813078250063, 0.001562843504973615,-129.6800658538650740E-6
};

However when I tried this I discovered there is a bug in the code of the poly update routine that only shows up if you are not using 4x oversampling. Specifically, in the inner loop the indexing into the audio block at input/output currently uses "i*4" but it should be "i*INTERPOLATION"

Below is the corrected code for the inner loop in poly mode.
Code:
for(int os=0; os < INTERPOLATION; os++) {
  float input = blockOS[i*INTERPOLATION + os];
  float u = input - (z1[3] - pbg * input) * Ktot * Qadjust;
  u = fast_tanh(u);
  float stage1 = LPF(u, 0);
  float stage2 = LPF(stage1, 1);
  float stage3 = LPF(stage2, 2);
  float stage4 = LPF(stage3, 3);
  blockOutOS[i*INTERPOLATION + os] = stage4;
}

Hopefully Paul can correct this in the Audio lib version. If there is general interest in having a 2x oversampling option to get lower CPU use, I will defer to Paul to figure out the best way to include the options I've described above.

As an aside, from the brief listening test I did it seems to sound surprisingly (to me) similar to the 4x OS version. So if Paul is happy to incorporate the 2x option into the library version of the filter, perhaps 2x should be the default. The max CPU usage measures I see with the various configurations are:
4x poly 7%
4x linear 5%
2x poly 5%
2x linear 4%
(is there a way to see fractional % ?)

When %CPU is less critical, I highly recommend trying the full model filter instead which has far lower 'tanh aliasing' at high resonance.

Cheers,
Richard
 
Last edited:
Back
Top