polyBLAMP anti-aliasing in C++

Question

How should 4 point polyBLAMP work. It does not really seem to dampen aliasing.

I have an floating point digital audio oscillator that makes a +1.0f jump at the end of every period.

My question is, have I misunderstood how it should work?

I have been using these as reference: http://www.dafx17.eca.ed.ac.uk/papers/DAFx17_paper_100.pdf https://www.researchgate.net/publication/307990687_Rounding_Corners_with_BLAMP

Here is how I thought it should be done: (the code gets run every sample)

float t = phase / period; //normalized position in period
float dt = deltaPhase / period; // normalized phasestep

float polyBLAMP = 0.0f;
if (t < dt) {
    float d = (t)/dt; // fractional delay
    polyBLAMP = (3.0f*powf(d,5.0f)-10.0f*powf(d,4.0f)+40.0f*powf(d,2.0f)-60.0f*d/2.0f+28.0f)/120.0f;
} else if (t < dt*2.0f) {
    float d = (t-dt)/dt;
    polyBLAMP = (-powf(d,5.0f)+5.0f*powf(d,4.0f)-10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)-5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt) {
    float d = (dt-1.0f+t)/dt;
    polyBLAMP = (-3.0f*powf(d,5.0f)+5.0f*powf(d,4.0f)+10.0f*powf(d,3.0f)+10.0f*powf(d,2.0f)+5.0f*d+1.0f)/120.0f;
} else if (t > 1.0f - dt*2.0f) {
    float d = (2.0f*dt-1.0f+t)/dt;
    polyBLAMP = powf(d,5.0f)/120.0f;
}
float output = osc(phase)+polyBLAMP;

Answer 1

First of all, if your phasor oscillator has a step discontinuity (as opposed to a ramp one) at $\phi\bmod 1 = 0$, then you should use the four-point, fourth order polyBLEP residual, not the fifth-order polyBLAMP one.

Second, the easiest way to simplify the piecewise polynomials is to see that the residuals are the successive integrals of the piecewise approximation to the band-limited pulse, or in the terminology of Esqueda, Välimäki & Bilbao [EVB2016], the “third-order B-spline basis function.” This basis function is needlessly complicated as given; it can be expressed as the difference of two overlapping piecewise polynomials:

$$ \frac{1}{6}\,\begin{cases} (2+d)^3 & \text{if}\ {-2} \le d < 0 \\ (2-d)^3 & \text{if}\ \phantom{-}0 \le d < 2 \\ 0 & \text{otherwise} \end{cases} -\ \frac{1}{6}\,\begin{cases} 4(1+d)^3 & \text{if}\ {-1} \le d < 0 \\ 4(1-d)^3 & \text{if}\ \phantom{-}0 \le d < 1 \\ 0 & \text{otherwise} \end{cases} $$

From this, the successive residuals are trivial integrals of polynomials.

Third, the residuals are centered on $d=0$, where $d$ is the offset from the discontinuity to the immediately previous sample point, and the function extends two samples on either side of $0$. In order to apply them to a function of a phasor, the positive portion of the piecewise definition must be shifted to the $\phi=0$ discontinuity while the negative portion must be shifted to the $\phi=1$ discontinuity. Let $0\le\phi<1$ be the phase and $d\phi = f_0/f_s$ the phase increment. The modified polyBLEP in single-sided form is:

$$ \frac{1}{12}\,\begin{cases} \left(2-\frac{\phi}{d\phi}\right)^4 & \text{if}\ 0\le \phi < 2d\phi \\ 0 & \text{otherwise} \end{cases} -\ \frac{1}{12}\,\begin{cases} 4\left(1-\frac{\phi}{d\phi}\right)^4 & \text{if}\ 0\le \phi < d\phi \\ 0 & \text{otherwise} \end{cases} $$

(for the double-sided form, add two further terms under the substitution $\phi := 1 - \phi$.) In Java:

public static double blep(double t, final double dt) {
    double y = 0;
    if (0 <= t && t < 2 * dt) {
        final double x = t / dt;
        double u = 2 - x;
        u *= u;
        u *= u;
        y += u;
        if (t < dt) {
            double v = 1 - x;
            v *= v;
            v *= v;
            y -= 4 * v;
        }
    }
    return y / 12;
}

and this must be added on the right of the discontinuity and subtracted on its left. Hence, the sawtooth is:

public static double saw(double t, double dt) {
    double v = 2 * t - 1;
    v += blep(t, dt);
    v -= blep(1 - t, dt);
    return v;
}

(note that this is a positive-going ramp. For a negative going one, built on $1-2\phi$, subtract and then add the BLEPs.) A square wave or pulse with variable duty cycle $2d\phi < d < 1 - 2d\phi$ can be built from a difference of phase shifted sawtooth waves, or directly by “polyBLEPping” on both discontinuities in the cycle:

public static double sqr(double t, double dt) {
    double v = Math.signum(2 * t - 1);
    v += blep(t, dt);
    v -= blep(1 - t, dt);
    v += blep(0.5 - t, dt);
    v -= blep(t - 0.5, dt);
    return v;
}

If your oscillator is a triangle wave, then you do need to use the polyBLAMP residual. This can be found by integrating the polyBLEP residual, and scaling appropriately to the slope of the ramp it will be applied to. In single-sided form:

$$ -\frac{d\phi}{15}\,\begin{cases} \left(2-\frac{\phi}{d\phi}\right)^5 & \text{if}\ 0 \le \phi < 2\phi \\ 0 & \text{otherwise} \end{cases} +\ \frac{d\phi}{15}\,\begin{cases} 4\left(1-\frac{\phi}{d\phi}\right)^5 & \text{if}\ 0 \le \phi < d\phi \\ 0 & \text{otherwise} \end{cases} $$

(for the double-sided form, add two further terms under the substitution $\phi := 1 - \phi$.) The polyBLAMP residual in Java:

public static double blamp(double t, double dt) {
    double y = 0;
    if (0 <= t && t < 2 * dt) {
        final double x = t / dt;
        double u = 2 - x, u2 = u * u;
        u *= u2 * u2;
        y -= u;
        if (t < dt) {
            double v = 1 - x, v2 = v * v;
            v *= v2 * v2;
            y += 4 * v;
        }
    }
    return y * dt / 15;
}

Then the triangle oscillator is:

public static double tri(double t, double dt) {
    double v = 2 * Math.abs(2 * t - 1) - 1;
    v += blamp(t, dt);
    v += blamp(1 - t, dt);
    t += 0.5;
    t -= Math.floor(t);
    v -= blamp(t, dt);
    v -= blamp(1 - t, dt);
    return v;
}