Sine hard sync

[Z1202]

We deconstruct the waveform into various integrals of Heaviside, convolve each part separately, swap the integrals with the convolutions (which is fine, they are linear) and drop the identity terms... then add it all back together and sample the result.

This is where the catch might lie. Because this decomposition into Heavisides might immediately expose convergence problems with your convolution, so I'd question the equivalence of such decomposition, unless you try to argue smth like generalized Cesaro convergence, but I wouldn't be immediately sure if even that holds - I'm not really having much experience with those depths of Fourier transform theory.

Now.. the way you actually design the low-pass filter is up to you, it's effectively a separate thing from the BLEP-convolution algorithm that just needs some suitable kernel... but if you do choose w(x)*sinc(x) for some window w as your low-pass filter...

Introducing a window is not a "correct" or equivalent transformation of what we're doing. And it seems to me the reasoning here gets quite handwavy. Especially since we're operating outside the usual convergence boundaries of Fourier transform. This was my whole point. Now, obviously, applying the window before the integration is much easier to reason about and to quantify the respective spectral errors, I guess this is what you're getting at. However there is nothing which immediately suggests that applying a window afterwards is wrong, they are both equally or at least comparably wrong. It's just that applying a window (to the residuals) afterwards might be more difficult to quantify.

[mystran]

Now, obviously, applying the window before the integration is much easier to reason about and to quantify the respective spectral errors, I guess this is what you're getting at.

You're missing the point. If you build the BLEP residues by integrating your kernel as-is, you theoretically have just ordinary linear convolution (provided the identity operator with respect to polynomials assumption holds) with said kernel and whether or not windowing is part of the design of said kernel is irrelevant; it's just a low-pass kernel with an extra constraint on the design.

In fact, when we're doing PolyBLEPs, the only thing that changes is that we use a piece-wise polynomial kernel. Everything else is exactly the same and just like in the case with a windowed sinc kernel, the Fourier-transform of the kernel predicts where the aliasing will be. We could probably use minimax designs too, whatever.

If you couple the windowing into the BLEP algorithm itself, you can no longer just treat it as another case of (linear) convolution, rather it becomes a black box. I don't like black boxes I can't analyze if I can help it, so I'll stick to the BLEP algorithm that implements linear convolution.

[mystran]

We deconstruct the waveform into various integrals of Heaviside, convolve each part separately, swap the integrals with the convolutions (which is fine, they are linear) and drop the identity terms... then add it all back together and sample the result.

This is where the catch might lie. Because this decomposition into Heavisides might immediately expose convergence problems with your convolution,

As far as I can see, if you are willing to treat Dirac-delta as a well-behaved signal with respect to convolution (and related operators), then there's no problem. Because treating Dirac as an "ordinary" signal tends to make everything in signal-processing a whole lot easier and because I have yet to encounter a situation where it wouldn't work in practice, I'm happy to take it as an axiom without proof.

[Z1202]

If you couple the windowing into the BLEP algorithm itself, you can no longer just treat it as another case of (linear) convolution, rather it becomes a black box. I don't like black boxes I can't analyze if I can help it, so I'll stick to the BLEP algorithm that implements linear convolution.

