SVF topology review

[catfact]

i'm taking another careful look at the topology for filter_svf after noticing some weird behavior, most prominent when given DC signals.

here's the supercollider patch generating the excitation signals. there are 5 varieties:

continuous white noise
white noise with 1ms exponential decay envelope
white noise with 0.1-1ms pulse rectangular envelope
DC with 1ms exponential decay envelope
DC with 0.1-1ms pulse rectangular envelope

if the rq is set to 0 and the filter is pinged with an envelope, sending it into oscillation, and the corner frequency subsequently swept, the frequency response is actually quite flat. the sine wave becomes distorted as CF approaches nyquist/2, as expected for this topology.

for RQ slightly >0, repeatedly sending a "decay DC" excitation has a predictable, smooth "ping". but repeatedly sending any of the others creates inconsistent amplitudes in the resulting sound. it is as if the system is not really LTI... feedback terms cancelling out new input. or something. i had noticed this before but attributed it to saturation during one or both of the summation steps.

however, softclipping isn't helping this particular issue (though it does i think improve the sound for typical uses cases when the input signal is hot.)

so i think you're right that it might be something totally wrong with the topology, event though it gives plausible results with most input. probably with the fixed-point implementation in particular. i'm going to try redoing it. i've come across a couple of improved methods since 2012 or whatever anyway.

[catfact]

ok well i looked into this yet again more closely, doing some "forensic" studies of what my original development process probably looked like.

i actually now believe the behavior is pretty much correct and as expected, but could of course be improved.

the algorithm i used is a classic digitization of chamberlin's 2-pole SVF published by jon dattorro back in 1997. the code can be found on ye olde musicdsp.org and the paper where he shows the derivation is still up on ccrma. (here is hal chamberlin's original schematic, by the way.)

the point of the derivation is to figure out the pole coefficients in terms of Q and freq such that those can be easily decoupled. dattoro goes about this by plugging the z-plane transfer function into mathematica, getting a simple equation for one coefficient and an incredibly huge horrible equation for the other, then applies approximation by mclaurin series, notices that some terms become vanishingly small in musical ranges, &c &c, and eventually (equations (28) snd (29) in the paper) comes up with a relatively simple new coefficient definition and transfer function.

to be honest, i implemented this years ago (long before the aleph) and ripped my own code. it seems like i tossed the "scaling" step, because it precludes self-oscillation. but maybe a compromise would be a better idea. you see, at high Q this thing has a huge amount of gain.... so for musical implementations, people often put some kind of nonlinearity in there (like our softclipper.)

back on musicdsp, you can see one of these attempts back in 2004 or something, by andrew simper from cytomic (http://www.musicdsp.org/showArchiveComment.php?ArchiveID=92). he adds a "drive" term that mixes in some band^3, and does a really basic 2x oversampling (same as ours; probably where i got the idea.)

in the last few years, there have been some new developments in analog filter modelling. it is becoming popular to use nodal analysis on the cirtuit, and then techniques like trapezoidal integration. andrew simper again (who is really quite smart with this stuff and makes excellent plugins) has a paper where he outlines the way to do this on the chamberlin SVF (and some other stuff.)

as an experiment, i implemented this newer approximation as a pd external: https://github.com/catfact/pd-externals/tree/master/svf

there's an example patch with white noise as input. so try it out and let me know what you think i guess... to me, it seems pretty similar in sound to the original, and has this same issue of being kind of "fluttery" with high Q.

the big benefit of this is that it has a much more natural-seeming resonance control. i believe it is a far closer approximation of the analog circuit.

the big drawback is that Q and frequency are no longer fully decoupled! several coefficients need to be recalculated each time either of them changes. the calculations are relatively simple and don't involve trig or square roots (it's not a biquad) but it does necessaitate a couple of divides on each parameter change.

but maybe it's worth porting to fixed point; i don't see a lot of issues with that, and it would be a good exercise to get our act together as far as providing an interface for efficient, arbitrary-radix operations (these need a couple of integer bits at least.)

so... i dunno. i'll put a little more time into this (it's come this far) but i'm not sure it's going to deliver massively superior results. it may be that some musical concerns would be better addressed by using a different type of filter entirely (sallen-key or cascaded 1-pole). and certainly by moving to the trapezoidal filter we would be adding some operations to the parameter change.

@electropourlaction , can you say something about how you derived your filter code? @rick-monster tagging you too

