Skflneartrap paper

DSP, Plugin and Host development discussion.
RELATED
PRODUCTS

Post

Marvinh wrote: Sun Aug 18, 2024 6:32 am 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.
Suppose we have a linear two-pole low-pass (=all-pole in continuous-time) with a transfer function of 1/(s^2+s/Q+1).

Now, let's look at the response at three frequencies of interest:

At DC we have w=0. Substituting s=i*w=i*0=0, we have 1/(0^2+0/Q+1). So clearly this filter has unity gain at DC.

Since this filter has "normalized" cutoff at w=1 (we could set another cutoff by substituting s/wc for each s where wc is the desired cutoff, but working with normalized responses is often easier).

So at cutoff we have 1/(i^2 + i/Q + 1) where i^2=-1 (by definition), so this simplifies to 1/(i/Q) which is the same thing as Q/i. So at cutoff we have an amplitude gain of Q and 90 degrees phase-shift (from the imaginary unit in the denominator). This is the most important part: Q is essentially acting as a "gain control at cutoff" in terms of the amplitude response. The situation is a bit more complicated if our filter has zeroes as well as poles, since the "zeroes" also can have a "Q" value that's acting the opposite and giving inverse gain, but things kinda still follow the same idea (eg. a peaking EQ is literally poles and zeroes at the same frequency and the gain at cutoff is Qp/Qz where Qp and Qz are the Q-values of poles and zeroes respectively).

If we take w into the limit at infinity, then 1/(s^2+s/Q+1) = 1/inf = 0, which is why this is a lowpass.

The situation with higher order filters is a little bit more complicated, because while we can factor each filters into second order sections, the poles might be at different frequencies and how Q is actually defined can be chosen in a few ways.. but such filters can still be factored into a cascade of second order sections (each with their own cutoff and Q) and then the above applies to each such section again. If there's a "dominant" pole-pair that gives the resonance (eg. moog-style cascades) then it is mostly the Q value of that pole-pair that's most significant in terms of gain at cutoff. If we build a 4th order filter by stacking two 2nd order sections with the same cutoff and Q that we control directly, then the effective gain at cutoff becomes Q^2.

The thing with saturation in "musical" filters is that this linear behaviour is only allowed at rather low amplitudes, because we want to use a high resonance (=high Q), but we don't want it to get stupidly loud if the cutoff frequency is strongly excited. As soon as the resonance amplitude tries to grow higher, the non-linearities act (in a well-behaved filter) to dynamically lower the Q until the amplitude has fallen back to sane levels, at which point the tail of the ringing is allowed to continue at high-Q again. One can make this less obvious by driving the filter at much lower gain (at which point non-linear filters usually sound more like their linear counter-parts with annoyingly loud resonance), but the point is, this "dynamic Q" is a thing in many musical filters even if you don't drive them so hard that you'd actually start getting significant amounts of actual distortion (which.. we obviously often do, because it can also sound nice, but that's a different thing).

Note that "self-oscillation" is "infinite Q" which translates to "infinite amplitude at cutoff" yet often the self-oscillation amplitude of a musical filter might be less than the input level and that should give you an idea of how strongly the resonance is typically limited by non-linearities even when the filter still sounds like it has high resonance that's fairly "pure" in terms of distortion.

Post

I based my resonance two HP2 responses multipled together

0 == (1 + 2R s + (2 + R^2)s^2 + 2R s^3 + s^4)

Thats what I am hearing just straight pure clipping almost at high resonance then self oscillation to infinity.

To acheive dynamic Q i have to add the saturators?
edit nothing but regular resonator: lowpass2hipass2 pole
n = v0 == vin -tanh(2 R v1) - ((2+ R^2)v2) -tanh(2 R v3) - v4

Post

Have a read on the various KVR posts on regular 2 pole SVF single resonance non-linearities, you can apply similar methods here.
The Glue, The Drop, The Scream - www.cytomic.com

Post

Thank you guys for you help i did my original goal and applied the saturator to the Sallen Key Filter. It not the best but I am still messing around .

Edit: Self oscillation is resonance*(2 + M_E) solved with tangential solution then newton 4 X oversampling runs nice