But that's exactly what I'm saying: post-integration windowing of residuals is just more difficult to analyse (I guess that's what you mean by "black box"), not that it's "worse" in any obvious way as a means in principle (except for it being possibly more difficult to find good windowing functions).

As far as I can see, if you are willing to treat Dirac-delta as a well-behaved signal with respect to convolution (and related operators), then there's no problem. Because treating Dirac as an "ordinary" signal tends to make everything in signal-processing a whole lot easier and because I have yet to encounter a situation where it wouldn't work in practice, I'm happy to take it as an axiom without proof.

I'm not 100% happy doing that, as direct/inverse Fourier transforms involving Dirac deltas (especially an infinite number of those) often fail to converge in the ordinary sense, especially after integration. So one needs to be somewhat careful doing anything more than the simplest possible operations with those. Summing up an infinite number of Dirac deltas doesn't fall into that simplest possible class in my book. E.g. upon integration (obtaining Heaviside steps) the sum converges only in Cauchy principal sense, if I'm not mistaken. With polynomials you're entering the realm of differentiated Dirac delta's, which are even "worse" than deltas.

[mystran]

With polynomials you're entering the realm of differentiated Dirac delta's, which are even "worse" than deltas.

Actually.. not really. We never need any deltas or it's derivatives explicitly, we just need them to exist and respect H*d'=d and d*x=x where * is convolution. We only need H as a "proper" signal.

You do not sum infinite deltas. Polynomials have finite number of non-zero derivatives. For signals that are not polynomials, if you evaluate a finite number of derivatives, you're effectively taking a Taylor expansion... which is a polynomial which might or might not converge to your signal, but once it's truncated it's a polynomial with finite degree, hence finite non-zero derivatives and hence can be filtered with BLEPs and whatever aliasing you're left with is the truncation error.

[mystran]

Let the signal be at^2+bt+c if t>0. We can rewrite 2aH^3+bH^2+cH where H^n is n-fold convolution.

If we convolve by some filter f, we have (with * for convolution)
f*(2aH^3+bH^2+cH)
= 2a(f*H^3)+b(f*H^2)+c(f*H)
= 2a(H^3+(f*H^3-H^3))+b(H^2+(f*H^2-H^2))+c(H+(f*H-H))
= 2aH^3+bH^2+cH + 2a((f*H-H)*H^2) + b((f*H-H)*H)+c(f*H-H)
= (t>0?at^2+bt+c:0) + 2a((f*H-H)*H^2) + b((f*H-H)*H) + c(f*H-H).

So f*H-H and it's integrals (convolutions by H) give us the BLEP kernels and to get a practical algorithm we just require that (f*H-H)*H^2 is compactly supported (other than that, f can be anything), which is what I was calling "indentity with respect to polynomials up to some degree" earlier.

No multiplications involved (except scalars and as an optimization to evaluate the trivial part), no deltas or their derivatives anywhere, just a bunch of additions and convolutions. It's a linear filter! To then bound the support from the other side, you add an inverse polynomial with the lower bound shifted in time, so they cancel out.

ps. If we take a+b*t+c*t^2+O(t^3) as a Taylor expansion of some random function at some random point, then we can filter the first part and it's the O(t^3) error term that dictates how much aliasing is left. If the series converges fast, we expect this to be small, if it converges poorly, then we need more terms.

[edit: Got the scalar multipliers wrong earlier because of a brain glitch... obviously the first H doesn't count as integration, fixed now.. H, H^2, 2*H^3, 6*H^4, 24*H^5 and so on.. right?]

[Z1202]

I was not sure which exact process you were referring to, therefore I suggested to exercise care in the vicinity of [strikethrough]tardium barrels[/strikethrough] Dirac deltas, which is what you effectively did, by avoiding Fourier transforms of monomials

Anyway, I guess we're sidetracking. The argument was whether post-integration windowing is wrong or just more difficult to analyse

[mystran]

Anyway, I guess we're sidetracking. The argument was whether post-integration is wrong or just more difficult to analyse

Let dK/dt = k be the lowpass kernel, w be the window, then dwK/dt = w'K + wk. So we have an additive term of the integral of the kernel windowed by the derivative of the window. If you plot this for a Hann window and sinc, https://www.desmos.com/calculator/aw9rrdideq we see that the additive term is kinda like a windowed sinc, but not really as there's (at least) a low-frequency component left from the window.

Is it possible to choose a window that somehow cancels out so we avoid this issue? No idea... but I'd say this is just wrong unless we place further constraints on the window or the kernel somehow.

ps. Not sure about proving this, but it appears that this "error" gets smaller as the window grows larger, so perhaps it tends to zero at the limit such that with a "long enough" kernel it's not necessarily a show stopper... but it still seems like an "error" term in any case.

[Z1202]

I need to meditate a bit on what you wrote (maybe you need to give a bit more detail), but my first question is: why are we windowing the Si function (which is what K would normally stand for, I guess, IIUC?) instead of windowing the respective residual?

[mystran]

I need to meditate a bit on what you wrote (maybe you need to give a bit more detail), but my first question is: why are we windowing the Si function (which is what K would normally stand for, I guess, IIUC?) instead of windowing the respective residual?

Does that matter? Oh.. it does.. https://www.desmos.com/calculator/ospvpzbijq

Well.. it's closer now, but still not quite a windowed sinc... and in fact we see that if we do subtract a windowed sinc, that there appears to be a discontinuity in first derivative at zero.

[rigatoni_modular]

Here is corrected version that computes analytic blep residuals.

https://www.desmos.com/calculator/j3jzbjc10m

Bleps 1, 2 and 3 look about right, but I am not 100% sure about oscillating nature of higher bleps.

Some things to notice: if you want to tabulate these, you need to use windowing and the higher the blep, the more smooth windowing it needs. I don't have the exact recipe how to window analytic bleps, but in theory it shouldn't be too hard to have something workable with some experimenting.

Windowing should counteract the infinite growth of higher order blep residuals too.

A couple things!

I decided to actually go ahead and try this approach. I got decent results for the tabulated residuals but had to apply a Blackman window each stage to keep them from growing a lot.

Screenshot 2023-11-08 at 8.22.21 PM.png

I also modified your Desmos link to scale x by 2*pi and got something SUSPICIOUSLY well-behaved: https://www.desmos.com/calculator/22wrdhswis

I'm not sure whether to trust that version of the residuals or not.

I'm a little frustrated about having to window every stage of residual calculation. My ∆h_0 table looks great but my higher order residuals grow a lot faster than they should. This is what happens if I don't window between each step:

Screenshot 2023-11-08 at 8.31.00 PM.png

Also, here's my code for calculating these tables (yep, I'm using doubles for extra precision), with the extra windowing commented out:

Code: Select all

// See https://www.native-instruments.com/fileadmin/ni_media/downloads/pdf/SineSync.pdf
template<int Z, int O, typename T>
struct HighOrderLinearBlep {
    static constexpr int TABLE_SIZE = 2 * Z * O + 1;
    static constexpr int BUF_SIZE = 2 * Z;
    static constexpr int MID_IDX = (TABLE_SIZE - 1) / 2;
    static constexpr double X_SCALE = (double)1 / (double)(O * 2);

    float resid0[TABLE_SIZE];
    float resid1[TABLE_SIZE];
    float resid2[TABLE_SIZE];
    float resid3[TABLE_SIZE];

    T buf[BUF_SIZE] = {};
    int readPos = 0;
    int writePos = Z;

    inline static double blackmanWindow(double piX) {
        return 0.58 + 0.5 * std::cos(2.0 * piX / Z) - 0.08 * std::cos(4.0 * piX / Z);
    }

    inline void populateResiduals() {
        // Calculate 0th order blep, starting by filling table with sinc
        double calcResid[TABLE_SIZE];
        calcResid[MID_IDX] = 1;
        for (int i = 1; i <= MID_IDX; i++) {
            double piX = M_PI * (double)i * X_SCALE;
            double windowedSincVal = blackmanWindow(piX) * std::sin(piX) / piX;

            // Sinc is symmetrical so put this on both sides of 0
            calcResid[MID_IDX + i] = windowedSincVal;
            calcResid[MID_IDX - i] = windowedSincVal;
        }

        // Integrate sinc
        calcResid[MID_IDX] = 0;
        for (int i = 1; i <= MID_IDX; i++) {
            double integral = calcResid[MID_IDX + i - 1] + calcResid[MID_IDX + i] * X_SCALE;
            calcResid[MID_IDX + i] = integral;
            calcResid[MID_IDX - i] = -integral;
        }

        // Normalize, subtract trivial step to get ∆h_0
        double norm = 1.0 / (calcResid[TABLE_SIZE - 1] - calcResid[0]);
        for (int i = 0; i < TABLE_SIZE; i++) {
            calcResid[i] = calcResid[i] * norm + ((i >= MID_IDX) ? -0.5 : 0.5);
            resid0[i] = calcResid[i];
        }

        // Calculate ∆h_1, ∆h_2, ∆h_3 based on ∆h_0
        constexpr double oneThird = (double)1 / (double)3;
        for (int i = 0; i < TABLE_SIZE; i++) {
            double piX = (M_PI * (double)(i - MID_IDX)) * X_SCALE;

            // ∆h_1
            calcResid[i] = piX * calcResid[i] + M_1_PI * std::cos(piX);
            //calcResid[i] *= blackmanWindow(piX);
            resid1[i] = calcResid[i];

            // ∆h_2
            calcResid[i] = 0.5 * (piX * calcResid[i] + M_1_PI * std::sin(piX));
            //calcResid[i] *= blackmanWindow(piX);
            resid2[i] = calcResid[i];

            // ∆h_3
            calcResid[i] = oneThird * (piX * calcResid[i] - M_1_PI * std::cos(piX));
            //calcResid[i] *= blackmanWindow(piX);
            resid3[i] = calcResid[i];

            DEBUG("%f, %f, %f, %f", resid0[i], resid1[i], resid2[i], resid3[i]);
        }
    }

    HighOrderLinearBlep() {
        populateResiduals();
    }
};

[Z1202]

Does that matter? Oh.. it does.. https://www.desmos.com/calculator/ospvpzbijq

Well.. it's closer now, but still not quite a windowed sinc... and in fact we see that if we do subtract a windowed sinc, that there appears to be a discontinuity in first derivative at zero.

Sorry, didn't have time to figure out the details of your desmos example yet. But the discontinuity sounds strange. I'm not sure what you're doing there though yet. Are you differentiating the obtained BLEP and comparing it to sinc?

The function should be the following:
x(t) = [ BLEP(t) - H(t) ] w(t) + H(t)
where BLEP(t)=(1/pi) Si t and H(t) is the odd-symmetric Heaviside. If w(0)=1, I'm not sure how a discontinuity can arise from an infinitely differentiable (at 0) w(t).

Notice that for higher-order BLEPs you need to match asymptotic derivatives at infinity (which is kinda happening for granted for (1/pi) Si t).

As your d(wK)/dt argument goes, I'm not sure it demonstrates anything, unless we assume for granted that wk is our perfect goal (which I'd question and which makes the whole argument circular).

