Skflneartrap paper

[Marvinh]

https://cytomic.com/files/dsp/SkfLinear ... mised2.pdf

I figured out a 6db hi pass from this filter but i am wondering which part of the equation describes the feedback loop in order to try out different saturators? Is it literally just saturate the input?

something like in = 2*tanh(in*drive/2);

Code: Select all

    void setParams(double newSampleRate, double newCutoff, double newRes) {
        sampleRate = newSampleRate;
        cutoff = newCutoff;
        g = tan(M_PI*cutoff/sampleRate);
        res = newRes;
        k = 2*res;
        a0 = 1.0/(pow((1.0 + g),2.0) - g*k);
        a1 = k*a0;
        a2 = (1 + g)*a0;
        a3 = g*a2;
        a4 = g*a0;
        a5 = g*a4;
    }
    
    double process(double in) {
        
        v1 = a1*ic2eq + a2*ic1eq + a3*in;
        v2 = a2*ic2eq + a4*ic1eq + a5*in;
        hi = in-v1;
        ic1eq = 2*(v1 - k*v2) - ic1eq;
        ic2eq = 2*(v2) - ic2eq;
        low = v2;
        return low;
    }

[mystran]

https://cytomic.com/files/dsp/SkfLinear ... mised2.pdf

I figured out a 6db hi pass from this filter but i am wondering which part of the equation describes the feedback loop in order to try out different saturators? Is it literally just saturate the input?

something like in = 2*tanh(in*drive/2);

Essentially the way this filter works, you have two one-pole low-pass stages on the top of Andy's beautiful drawing.

Andy's description is perhaps a bit confusing, 'cos it doesn't really explain the derivation too well, but basically each one-pole lowpass is an OTA (operational transconductance amplifier; it outputs a current proportional to the control current times the difference between the input voltages) which feeds a capacitor which integrates the current, which is then buffered (so we don't draw current from the integration capacitor) and the output is fed back into the negative input of the OTA, so the integrator becomes a low-pass filter instead.

To add resonance, we take the output of the 2nd stage, multiply it by the resonance control "k" and feed it back to "ground side" of the first integration capacitor. What this effectively does is that as far as the feedback is concerned, the capacitor (together with the local feedback loop around the OTA) is acting as a high-pass filter instead. The 2nd stage then acts as low-pass as usual in terms of the signal coming from the 1st stage, so the overall response around the feedback loop is a bandpass, which is what we need for a 2-pole filter.

The usual way to add saturation (eg. to allow self-oscillation, which is usually the fun part of these filters) is to saturate the feedback path after multiplication by "k." In an MS20 rev2 for example, the "k" is a potentiometer and saturation is done by a buffer with diodes in negative feedback which essentially gives a+tanh(b*x) for some values of a,b that depend on the diodes). If I recall correctly, in MS20 rev1 (which also uses slightly different kind of one-poles not based on OTAs) the saturation is done for the output signal before the feedback is taken, so that's another possibility.

However.. here's the catch: becomes this is a trapezoidal solver, you can't just add the tanh() as such somewhere... 'cos it's part of the feedback loop that needs to be solved to eliminate any delays. That's why trapezoidal filters are sometimes called "zero delay filters" 'cos we actually build them out of differential equations which we integrate using the trapezoidal rule and trapezoidal rule is an implicit method, so this essentially requires solving a simultaneous set of equations (=matrix inversion, but before you panic that's exactly what Andy's code is already doing).

So your choices are basically to either use Newton or cheat. I won't go into detail of how to do Newton here ('cos explaining it well is a bit too much effort; the short version is that you "linearize" the system by fitting a tangent-hyperplane and then iterate until said plane no longer changes), but basically to do the MS20-style thing "the cheap way" you can replace "k" with a time-varying gain-term (see the thread I linked; the example filter is a transistor ladder, but same basic idea) that depends on the output previous sample... but keep in mind that you'll then need to also recompute the coefficients that depend on k every sample (so the way Andy has optimized the linear case might no longer be ideal; not sure to be honest).

[Marvinh]

thank you for the info I will try your filter and try to replace with ota. Just one question.

Is pole mixing still zdf? If I have a bunch of one pole and just rearrange them.

a = lp.process(in);
b = lp.process(a);

a + b;


