Porting moog ladder filters to audio objects?

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.
Sure, here is a bit copy&paste:

The program calculates the filter on the fly, so if you want to that, you'll need this file: output_fm_math.cpp
You can do that one time and print the values and fill a const array with that. I do that dynamically, because I played a lot with the parameters, and it was faster this way.
So, for the coeffs, use can use a filter design tool, or this:

INTERPOLATION = 4
I_SAMPLERATE = SAMPLERATE * INTERPOLATION
NUM_SAMPLES = 128
I_NUM_SAMPLES = NUM_SAMPLES * INTERPOLATION


Code:
const unsigned interpolation_taps = 64;
const float n_att = 90.0; // desired stopband attenuation
const float interpolation_cutoff = 15000.0f;  // AUDIO_SAMPLE_RATE_EXACT / 2.0f; //
float interpolation_coeffs[interpolation_taps]; // <- coeffs
[...]
extern void calc_FIR_coeffs (float * coeffs_I, int numCoeffs, float fc, float Astop, int type, float dfc, float Fsamprate);
[...]
//calculate the coeffs:
  //calc_FIR_coeffs (float * coeffs_I, int numCoeffs, float fc, float Astop, int type, float dfc, float Fsamprate)
  calc_FIR_coeffs(interpolation_coeffs, interpolation_taps, interpolation_cutoff, n_att, 0, 0.0, I_SAMPLERATE);
This code in output_fm_math is not by me, and i can't says much about it.
Cutoff is 15kHZ here. You may want to use half of the samplerate (22kHZ) or way less.
Using less taps will be faster.

Variables for interpolation:
Code:
float interpolation_L_state[(NUM_SAMPLES-1) + interpolation_taps / INTERPOLATION];
arm_fir_interpolate_instance_f32 interpolationL;
Initialize the interpolation:
Code:
  // we initiate the number of input samples, NOT the number of interpolated samples
  if (arm_fir_interpolate_init_f32(&interpolationL, INTERPOLATION, interpolation_taps, interpolation_coeffs, interpolation_L_state, NUM_SAMPLES)
  {
    Serial.println("Init of interpolation failed");
    while (1);
  }

do the interpolation:
I_NUM_SAMPLES = interpolation factor * 128
bL float array[128] of source,
iL interpolated output
Code:
 //interpolation:

  float iL[I_NUM_SAMPLES];

  arm_fir_interpolate_f32(&interpolationL, bL, iL, NUM_SAMPLES);
  //Scaling after interpolation
  //This is done later, so we can save some multiplications
  //arm_scale_f32(iL, (float)INTERPOLATION, iL, I_NUM_SAMPLES);

Done :)

Decimation is similar.
 

Thanks Paul,
I took that code and added simple linear interpolation of the input data in the version I show below. It extends the freq range before aliasing becomes audible with a bandlimited sawtooth wave running into the filter at high resonance. I think it's very usable as is now, but I'll also experiment a bit with alternative interpolation methods along the lines of the discussion above.
I also removed the "#ifdef lfq" related sections out of the code as I'm pretty sure that option won't be revisited in further developments.
 