Code: Select all

    class MSKF {
    public:
        double sampleRate;
        MSKF() {
            
        }
        
        void setSampleRate(double newSampleRate) {
            sampleRate = newSampleRate;
            onp.setParams(240, sampleRate);
        }
        
        onePoleLowpass onp;
        
        
        double ic1eq = 0.0;
        double ic2eq = 0.0;
        double ic3eq = 0.0;
        
        void filter(double cutoff, double resonance, float *in, float * out, int blockSize) {
            //WARNGING: do not pass nyquist unless you oversample
            //CUTOFF = HZ/(SAMPLE_RATE*OVERSAMPLING_FACTOR)
            //RESONANCE = 0.0 - 1.0
            //std max will prevent self oscillation
            
            
            double g = tan(M_PI * cutoff);

            double k = (2+M_E)*resonance;
            double k2 = 2*resonance;
            
            //smooth cutoff because analog circuits are not discrete
            
            
            g = onp.tickLinear(g);
            
            for (int i = 0; i < blockSize; i++) {
                
                g = onp.tickLinear(g);
                
                double tg = 2*g;
                
                double g2 = g*g;
                
                double gd = 1.0/(-1 - tg - g2 + g*k2);
                
                double v0 = in[i];
                
                double error = 1;
                
                double tol = 1e-6;
                
                double v1 = 0.0;
                
                double v2 =  0.0;
                
                double gk = 1.0;
                               
                double g2v0 = g2*v0;
                double gv0 = g*v0;
                
                //LINEAR VOUT ESTIMATE
                v2 = -((g*ic1eq + ic1eq + g*ic2eq + g2v0)*gd);
                v1 = -(ic1eq - g*ic1eq - ic2eq*k2 - gv0 - g2v0)*gd;
                
                //Tangential Estimate
                //Use linear resonance to set tangential sounds better because we dont have ic3eq yet
                
                double base = k2*v2;
                float tBase = tanh(base);
                
                double a = sech2_with_tanh( tBase );
                double b = tBase - base;
                
                v2 = -(g*ic1eq + ic1eq + g*ic2eq + g2v0)/(-1.0 - tg - g2 + g*(a+b));
                
                
                double ggk = 1.0;
                double gic1eq = g*ic1eq;
                double gic2eq = g*ic2eq;
                double next = 0.0;
                for(int j = 0 ; j < 50; j++) {
                
                    double x = k * v2;
                    double fk = tanh(x);
                    gk = 1 - (fk*fk);
                    ggk = g*gk;
                    ic3eq = fk - gk*x;
                    next = -((gic1eq + ic2eq + gic2eq + g*ic3eq + g2v0)/(-1-tg-g2+ggk));
                    error = std::abs(next-v2);
                    if(error < tol){
                        v2 = next;
                        break;
                    }
                    v2 = next;
                }
                
                
                double omg = (-1-g);
                
                v1 = -((ic2eq * gk - omg * (ic1eq + ic3eq + gv0))/(-(omg*omg) + ggk));
                
                ic1eq = 2 * (v1 - (gk*v2 + ic3eq)) - ic1eq;
                
                ic2eq = 2 * v2 - ic2eq;
                
                out[i] = v2;
            }
            
        }
        
    };
    
Last edited by Marvinh on Thu Aug 22, 2024 3:43 am, edited 5 times in total.

Post

One very simple question that I am almost embarrassed to ask as I am sure it is obvious - why use the trapezoidal rule rather than the midpoint rule, which is both simpler and has lower error?

I understand (or rather, I think I understand) that the appeal of this kind of ZDF filter design is that you think of this all as continuous time circuit analysis. Still, naively, having a stream of discrete input samples seems to make the midpoint rule very attractive as the samples are simply the midpoints and the integral approximation is simply a summation. I feel like I am either missing or getting something very fundamental here :lol:

Post

The Glue, The Drop, The Scream - www.cytomic.com

Post

Thank you very much, this one has been bugging me :)

Post

stoopicus wrote: Tue Aug 20, 2024 11:43 am One very simple question that I am almost embarrassed to ask as I am sure it is obvious - why use the trapezoidal rule rather than the midpoint rule, which is both simpler and has lower error?
The short version is that trapezoidal rule and implicit midpoint are equivalent in the linear case, the implementation cost is basically exactly the same (at least in the fixed time-step case) and the difference is basically the evaluation of non-linearities with half time-step shift one way or another (see section 4.2 of the paper Andy linked).

Practically speaking averaging the non-linear contribution from two points for each segments smooths them out a bit. That's perhaps "higher error" but it might or might not be better for aliasing... which can be further reduced if we further apply an approach such as ADAA.

But ... like realistically you can pretty much mix and match trapezoidal vs. midpoint as you see fit, whatever sounds better. There's little point discussing both 'cos once you optimize the implementation the difference really is rather tiny.

Post

Thank you, interesting.
mystran wrote: Tue Aug 20, 2024 1:18 pm But ... like realistically you can pretty much mix and match trapezoidal vs. midpoint as you see fit, whatever sounds better. There's little point discussing both 'cos once you optimize the implementation the difference really is rather tiny.
Yeah I was thinking it would save an addition and division per sample but that's only in the naive implementation in my head, I am sure optimizing it removes that advantage. (It probably is trivial to optimize to just one division total for example.)

Cool. This is really helpful for me, as simple as it sounds. The background is I am pretty new to audio processing but I have a lot of signal processing background from back in college - I was an analytical chem major and did my undergrad research in signal acquisition and processing. Back then I just approximated integrals with implicit midpoint and it worked great, so all the discussion of trapezoidal integrators has had me thinking "...but why?"

Post

stoopicus wrote: Tue Aug 20, 2024 7:44 pm Yeah I was thinking it would save an addition and division per sample but that's only in the naive implementation in my head, I am sure optimizing it removes that advantage. (It probably is trivial to optimize to just one division total for example.)
The text-book trapezoidal rule (for unit time-step) is y[n+1] = y[n] + (1/2)*(x[n] + x[n+1]), where x[n+1] is what we solve using whatever method such as LU + Newton.