Also integration of windowed sinc should deliver a wrong step height, unless you scale the result to compensate for it, shouldn't it?

[mystran]

As your d(wK)/dt argument goes, I'm not sure it demonstrates anything, unless we assume for granted that wk is our perfect goal (which I'd question and which makes the whole argument circular).

The discontinuity can arrive when we take [B(t)-H(t)]*w(t) and differentiate: d(B-H)w/dt = dBw/dt - dHw/dt = B'w+Bw' - H'w - w'H. I think it's the w'H term.

As far as "perfect goal" I just want a frequency response I can control.

[Z1202]

The discontinuity can arrive when we take [B(t)-H(t)]*w(t) and differentiate: d(B-H)w/dt = dBw/dt - dHw/dt = B'w+Bw' - H'w - w'H. I think it's the w'H term.

Fair point. Except that w'=0 at t=0, so that doesn't explain your observation. Although maybe there would still be a second derivative discontinuity, which is not too nice either. I wonder however, whether this can become an issue with higher-order BLEPs (although this also can be compensated by an appropriate choice of window function).

[Z1202]

Further thinking. Let u=1-w be the complementary window. Then x(t) = BLEP(t) + u(t)H(t).
Effectively, we have a term u(t)sgn t added to BLEP(t) (the 0.5 coefficient between H(t) and sgn t doesn't affect differentiability and can be ignored). Assuming BLEP(t) is infinitely differentiable, we're interested in the differentiability of u(t)sgn(t). If w(t) is Hann window, then u(t)sgn t looks like two glued cosine branches, which obviously has a 2nd derivative discontinuity. This doesn't look too good, although it's also about the amplitude of the 1/f^3 decay arising from it, maybe it's sufficiently small.