/* 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. 18 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 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=269755&viewfull=1#post269755
// https://forum.pjrc.com/threads/60488?p=269609&viewfull=1#post269609

#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))

float AudioFilterLadder::LPF(float s, int i)
{
float ft = s * (1.0f/1.3f) + (0.3f/1.3f) * z0 - z1;
ft = ft * alpha + z1;
z1 = ft;
z0 = 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.5f) pbg = 0.5f;
if (pbg < 0.0f) pbg = 0.0f;
input_drive(host_overdrive);
}

void AudioFilterLadder::input_drive(float odrv)
{
host_overdrive = odrv;
if (host_overdrive > 1.0f) {
if (host_overdrive > 4.0f) host_overdrive = 4.0f;
// max is 4 when pbg = 0, and 2.5 when pbg is 0.5
overdrive = 1.0f + (host_overdrive - 1.0f) * (1.0f - pbg);
} else {
overdrive = host_overdrive;
if (overdrive < 0.0f) 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 / ((float)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) > 0.0001f) return true;
if (fabsf(z1) > 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, Kmax;
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 = 0;
}
}
if (!blockb) {
FCmodActive = false;
}
if (!blockc) {
QmodActive = false;
Ktot = K;
}
for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
if (FCmodActive) {
float FCmod = blockb->data * octaveScale;
float ftot = Fbase * fast_exp2f(FCmod);
if (ftot > MAX_FREQUENCY) ftot = MAX_FREQUENCY;
compute_coeffs(ftot);
}
if (QmodActive) {
float Qmod = blockc->data * (1.0f/32768.0f);
Ktot = K + 4.0f * Qmod;
}
Kmax = MAX_RESONANCE * 4.0f;
if (Ktot > Kmax) {
Ktot = Kmax;
} else if (Ktot < 0.0f) {
Ktot = 0.0f;
}
float total = 0.0f;
float interp = 0;
float input = blocka->data * (1.0f/32768.0f) * overdrive;
for(int os = 0; os < osTimes; 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 * (1.0f / (float)osTimes);
interp += float(1.0f/osTimes);
}
oldinput=input;
blocka->data = total * 32767.0f;
}
transmit(blocka);
release(blocka);
if (blockb) release(blockb);
if (blockc) release(blockc);
}
 
/* 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, v. 1.04, Feb. 18 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 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_

#include "Arduino.h"
#include "AudioStream.h"

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);
virtual void update(void);
private:
float LPF(float s, int i);
void compute_coeffs(float fc);
bool resonating();
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 = 2.8;
float Fbase = 800;
float Qadjust = 1.0f;
float octaveScale = 1.0f/32768.0f;
float pbg = 0.0f;
float overdrive = 1.0f;
float host_overdrive = 1.0f;
float oldinput = 0;
audio_block_t *inputQueueArray[3];
};

#endif
 
OK, Just downloaded and installed it in my test Teensy and I'm really starting to hear it. The self oscillation is there at a level that to me is actually preferable to the real beast. That's just a personal preference though. The rest of the range is starting to get the warmth that I was hoping for but at this stage the test rig is still just based on the Teensy and three pots. If I can get myself moving on it over the weekend I'm looking at rigging it to some op-amps so that I can get it leveled up into some semi-modular situations and see what it does, possibly with some analog sources containing some harmonics for it to get it's teeth into. At the moment it is a little too pristine and I really suspect that it is due to the very simple waveforms we're trying it with at the moment. I think I also have to get a keyboard patched into it for a better sense of the effect but as it stands it is definitely moving in the right direction if it hasn't gotten there already. From what I have seen most clones of the ladders avoid driving to the same degree as the Minimoog and I don't think that that is necessarily a bad thing. It was a haphazard result of a mathematical mistake after all and wasn't fine tuned for production. I've always thought my amateurish self that it didn't need all of that reach into the self oscillating range. As yours stands it has enough of this to experiment with in potential designs. I have an idea or two already. Kudos folks. You are doing things I didn't think you would get to. Even if it doesn't emulate the original perfectly suspect that it will still offer a very nice option when finished (if we can manage to avoid driving the resource tax too high).
 
What is CPU/Memory usage like on this? Do we see a polysynth using this? Or more like one teensy in a eurorack module?

Yes, a polysynth is easily possible. I have a preliminary 8-note polyphonic synth running with two bandwidth-limited audio oscillators, two envelope generators (vca/vcf), and one of these filters per voice, and the total comes in at about 40% CPU usage on a teensy 4.0.
 
@RVH,
do you need help with the oversampling?

@Frank, thanks very much for your offer. I plan to give it a go this weekend, and will contact you if I run into problems. One question - when I did a test run of the coefficient generating code for 16 coefficients and cutoff of 17kHz, all the coefficients were positive, whereas I was expecting a sinc function of some sort. Also, despite the length being even, the samples were offset rather than symmetric. Does all this sound right to you?
 
Sure, here is a bit copy&paste:
...

Decimation is similar.

OK, so I think I have the upsampling side up and running. I decided for now to just try an external design software (using the Iowa Hills link you sent) and it produced coefficients that were more like I was expecting (windowed sinc function). I decided I would 32 points initially to see how that went.

The structure has been initialised with:

#define INTERPOLATION 4
#define I_SAMPLERATE AUDIO_SAMPLE_RATE_EXACT * INTERPOLATION
#define NUM_SAMPLES 128
#define I_NUM_SAMPLES NUM_SAMPLES * INTERPOLATION

const unsigned interpolation_taps = 32;
float interpolation_L_state[(NUM_SAMPLES-1) + interpolation_taps / INTERPOLATION];
arm_fir_interpolate_instance_f32 interpolationL;
arm_fir_interpolate_init_f32(&interpolationL, INTERPOLATION, interpolation_taps, interpolation_coeffs, interpolation_L_state, NUM_SAMPLES))


Here is the code I have in place in the audio object update routine:
...
/*---------------- upsample (block rate) --------------------*/
float iL[I_NUM_SAMPLES], blockIn[AUDIO_BLOCK_SAMPLES];
for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++)
blockIn = blocka->data * overdrive * (1.0f/32768.0f);
arm_fir_interpolate_f32(&interpolationL, blockIn, iL, NUM_SAMPLES); // iL is interpolated buffer