this seems like a good way for newton raphson afterwards
https://urs.silvrback.com/one-pole-monster

[mystran]

thank you for the info I will try your filter and try to replace with ota. Just one question.

Rather try to understand the technique. This method (trapezoidal integration in general) can be bit of a mind-bender at first, because you really don't design anything in code, you design it mathematically and then you write a solver... and it's that solver that becomes your code.

We can literally put any of these filters into matrix form and just use a generic ODE solver (eg. LU solver with Newton-Raphson to linearize; that's what Spice does.. theoretically anyway). This works perfectly fine, but it's slow... so instead it makes sense to essentially do partial evaluation of the generic solver and just write the actually useful steps into straight-line code for compiler to optimize properly, but there's no magic in the code because it's just solving a system of equations.

All the magic happens in the equations.

Is pole mixing still zdf? If I have a bunch of one pole and just rearrange them.

So called "ZDF" is for all intents and purposes just trapezoidal integration. Pole-mixing is valid in analog and trapezoidal preserves that. In fact since trapezoidal (at least once you "prewarp" the cutoff, but that's where the tan(pi*f/fs) comes from) gives you bilinear transform, which maps the frequency response from the analog domain "as-is" except for warping of the frequency axis.

In fact, most of the time if you're fine with the frequency warping (for non-linear filters that's usually not your biggest problem, 'cos mitigating aliasing generally requires some oversampling anyway), you can design your pole mixing in analog and use the same coefficients with a trapezoidal filter and everything will work as expected.

this seems like a good way for newton raphson afterwards
https://urs.silvrback.com/one-pole-monster

There's a few (mathematically equivalent) ways of doing it, that's one way... but the key to this whole trapezoidal filter business is that the exact code used to solve the system no longer matters, it's not what defines the filter. You have a bunch of differential equations, you need a linear solver and a way to linearize those equations... but how you choose to do that mostly only matters as far as some approaches are faster and some approaches are numerically better behaved.

[Marvinh]

For some reason i thought newton raphson was almost always necessary for some circuits so i avoided it for a while until I tried this ms 20 clone https://www.threetom.com/products/steves-ms22/ wanted to get that feel into a project im working on.

thank you for your help.

[mystran]

Well.. Newton (or some other iterative scheme) is typically necessary for non-linear circuits if you want accurate results. With non-iterative "short-cuts" you're basically making a "best guess" and then you live it, whatever the accuracy. How well this works depends a whole lot on both on the circuit and the nature of the non-linearities. Simple "lumped" saturation of some signal (eg. a feedback path) typically works well enough with the one-pass origin-pivot for example, while something like individual pn-junctions really do require Newton to get even ball-park reasonable results.

The beautiful thing about Newton though is that once you can get your estimate close enough to get into the quadratic convergence region, precision (ie. number of correct digits) doubles every iteration, so with Newton the first goal is to make sure it always converges (which is not always given) and then try to make an initial guess every iteration that is good enough to get into the quadratic convergence quickly. Fortunately, most of the time (especially with oversampling to mitigate aliasing; note that once you get significant aliasing, your simulation results aren't going to be all that accurate anyway) audio signals are actually quite "slowly moving" in terms of voltages from one iteration to the next, so even just taking the previous time-step operating point as the "initial guess" for Newton often allows a quick convergence (even just an iteration or two) for most time-steps. Some time-steps might take a lot more (eg. 20+) but as long as the average stays low, it's quite workable.

[andy-cytomic]

If you want a rundown on the non-linearity placement and how to solve for them have a look at my ADC 2020 slides and video which are the first links on my technical papers page:

https://cytomic.com/technical-papers/

I show both MNA and State Space solving, as well as how to linearise the non-linearities and then iterate to a solution, and also how to pre-compute solutions via the DK method. You can apply the same method I show how to solve the diode clipper to solving the SKF non-linear filter. As Mystran has suggested I would first use a single non-linearity around the feedback "k" term, so going on the equations on page 63 we have:

Code: Select all

0 = -g.(v1 - v2) + gc1.(v2 - f(k.v3)) - ic1eq
0 = -g.(v2 - v3) + gc2.(v3 - 0) - ic2eq

where f(x) is whatever non-linear function you want, eg f(x) = tanh(x)

Code: Select all

x = k.v3
fk = f(x)
gk = f'(x) (derivative)
gkeq = fk - gk.x

then solve the following, and iterate to a solution by re-evaluating the above terms each time:

Code: Select all

0 = -g.(v1 - v2) + gc1.(v2 - (gk.v3 + gkeq)) - ic1eq
0 = -g.(v2 - v3) + gc2.(v3 - 0) - ic2eq

[mystran]

As Mystran has suggested I would first use a single non-linearity around the feedback "k" term, so going on the equations on page 63 we have:
[...]
where f(x) is whatever non-linear function you want, eg f(x) = tanh(x)

I want to elaborate a little bit more on the non-linearity itself in MS20-style designs.

Ignoring the diodes, the opamp is in an elementary non-inverting amplifier configuration, with the resistor divider from the output to the negative input setting the gain. I can't seem to find original schematic for rev2 right now, but clones and original rev1 all generally choose resistor values such that there's roughtly 5.5:1 ratio (eg. 10k:1.8k, or 270k:47k) with feedback to ground (neglecting whatever capacitors might be on the path), which gives a nominal gain of about 6.5.

That's basically the maximum feedback for small signals with resonance cranked and puts self-oscillation (gain around 4) at around 2/3rds of the range, so the chosen ratio actually makes sense. Now, going back to the diodes, for large signals these bypass the feedback resistor, essentially turning the opamp into a unity-gain non-inverting buffer. So this isn't truly a clipper at all, it just drops the "small-signal" gain of the feedback down to unity (ie. significantly below self-oscillation) when it exceeds the diode threshold.

So, assuming symmetrical diodes (ie. lumped iD=sinh(vD)), the feedback function is well-approximated by f(x)=x+5.5*(tanh(x)) which applies the gain of 6.5 for low signal levels and gain of 1 for high signal levels, so external resonance feedback should now be [0,1].

This is obviously approximate... but my main point is that it isn't a true clipper in MS-20-style designs, rather it just clips "most of" the resonance.

[Marvinh]

Nvm I think i am mistaken LP2 is still a hi pass i have to research some more it might be a 4 pole HP but i doubt it. I have to plot this out


Really excited about this.
I think I was able to piece together a 4 pole SVF LowPass with some help from you guys and vadims book.
Its super fast too as long as you process in blocks of 16 or so

g = tan(pi*cutoff/samplerate)
R = 1-r;
LP2 = (x - (2g^2 R s1 + 2R s1 + g s1 + s2)) / (1 + 2g^3 R + 2g^2 R + 2gR + 2R + g^2);
v1 = g*lp2;
BP = v1+s1;
s1 = BP + v1;
v2 = g*BP;
lp4 = v2+s2;
s2 = lp4 + v2;

SVF 2 pole but instead of taking the resonance from the hi pass take it from the output of the two pole low pass feed it in to another SVF 2 pole

this is original equation I put together :

LP2 = x - (g(g(2R(g LP2 + s1))) + g(2 R(g LP2 + s1)) + 2 R(g LP2 + s1) + g(g LP2 + s1) + s2)

[andy-cytomic]

Here is a reference for an analog 4 pole SVF for use as a crossover:
https://jahonen.kapsi.fi/Audio/Papers/Statevariable.pdf

From that the SVF4 has these Laplace equations, where d the damping factor = 1/Q, the "Quality" factor:

Code: Select all

v0 == vin - 2 d v1 - (2 + d^2) v2 - 2 d v3 - v4
0 ==  -g v0 + s v1
0 ==  -g v1 + s v2
0 ==  -g v2 + s v3
0 ==  -g v3 + s v4

v0 is hp4
v1 is lp1 hp3
v2 is lp2 hp2
v3 is lp3 hp1
v4 is hp4

[Marvinh]

n1 = 0 == vin - 2 d v1 - (2 + d^2) v2 - 2 d v3 - v4 - v0;
n2 = 0 == -g v0 + s v1;
n3 =0 == -g v1 + s v2;
n4 = 0 == -g v2 + s v3;
n5 = 0 == -g v3 + s v4;
Solve[{n1,n2,n3,n4,n5},{v0,v1,v2,v3,v4}]
edit: no s in the solver

s is the same for all equations just one state? Is this the correct way to input into mathematica?

[andy-cytomic]

s is the laplace variable, which is s = i w, or s = i 2 pi fc, where i*i = -1

Have a read of my technical papers on how to do numerical integration to turn the s into a gc and iceq using trapezoidal integration:

https://cytomic.com/technical-papers

[Marvinh]

it works but not sure why but happy with it. im not so certain about the state updates i just followed Art of Va f0 is a super HighPass f4 is Low Pass.
everything else seems like Band Pass

edit: it resonates a little to loudly i am wondering if this is a reason for a saturator? I know we can add smoothers to the g parameter.

Code: Select all

n = v0 == vin -2 R (g v0+s0)- (2+ R^2)(g (g v0 +s0)+s1) -2 R(g(g (g v0 +s0)+s1)+s2) - (g(g(g (g v0 +s0)+s1)+s2)+s3)

Code: Select all

    struct MSVF {
        double s0;
        double s1;
        double s2;
        double s3;

        MSVF() : s0(0.0), s1(0.0), s2(0.0), s3(0.0) {}

        void filter(double cutoff, double resonance, float *in, float *out, unsigned nsamples) {
            //WARNGING: do not pass nyquist unless you oversamples
            //CUTOFF = HZ/(SAMPLE_RATE*OVERSAMPLING_FACTOR)
            //RESONANCE = 0.0 - 1.0
            //std max will prevent self oscillation to inifity
            
            double g = tan(M_PI * cutoff);
            double g2 = g * g;
            double g3 = g2 * g;
            double R = std::max(2 - 2 * resonance, 0.01);
            double R2 = R * R;
            double denominator = std::pow(1 + g2 + g * R, 2);

            for (unsigned i = 0; i < nsamples; i++) {
                double vin = in[i];

                
                double g_s0 = g * s0;
                double g_s1 = g * s1;
                double g_s2 = g * s2;
                double g_s3 = g * s3;
                double pow_g3_s0 = g3 * s0;
                double pow_g2_R_s0 = g2 * R * s0;
                double pow_g2_s1 = g2 * s1;
                double pow_R2_s1 = R2 * s1;
                double pow_R2_s2 = R2 * s2;

                
                double numerator = -2 * g_s0
                                   - pow_g3_s0
                                   - 2 * R * s0
                                   - 2 * pow_g2_R_s0
                                   - g * R2 * s0
                                   - 2 * s1
                                   - pow_g2_s1
                                   - 2 * g * R * s1
                                   - pow_R2_s1
                                   - g * s2
                                   - 2 * R * s2
                                   - s3
                                   + vin;

                
                double y = numerator / denominator;
                double f0 = y;       //HIPASS 4
                double v1 = g * f0;
                double f1 = v1 + s0; //BP
                s0 = f1 + v1;

                double v2 = g * f1;
                double f2 = v2 + s1; //BP
                s1 = f2 + v2;

                double v3 = g * f2;
                double f3 = v3 + s2; //BP
                s2 = f3 + v3;

                double v4 = g * f3;
                double f4 = v4 + s3; //LOWPASS 4
                s3 = f4 + v4;

                
                out[i] = f4;
            }
        }
    };

[Marvinh]

s is the laplace variable, which is s = i w, or s = i 2 pi fc, where i*i = -1

Have a read of my technical papers on how to do numerical integration to turn the s into a gc and iceq using trapezoidal integration:

https://cytomic.com/technical-papers

https://cytomic.com/files/dsp/OnePoleLinearLowPass.pdf the numerical integration is here correct?

this is the correct idea or its different for almost everything just chaining like this ? when it comes to any type of filter

n = y == x - (2/Q) (g y+s0) - ((1+2Q^2)/Q^2) (g(g y+s0)+s1) - (2 /Q)(g(g(g y+s0)+s1)+s2) - (g(g(g(g y+s0)+s1)+s2) + s3)

[andy-cytomic]

You can use whichever integration scheme and formulation you prefer, normally trapezoidal and nodal / state space are used.

Also the general form of a 4 pole state variable filter has v0 as:

v0 = vin - d1 v1 - d2 v2 - d3 v3 - v4

so feel free to change d1 to d3 to whatever you prefer, I just provided you a link with a particular type of damping arrangement that was used from the paper I linked to, you can use whatever you want. I would recommend solving for values of d1 to d3 that place one pair of complex conjugate poles on the imaginary axis in Laplace space, as these are the ones that result in self oscillation.