[ghost]

Haven't looked at your pd external yet, or the SC code - will do so bit later!

just to note one really easy optimisation to the non-interpolated oversampling (IIRC that's what's in there currently) is a simple linear-interpolating oversample:

fract32 run_two_frames (svf* s, fract32 next_sample) {
  fract32 halfway_sample = add_fr1x32(shr_fr1x32(next_sample, 1),
                                      shr_fr1x32(s->last_sample, 1));
  fract32 res = shr_fr1x32(run_one_frame(s, halfway_sample), 1);
  fract32 res = add_fr1x32(res, shr_fr1x32(run_one_frame(s, next_sample), 1));
  s->last_sample = next_sample;
}

providing an interface for efficient, arbitrary-radix operations

yeah I already have at least one half-baked, crappily named thing along these lines wafting round the new modules like eggy guffs, let's sort it out before bit-rot sets in. Any thoughts on how to treat fract division systematically? Is the answer that 'real' DSP engineers don't use divide!?

the big drawback is that Q and frequency are no longer fully decoupled!

Yeah, that's the magic of your current svf code! - pretty awesome hack I thought. Biquads are useful for static filters in modules, but totally agree with your analysis here. I used some random lisp code to calculate coefficients for voder's crossover then baked them right into the source code - not exactly 'maintainable'!

(sallen-key or cascaded 1-pole). and certainly by moving to the trapezoidal filter we would be adding some operations to the parameter change.

This sounds like something I'd like to learn - either of you know of any coefficient hacks for sallen key or gems from an old forum? Now correct me if I'm wrong, but isn't cascaded 1-pole pretty different beast from a 2nd order filter? If I understand right, a discontinuity in 1st derivatives will persist, no matter how many 1-pole filters you run in series? I guess at worst it's another thing to chuck in the toolkit...

Reminds me, I need to experiment with a 2-pole filter on the zero-crossing adaptive pitch detection. It currently uses 1-pole lowpass in series with 1-pole highpass, which is all I knew/understood earlier this year!

(like our softclipper.)

Btw I know a way better softclipper design using hyperbolic fns, just not quite sure how to implement on blackfin yet...

[electropourlaction]

I have just tried the new pd-example, by calculating the g-parameters (floats) on the avr32 and just doing the update function on the bfin, here is my naive test!

scene->chn[n].a is cutoff and scene->chn[n].b is resonance, n is just the channel number.

CUTOFF tmp = scene->chn[n].a; tmp += val; if (tmp < 0) tmp = 0; scene->chn[n].a = tmp; scene->g[n] = tan(M_PI * tmp / 48000); scene->g1[n] = scene->g[n] / (1.0 + scene->g[n] * (scene->g[n] + scene->chn[n].b)); scene->g2[n] = 2.0 * (scene->g[n] + scene->chn[n].b) * scene->g1[n]; scene->g3[n] = scene->g[n] * scene->g1[n]; scene->g4[n] = 2.0 * scene->g1[n]; ctl_param_change(n, eParamG, float_to_fr32(scene->g[n])); ctl_param_change(n, eParamG1, float_to_fr32(scene->g1[n])); ctl_param_change(n, eParamG2, float_to_fr32(scene->g2[n])); ctl_param_change(n, eParamG3, float_to_fr32(scene->g3[n])); ctl_param_change(n, eParamG4, float_to_fr32(scene->g4[n]));

Q tmp = scene->chn[n].b; tmp += val; if (tmp < 0) tmp = 0; scene->chn[n].b = tmp; scene->g[n] = tan(M_PI * scene->chn[n].a / 48000); scene->g1[n] = scene->g[n] / (1.0 + scene->g[n] * (scene->g[n] + tmp)); scene->g2[n] = 2.0 * (scene->g[n] + tmp) * scene->g1[n]; scene->g3[n] = scene->g[n] * scene->g1[n]; scene->g4[n] = 2.0 * scene->g1[n]; ctl_param_change(n, eParamG, float_to_fr32(scene->g[n])); ctl_param_change(n, eParamG1, float_to_fr32(scene->g1[n])); ctl_param_change(n, eParamG2, float_to_fr32(scene->g2[n])); ctl_param_change(n, eParamG3, float_to_fr32(scene->g3[n])); ctl_param_change(n, eParamG4, float_to_fr32(scene->g4[n]));