for (int i=0; i < AUDIO_BLOCK_SAMPLES; i++) {
//---- this part runs at 44kHz (modulators update rate)

if (FCmodActive) {
float FCmod = blockb->data * octaveScale;
float ftot = Fbase * fast_exp2f(FCmod);
if (ftot > MAX_FREQUENCY) ftot = MAX_FREQUENCY;
compute_coeffs(ftot);
}
if (QmodActive) {
float Qmod = blockc->data * (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;

//---- this part runs at 4x 44kHz (audio filter)

for(int os = 0; os < osTimes; os++) // calc OS number of samples from upsampled buffer now
{
float input = iL[i*4 + os] * overdrive * (1.0f/32768.0f);
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);
total += stage4 * float(1.0/osTimes);
}
blocka->data = total * 32767.0f ;
}

The CPU usage with this 32pt input interpolation increased to about 6%, whereas with linear interpolation it was more like 4%.

At this stage I'm still using a linear average of the 4 interpolated filtered samples to get back to 44kHz. If instead I write the filter output back to the iL buffer as each sample is processed, is the decimation call I need as follows: arm_fir_interpolate_f32(&interpolationL, block->data, iL, I_NUM_SAMPLES) ?

Thanks.
 
is the decimation call I need as follows: arm_fir_interpolate_f32(&interpolationL, block->data, iL, I_NUM_SAMPLES) ?
Even allowing for the fact that I meant to use the decimate rather than interpolate functon, this turns out to be incorrect anyway. It looks like I need to create a separate structure for the decimation. I am assuming however that I can use the same filter coefficients. Please let me know if that's not the case.
 
I now have both a polyphase FIR interpolator and decimator up and running. CPU use for the filter on teensy 4 is still not bad at 7% for a 32-tap FIR. I'll do a bunch of listening tests for the various approaches, including different FIR designs and filter lengths, and will report back with my findings in due course.
 
I think 16 taps will be a good compromise between performance and CPU demands, and it sounds much better than linear interpolation at high frequencies.

However, I can't seem to run more than one filter at a time without severe audio problems, even though CPU use is only about 5% per filter, and increasing AudioMemory doesn't help either.
My guess is that the arm interpolation structures are not properly set up when more than one filter is created. Here is the related code at the start of my filter cpp filter file:

#define INTERPOLATION 4
#define I_SAMPLERATE AUDIO_SAMPLE_RATE_EXACT * INTERPOLATION
#define NUM_SAMPLES 128
#define I_NUM_SAMPLES NUM_SAMPLES * INTERPOLATION

const unsigned interpolation_taps = 16;
float interpolation_coeffs[interpolation_taps] = {
-0.003663756664110023,-0.013269359520045369,-0.021090072695910182,-0.011622057331749943, 0.029559652443450313, 0.101695838499519950, 0.182097194316226174, 0.235527892430298036,
0.235527892430298036, 0.182097194316226174, 0.101695838499519950, 0.029559652443450313,-0.011622057331749943,-0.021090072695910182,0.013269359520045369,-0.003663756664110023
};

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

void myAudioFilterLadder::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);
}
else
Serial.println("Init of interpolation OK");

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);
}
else
Serial.println("Init of decimation OK");
}

void myAudioFilterLadder::interpMethod(int imethod)
{
if(imethod==FIR_POLY) {
initpoly();
polyOn=true;
}
}

From my main code, I have assumed this sets up separate structures for each filter #N by calling flterN.interpMethod(FIR_POLY). The serial monitor messages seem to indicate the initializations succeed. But even for just two filters, if both are included in my Audioconnections, I get severe distortion in the output. When I disconnect one, the audio is fine.
For completenes, hereis the audio update routine's filter section:

