Aliasing suppression in filter feedback loop in a ZFD filter

[xhy3]

Hi together

I like to suppress aliasing within a zfd filter in the feedback loop. I have found following interesting paper about the topic:

https://www.researchgate.net/publicatio ... onvolution

The solution looks promising and easy to implement. The only problem is the delay the filter kernel of this method introduces. I see fig. 9 in the paper, but i'm not sure how to implement that and had no success so far. I also want the tanh only for the feedback signal and not for the whole filter pole when possible.

Here is the filter code:

Code: Select all

m_wcD = 2.0f * fasttan(m_T * Pi * freq) * sampleRate;
float TwcD = m_T * m_wcD;
m_b0 = (TwcD / (TwcD + 2.0f));
m_a1 = ((TwcD - 2.0f) / (TwcD + 2.0f));

k = resonance * 4.0f; 

oL = m_s4L + m_b0 * (m_s3L + m_b0 * (m_s2L + m_b0 * m_s1L));
const float b4 = m_b0 * m_b0 * m_b0 * m_b0;
const float divisor = 1.0f / (1.0f + b4 * k);
outputL = (b4 * inputWithoutRefL + oL) * divisor;

float clippedL = k * outputL;
antialiasedClipper.processTanh(&clippedL); // introduces delay (rect filter kernel, maybe 0.5 sample delay)
        
uL = inputWithoutRefL - clippedL;

float future1L = (m_b0 * uL + m_s1L);
m_s1L = (m_b0 * uL) - (m_a1 * future1L);

float future2L = (m_b0 * future1L + m_s2L);
m_s2L = (m_b0 * future1L) - (m_a1 * future2L);

float future3L = (m_b0 * future2L + m_s3L);
m_s3L = (m_b0 * future2L) - (m_a1 * future3L);

m_s4L = (m_b0 * future3L) - (m_a1 * outputL);

Would be great if someone can help here. Thanks!

[stratum]

The link doesn't work

[Miles1981]

The link doesn't work

It works in my browser.

[xhy3]

Sorry for the slow progress. I don't have much time for this at the moment...

For simplifying things i take the zdf filter posted by mystran (viewtopic.php?t=368466) without the non linearities:

Code: Select all

        float input = *inputL;
        if (cutoff > 1.0f) cutoff = 1.0f;
        if (cutoff < 0.0f) cutoff = 0.0f;
        float f = tan(M_PI * cutoff * 0.5f);
        float r = 4 * resonance;
        
        float g = 1 / (1 + f);
        float f3 = f*g, f2 = f*g*f3, f1 = f*g*f2, f0 = f*g*f1;

        float y3 = (g*s[3] + f3*g*s[2] + f2*g*s[1] + f1*g*s[0] + f0*input) / (1 + r*f0);
        
        float xx = input - r*y3;
        float y0 = g*(s[0] + f*xx);
        float y1 = g*(s[1] + f*y0);
        float y2 = g*(s[2] + f*y1);
        
        s[0] += 2*f * (xx - y0);
        s[1] += 2*f * (y0 - y1);
        s[2] += 2*f * (y1 - y2);
        s[3] += 2*f * (y2 - y3);

The paper says we can fold the antialiased non-linearity into the first pole. I get this when i do that:

Code: Select all

        float y0 = g*(log(cosh(f*(in - r*y3))) - log(cosh(s[0])))/(f*(in - r*y3) - s[0]) * 2;
        float y1 = g*(s[1] + f*y0);
        float y2 = g*(s[2] + f*y1);

The problem is that maxima shows no result for this when i try to solve it for y3.
Looks like there is no solution. Any input is welcome. I keep trying to make it work...

After that i also have to take a look into the x = x-1 case.

edit:
I can solve the equation when i move the y3 outside of the nonlinearity, but the filter resonance is unstable this way and it does not help much, because i need especially the feedback part to limit.

Code: Select all

float y0 = g*2*((log(cosh(f*input)) - log(cosh(s[1])))/(f*input - s[1]) - f*k*y3);

Not sure if i also have to modify the integrator part:

Code: Select all

        s[0] += 2*f * ((input - k*y3) - y0);

[xhy3]

I can't find any solution. Anyone an idea how to fold the antialising filter into the first filter pole and eliminate the delay?

Edit: Anyone already used a Jacobian matrix to solve non linear equations?

[Smashed Transistors]

Maybe you can use some step invariant transform for the linear "ZDF" part and combine it with DPW (also known as 'antiderivative anti aliasing' or 'pre-integrated wavetables'). The maths should be fun to play with

[xhy3]

Maybe you can use some step invariant transform for the linear "ZDF" part and combine it with DPW (also known as 'antiderivative anti aliasing' or 'pre-integrated wavetables'). The maths should be fun to play with

Interesting idea. This could work. I also think i should consider other solutions like the one you mentioned. BLAMPs maybe work too.

[Smashed Transistors]

Another idea that can limit distorsion induced by aliasing is to consider the non linearity as a non linear gain/integration factor.

I have combined it to the step invariant method in a 2 pole SVF filter JSFX plugin (for Reaper) :

post:
https://forum.cockos.com/showpost.php?s ... stcount=18

jsfx thread:
https://forum.cockos.com/showthread.php?t=166480

and it sounded nice to my ears.

[xhy3]

Thanks for helping. I will have a look at the step invariant method. Looks interesting.
In the meantime i have used mystran's idea for folding in the nonlinearities and it seems to work well. It's not perfect, because i do nothing about he 0.5 sample delay of the antiderivate, but it's still stable and does a good job in antialiasing.

You maybe have to correct the instabilities the filter has with full resonance at higher frequencies and i would use it with 2x oversampling.