There's no division. In practice, there's also kinda not a multiplication either, because the 0.5 factor can simply be merged into the tuning (or conductance in the circuit simulation case) coefficients, that is, we solve for (1/2)*x[n+1] directly (eg. by using tan(pi*f/fs) as tuning rather than 2*tan(pi*f/fs) which would be the text-book coefficient before averaging).

Finally, in the fixed time-step case we can write the integrator in "transposed direct form 2" form eliminating the second state variable that we'd need for x[n] in the text-book (= direct form 1) version:

y[n+1] = s[n] + (1/2)*x[n+1]
s[n+1] = y[n+1] + (1/2)*x[n+1]

That is, we add the contribution of (1/2)*x[n+1] that will be the (1/2)*x[n] term for the next time step already into the state variable as we update it.

The linear system we solve can also be written in terms of s[n] rather than y[n] so only s[n] needs to be preserved to the next time-step, so the most you will lose is perhaps one add per dimension... and that's assuming we won't need any extra multiplications to solve for a half time-step in the midpoint case (no idea if that can be optimized out; either way, the difference isn't exactly huge one way or another).

ps. In a sense we're solving for both at the same time: y[n] is the trapezoidal and s[n] is the mid-point shifted half time-step.

Post

Have any of you ever tried building a model of a differential equation in tensorflow and running the prediction in real time ? I know its mentioned in ADC video for analog modeling is it less taxing then newton with a good guess ? do state updates still need to happen?

Post

Figured out SVF! Suppose to use anti saturation. Not optimized of course.

Clear["Global`*"]

n1 = 0 == g (v0 - 2((R-1) v1+(gk v1+ic3eq)) - v2 ) - (v1 - ic1eq)
n2 = 0 == g(v1) - (v2-ic2eq)

s1 = Solve[{n1,n2},{v1,v2}]
FullSimplify[{s1}[[1]]]

Code: Select all

    struct SVF_ANTI_SAT {
        double ic1eq = 0;
        double ic2eq = 0;
        double ic3eq = 0;
        double cutoff = 40;
        double g = tan(cutoff * M_PI / 44100.0);
        double sampleRate = 44100.0;
        double fs = 1.0/44100.0;
        double R = 0;
        double impulse = 0;
        
        
        void setSampleRate(double newRate) {
            sampleRate = newRate;
            fs = 1.0/sampleRate;
        }
        
        void setCutoff(double newCutoff) {
            cutoff = newCutoff;
            g = tan(cutoff * M_PI / sampleRate);
        }
        
        void setRes(double res) {
            R = 1-res;
        }
        
        const double tol = 1e-8;
        
        double process(double v0) {
            
            double gk = 0.0;
            double v1 = 0.0;
            
            for(int j = 0 ; j < 50; j++) {
            
                double x = v1;
		// random diode clipper 
		//(1 - a/b) * x + a*sinh(x/b) = fk
		// gk is derivative
                double fk = (1-2/4.0)*x + (2)*sinh(x/4.0);
                gk = 1/2.0 * (1 + cosh(x/4.0));
                
                ic3eq = fk - gk*x;
                
                const double next = (ic1eq + g * (-ic2eq - 2 * ic3eq + v0)) / (1 + g*g + 2 * g * (-1 + gk + R));
                const double error = std::abs(next-v1);
                
                if(error < tol){
                    v1 = next;
                    break;
                }
                v1 = next;
            }
        
            const double v2 = (ic2eq + g * (ic1eq + 2 * ic2eq * (-1 + gk + R) + g * (-2 * ic3eq + v0))) / (1 + g*g + 2 * g * (-1 + gk + R));

            
            ic1eq = 2*v1 - ic1eq;
            ic2eq = 2*v2 - ic2eq;
            
            return v2;
        };
        
        void init() {
            ic1eq = 0;
            ic2eq = 0;
            ic3eq = 0;
        }
        
    };

Post

Marvinh wrote: Fri Sep 13, 2024 1:03 am Figured out SVF! Suppose to use anti saturation. Not optimized of course.
You can optimize this a lot simply replacing sinh with a cheaper "anti-saturation" function, because the exact conductance of the diode path is not all that important once you're above the "knee" and no real-world diodes would give you true sinh anyway (they necessarily have some series resistance).

In fact, this is a case where I'd argue a simple Taylor-expansion around zero (up to an order where the knee is hard enough for you) would perhaps be a sensible choice.

Post

Thanks for the info. I have been trying to figure out your pivotal method as well so I can avoid the the Newton Raphson loop is this a case as you mentioned would not work ? When I look up info about diodes it mentions PN junctions?

Post

Marvinh wrote: Fri Sep 13, 2024 6:17 pm Thanks for the info. I have been trying to figure out your pivotal method as well so I can avoid the the Newton Raphson loop is this a case as you mentioned would not work ? When I look up info about diodes it mentions PN junctions?
The method works fine for lumped diodes in SVF (eg. approximated by some polynomial).

Post Reply

Return to “DSP and Plugin Development”