/*----------------------- 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 = blocka->data * overdrive * (1.0f/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 * octaveScale;
float ftot = Fbase * fast_exp2f(FCmod);
if (ftot > MAX_FREQUENCY) ftot = MAX_FREQUENCY;
compute_coeffs(ftot);
}
if (QmodActive) {
float Qmod = blockc->data * (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;
float total = 0;
for(int os = 0; os < INTERPOLATION; os++) // calc OS number of samples from upsampled buffer now
{
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 = blockOut * float(32768.0) * INTERPOLATIONscaler;

//--------------------------------------------------------------------------------------------------

Please advise if you can see what I am doing wrong. Thanks.

p.s. if I use linear interpolation and don't set up the polyphase structures, I can run multiple filters just fine.
 
From what I can tell from the docs for the interpolation and decimation functions you are doing things right.

Have you tried without any of your intermediate processing, just interpolation straight back to decimation?

Have you tried adding some guard words to the state arrays in case the lengths documented for them aren't quite right?

[ oh, just had a look at your LPF - its not good enough, you need to ensure the stop band is adequate - perhaps
you are just hearing aliasing? I'd go for -80dB of stop band from Nyquist/4 to be nice and clean - think more like
40 taps than 16...
 
From what I can tell from the docs for the interpolation and decimation functions you are doing things right.

Have you tried without any of your intermediate processing, just interpolation straight back to decimation?

Have you tried adding some guard words to the state arrays in case the lengths documented for them aren't quite right?

[ oh, just had a look at your LPF - its not good enough, you need to ensure the stop band is adequate - perhaps
you are just hearing aliasing? I'd go for -80dB of stop band from Nyquist/4 to be nice and clean - think more like
40 taps than 16...

Thanks for your reply. Yes, I know the 16-tap filter is not great, but it's a whole lot better than linear interpolation. The 32-tap filter certainly is significantly better again, but it uses up quite a bit more processing (7% CPU rather than 5% per filter) which may be an issue in polyphonic applications. It's not a final filter yet, and maybe I'll add in multiple options in the final one.

But in any case, no the 'breaking up' is much more severe than aliasing, and occurs with even very low freq signals. It's a bit like the audio is rapidly switching between two different buffers, or buffer indexing is messed up when more than one filter is running. I'll try your suggestion of taking out intervening processing, but I can't see why it would work with just one filter running if that was the problem.
 
I moved the state array and instance declarations into the .h file as private, and I think that's sorted the problem out.
 
Last edited:
@rvh - is there a useable version of the code that I could try. I'd like to get a idea of latency and also to try with small block sizes. cheers Paul
 
@rvh - is there a useable version of the code that I could try. I'd like to get a idea of latency and also to try with small block sizes. cheers Paul

Paul, I'm still working on the FIR LPF for the polyphase approach, but hope to have it ready soon. At this stage I'm not thinking of making that filter minimum phase, and in a mock up polysynth I haven't noticed any latency issues. Does the version you currently have (without interpolation) have latency that bothers you? Note that changing the block size in the next version will require you to adjust some of the polyphase filtering parameters.
 
After further tsts, my current thoughts are to use 32 taps in the FIR after all. The CPU load goes up to ~7%, compared to ~5% for linear interpolation.
Below is the LPF I'm currently trying out (designed for 176kHz). Any comments are most welcome.

-0.000107320070740979531,
0.000747985193836334474,
0.002361583812865433,
0.004088989908423131,
0.004449967446630523,
0.001648422849822054,
-0.005257233263558728,
-0.014998389911544174,
-0.023380373153265226,
-0.023967831375363936,
-0.010344248425402229,
0.020767330578513626,
0.066936498872738684,
0.119469109852107777,
0.165552147905831670,
0.192517763075602144,
0.192517763075602144,
0.165552147905831670,
0.119469109852107777,
0.066936498872738684,
0.020767330578513626,
-0.010344248425402229,
-0.023967831375363936,
-0.023380373153265226,
-0.014998389911544174,
-0.005257233263558728,
0.001648422849822054,
0.004449967446630523,
0.004088989908423131,
0.002361583812865433,
0.000747985193836334474,
-0.000107320070740979531
 
Richard - no, interpolating version is fine.
I don't understand the maths well enough on the latest (FIR LPF) to understand what latency would be.
I'm interested in the filter for and effect unit (where it will be re-combined with a dry signal) so latency becomes an concern quite quickly to avoid unwanted comb effects
Cheers Paul.
 
After further tsts, my current thoughts are to use 32 taps in the FIR after all. The CPU load goes up to ~7%, compared to ~5% for linear interpolation.
Below is the LPF I'm currently trying out (designed for 176kHz). Any comments are most welcome.

Well the stopband is about 1/3rd of Nyquist, not 1/4th... But its about -60dB which is much better.

BTW I use Python and numpy / scipy.signal / matplotlib.pyplot packages for all my experimenting with
signal processing, makes this sort of thing very easy:
fir_interpolate.png
 
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:

FIR32.jpg
 
Back
Top