Code: Select all

        // tuning and feedback
        float input = *inputL;
        float f = tan(M_PI * cutoff * 0.5f);
        float r = 5.5 * resonance;
        
        // delay input
        float ih = 0.5 * (input + z); z = input;
        float ihf = ih - r * s[3];
        
        float s0 = s[0];
        float t0 = antialiasedClipper.processTanh(ihf, s0);

       // tanh(x) / x
        t0 /= 0.5 * (ihf + s0);
        
        // g# the denominators for solutions of individual stages
        float g = 1 / (1 + f);
        
        // f# are just factored out of the feedback solution
        float f3 = f*g, f2 = f*g*f3, f1 = f*g*f2, f0 = t0*f*g*f1;
        
        // solve feedback
        float y3 = (g*s[3] + f3*g*s[2] + f2*g*s[1] + f1*g*s[0] + f0*input) / (1 + r*f0);
        
        // then solve the remaining outputs (with the non-linear gains here)
        float xx = t0*(input - (r*y3));
        float y0 = g*(s[0] + f*xx);
        float y1 = g*(s[1] + f*y0);
        float y2 = g*(s[2] + f*y1);
        
        // update state
        s[0] += 2*f * (xx - y0);
        s[1] += 2*f * (y0 - y1);
        s[2] += 2*f * (y1 - y2);
        s[3] += 2*f * (y2 - y3);
        
        *inputL = y3; 

Code: Select all

    float processTanh(float x, float x1)
    {
        if (fabs(x - x1) > 0.0000001f)
        {
            return (calcAntiderivate(x) - calcAntiderivate(x1)) / (x - x1);
        }
        else
        {
            return calcAntiderivateSmallValue(x, x1);
        }
    }

    float calcAntiderivate(float x)
    {
        return log(cosh(x));
    }
    
    float calcAntiderivateSmallValue(float x, float z1)
    {
        // edited
        return tanh((x + z1) / 2);
    }

[Richard_Synapse]

I tried the above code but it doesn't seem to work right, I see hardly any change in alias noise?

The method works fine though, just for fun I tried it with an iterative solver as described in the paper, and there was a strong reduction in aliasing, perhaps 20-30 dB or so. What I did was simply to replace the current and past non-linearity of the trapezoidal solver with your processTanh() function.

Richard

[xhy3]

I tried the above code but it doesn't seem to work right, I see hardly any change in alias noise?

The method works fine though, just for fun I tried it with an iterative solver as described in the paper, and there was a strong reduction in aliasing, perhaps 20-30 dB or so. What I did was simply to replace the current and past non-linearity of the trapezoidal solver with your processTanh() function.

Richard

20 - 30dB is a lot, great to hear that it worked.

I fixed the calculation for very small x - x1 in the code above (ill condition). I also came to the conclusion that the alias reduction is very small or even not notable with the cheap non linear filter approach above. But i'm pretty sure it should work too with the right implementation.
Another problem is the computation power it needs, but its also possible to use an antiderivate of a tanh approximation instead of tanh. This saves some CPU...

I will make some more tests and measuring. Maybe i also try an iterative solver.

Edit: Maybe we just have to delay the z1 value instead of taking a delayed filter pole output to make anti aliasing work with the cheap non linear filter approach too.

[Z1202]

The paper says we can fold the antialiased non-linearity into the first pole. I get this when i do that:

Code: Select all

        float y0 = g*(log(cosh(f*(in - r*y3))) - log(cosh(s[0])))/(f*(in - r*y3) - s[0]) * 2;
        float y1 = g*(s[1] + f*y0);
        float y2 = g*(s[2] + f*y1);

I don't have time to look into the details of your code now and I'm a bit out of this subject at the moment, however, if you don't want to use the iterative solver but the linearized ZDF equation instead, you need some appropriate approximation of the antialiased nonlinearity. E.g. a direct analogy to replacing tanh(x) with x would be taking the formula (10) where f(x) is replaced with x, resulting in "restoring back" the half-sample delay we folded into the nonlinearity. You can also obtain this result explicitly by letting f(x)=x instead of f(x)=tanh(x) and then applying the "nonlinearity".

[xhy3]

The paper says we can fold the antialiased non-linearity into the first pole. I get this when i do that:

Code: Select all

        float y0 = g*(log(cosh(f*(in - r*y3))) - log(cosh(s[0])))/(f*(in - r*y3) - s[0]) * 2;
        float y1 = g*(s[1] + f*y0);
        float y2 = g*(s[2] + f*y1);

I don't have time to look into the details of your code now and I'm a bit out of this subject at the moment, however, if you don't want to use the iterative solver but the linearized ZDF equation instead, you need some appropriate approximation of the antialiased nonlinearity. E.g. a direct analogy to replacing tanh(x) with x would be taking the formula (10) where f(x) is replaced with x, resulting in "restoring back" the half-sample delay we folded into the nonlinearity. You can also obtain this result explicitly by letting f(x)=x instead of f(x)=tanh(x) and then applying the "nonlinearity".

Thanks for that. Will have a look at it. I will try to fix the linearized version as a first step, but can't wait to try the iterative thing. Never had to use an iterative solver so far.

[Urs]

Edit: Anyone already used a Jacobian matrix to solve non linear equations?

/w source code:

http://www.kvraudio.com/forum/viewtopic ... 3&t=459955

(there is a sign error, but it doesn't matter - the code works anyway)

[xhy3]

Edit: Anyone already used a Jacobian matrix to solve non linear equations?

/w source code:

http://www.kvraudio.com/forum/viewtopic ... 3&t=459955

(there is a sign error, but it doesn't matter - the code works anyway)

Hi Urs, thanks for the link. I will have a closer look at it.