BFIN (re-named the v-parameters to c-something just because I already was naming that kind of parameters in that way..) c-> is the selected channel, only output is the bandpass filter.

c->c0 = c->input(c); c->c1z = c->c1; c->c2z = c->c2; c->c3 = sub_fr1x32(add_fr1x32(c->c0, c->c0z), mult_fr1x32x32(0x7fffffff, c->c2z)); c->c1 += sub_fr1x32(mult_fr1x32x32(c->g1, c->c3), mult_fr1x32x32(c->g2, c->c1z)); c->c2 += add_fr1x32(mult_fr1x32x32(c->g3, c->c3), mult_fr1x32x32(c->g4, c->c1z)); c->c0z = c->c0; return c->c1;

some trouble with the resonance parameter in that the volume drops significantly as I increase it, probably an error somewhere on my side!

the big question though is how to modulate this on the bfin-side, I think that using that tan-function on the bfin would really slow things down, but maybe there is a fast fixed-point bfin-library function that can be used instead?

sounds promising so far though, will make some more tests and try to find the cause of the volume drop.

[catfact]

c->c3 = sub_fr1x32(add_fr1x32(c->c0, c->c0z), mult_fr1x32x32(0x7fffffff, c->c2z));

that line seems not quite right? mult_fr1x32x32(0x7fffffff, x) returns x unchanged? maybe should be somrthing like c->c3 = sub_fr1x32(add_fr1x32(c->c0, c->c0z), shl_fr1x32(c->c2z, 1)); and of course make sure the input is always small. maybe right shift 1 or 2 places and left shift the filter output. or use a different radix for all the calculations...

another thing that might be worth trying is just using the blackfin's fast-float library. there are intrinsics for non-IEEE floats that use 16b mantissas and should actually be pretty fast. not as fast as fixed-point but certainly doable if you're not trying to do a whole lot of stuff. here's an appnote

the appnote uses VDSP intrinsics, but they exist in the gcc toolchain too under slightly different names e.g. fract16_to_ff16 in VDSP is fr16_to_fl16 ... place to look is bfin-elf/include/float16.h in your toolchain directory.

[catfact]

isn't cascaded 1-pole pretty different beast from a 2nd order filter?

ach, sorry, my bad. i meant cascaded first-order, e.g. a ladder filter. i guess those stages have 1 pole and 1 zero.

know of any coefficient hacks for sallen key or gems from an old forum?

not really. sallen-key is popular generic topology in analog because it is so very simple. but i haven't seen a lot of digital implementations, except (again) andrew simper's at cytomic, doing a nodal analysis model.

both the ladder and the sallen-key require a tan in the coefficient calculation, just like the nodal-analysis SVF above... regarding that...

I think that using that tan-function on the bfin would really slow things down, but maybe there is a fast fixed-point bfin-library function that can be used instead?

there are fast_tan functions for floats in the toolchain. but yeah, that is probably not a good way to go. i would use a lookup table for frequency on the avr32 side, just as i did with the other implementation.

[electropourlaction]

shl_fr1x32(c->c2z, 1)

aha, that is how to multiply with 2, makes total sense now that I see it in written, thanks! lots of ideas in this thread now, need time to digest!

[ghost]

ach, sorry, my bad. i meant cascaded first-order, e.g. a ladder filter. i guess those stages have 1 pole and 1 zero.

Now I'm a little confused (don't understand the poles/z-transform stuff well at all). So just to confirm we're talking the same language - 1st order filter means 1 state variable (e.g exponential decay, simple slewing filters, RC voltage divider). 2nd order means 2 state variables (e.g damped simple harmonic oscillator, svf/biquad filters, RLC network). Pretty sure that's correct!? I'm not gonna try and talk poles till I've sat down & wrapped my head around it...

Kind of surprised to hear the 'ladder filter' topology is effectively cascaded 1st order (Just pored over the circuits briefly for diode/transistor ladder filters, didn't really get how they work) - so as I understand your description of ladder, the lowpass version never has a resonance bump before tailing off, but as you add more & more stages, it becomes steeper dropoff going above the break frequency?

Then reading through this (http://www.timstinchcombe.co.uk/synth/Moog_ladder_tf.pdf) a bit more, seems you can achieve variable Q by feedback round the whole loop. I guess this works because the phase shift is zero at break frequency for each 1st order stage?

Feel like I just learned something new - hopefully that doesn't sound like total bs!?

If the LTI version works, I guess it's easy to just softclip a bit at different stages in this cascaded 1st order architecture, then we've got another flavour of filter to choose from...

[catfact]

1st order filter means 1 state variable... 2nd order means 2 state variables

not exactly, it's the number of poles and/or zeros, which corresponds to the order of the polynomial in the numerator (for zeros) or denominator (for poles) of the transfer function... or something like that.

in discrete signal processing terms, that translates to how many past samples we end up using.

without getting too deep into z-transforms (which i'm rusty on anyways,) we can say in a (very small) nutshell that feed-forward terms create zeros (that is, mixing in old input) and feedback terms create poles (mixing in old output.)

the SVF we're using so far is 2nd order (there are two samples of delay, though it's not immediately obvious from just looking at the equations), and it's all poles (delays always come after summation.)

now that i'm revisiting at papers on discretization of the moog vcf, i see that there are all kinds of ways of dealing with it. the continuous-time transfer function is all-pole, but with delay-free feedback. anyways this paper (http://www.acoustics.ed.ac.uk/wp-content/uploads/AMT_MSc_FinalProjects/2012__Daly__AMT_MSc_FinalProject_MoogVCF.pdf) is a good rundown of different approaches, with bibliography... also worth checking out the implementations in supercollider, which is relatively easy to follow, and csound, which... uh... isn't.

variable Q by feedback round the whole loop

yes, that's exactly how resonance is added in a ladder filter. the cascade is just to sharpen the rolloff...

anyways, back to my assigned task! which is issue 270. i'm "unassigning" myself from this one for now.

[ghost]

ok well - this should be pretty much the ladder filter discussed above:

https://github.com/rick-monster/aleph/blob/dev/dsp/filter_ladder.c

Felt like the feedback didn't really work until I added a DC blocker in feedback path. Whole thing just sounds a bit wrong... Gonna try and check that MSc project, also try doing scaling multiplies see if that makes a difference to stabiltiy/character...

[ghost]

OK just read through Daly's thesis - looking at his difference equations it seems that there are two things totally different in filter_ladder.c:

1. Houvilainen applies negative feedback to achieve Q control, whereas I'd randomly reached for +ve feedback - probably hence the need for DC blocker...
2. Houvilainen optimised his algorithm by running the feedback with an extra half sample delay (linear interpolated)
3. The difference equations can be written a bit more neatly, requiring only a single multiply per stage
4. Seems there's some value in over-sampling for a clean high-frequency response (gonna try the linear-interpolated samplerate doubler)

So yea - should be easy enough to nail down the linear version tonight, then try chucking our existing softclipper in there. Also going to review the lookup table implementation in waves, get a better sense how CPU intensive to implement tanh soft-clipping (I know a tanh spell allowing +ve/-ve limit/knee controls, which seems pretty optimal to me...)

[ghost]

Ha - that turned out not difficult at all. Just pushed oversampling, clipping & oversampling clipping versions of the filter to my dev branch. Sounds like a kind of squelchy, resonant digital low-pass to these ears!

This ladder filter topology is pretty understandable/hack-able, gonna try and come up with some extra twists... to finish off the filter_ladder class:

• asym clipping version (oversampling and non-oversampling)
• tweakable clipping coefficients at different filter stages
• variable number of stages (restrict to even multiples so breakpoint attenuation scales in factors of 2)
• bandpass & highpass filter versions (bandpass can be achieved by alternating highpass/lowpass stages

[electropourlaction]

is Q before or after the output in our original SVF?

[electropourlaction]

probably not any better but this sounds ok to me!

//  parameter a sets cutoff frequency with cv modulation
//  parameter b sets feedback
fract32 mod;

filter_1p_lo_in(&(c->pSlew), c->a);

if (filter_1p_sync(&(c->pSlew)))
{
    mod = add_fr1x32(c->a, mult_fr1x32x32(cv[c->amod], c->amodl));
}
else
{
    mod = add_fr1x32(filter_1p_lo_next(&(c->pSlew)), mult_fr1x32x32(cv[c->amod], c->amodl));
}

c->c0 = add_fr1x32(c->c0, mult_fr1x32x32(mod, c->c2));
c->c1 = sub_fr1x32(sub_fr1x32(negate_fr1x32(c->input(c)), c->c2), c->c0);
return c->c2 = add_fr1x32(mult_fr1x32x32(c->b, c->c2), mult_fr1x32x32(mod, c->c1));