Tutorial: BLEPs (using PolyBLEPs but extensible)

[mystran]

AFAIK it's a new method, with good aliasing cancellation and (since it's min-phase) no lookahead. The example implementation supports up to 3rd-order BLEPs (where 1st is a step discontinuity, 2nd is a ramp).

How do you solve the issue where min-phase normally leads to frequency dependent DC offsets due to time-alignment mismatch between the naive and the BLEP?

[signalsmith]

I'm not sure what you mean, actually.

My understanding of BLEP methods in general is that you add in an anti-aliasing signal which is the difference between the continuous-time "feature" (step/ramp/whatever) and the (approximately) lowpassed equivalent, such that sampling/calculating it at discrete points isn't as bad. This difference should have no DC component in a sensibly-designed setup, so I'm not sure how adding it in would produce an offset.

(This MinBLEP description seems to run into DC-offset trouble because of (a) truncating the lowpass kernel, and (b) integrating using a cumulative sum which produces a half-sample time-shift. Both of those seem possible to avoid, but they're also not relevant when analytically integrating a continuous-time IIR filter.)

[mystran]

My understanding of BLEP methods in general is that you add in an anti-aliasing signal which is the difference between the continuous-time "feature" (step/ramp/whatever) and the (approximately) lowpassed equivalent, such that sampling/calculating it at discrete points isn't as bad. This difference should have no DC component in a sensibly-designed setup, so I'm not sure how adding it in would produce an offset.

Let me rephrase: do you take the phase-mismatch into account somehow when differencing the signals?

With ordinary tabulated BLEPs, that's the real problem with minimum-phase: the "scale by trivial step" is only correct if both the naive and the filtered signal are in phase and the only obvious way to force this is to make the filter into a linear-phase one and then delay the naive until it's in alignment.

But .. it's entirely possible I'm just missing some nuance with the analytical IIRs.. I've looked into the theory, but not spent much time with it.

[signalsmith]

[...] the "scale by trivial step" is only correct if both the naive and the filtered signal are in phase [...]

I think the naive feature you're subtracting has to be in phase with the input, but the band-limited one you're replacing it with doesn't have to be - for signals made entirely from polynomial segments, anyway. (And the band-limited features of different order need the same phase response.)

When it's a mixture of BLEP-able discontinuities and other stuff (e.g. hard-sync sine) then you do get phase issues. It doesn't result in aliasing, but things can add up wrong in the passband. Is that what you're talking about? They'll always match phase at DC anyway, so the difference wouldn't have low frequencies.

[mystran]

When it's a mixture of BLEP-able discontinuities and other stuff (e.g. hard-sync sine) then you do get phase issues. It doesn't result in aliasing, but things can add up wrong in the passband. Is that what you're talking about? They'll always match phase at DC anyway, so the difference wouldn't have low frequencies.

So you match the naive against the minimum-phase kernel's latency at DC?

[signalsmith]

No, but I still don't see why it would be necessary.

Here's a 0-24kHz sawtooth sweep using the example oscillator from my GitHub repo, along with a simple 100Hz lowpass of the same waveform, showing that there's no DC issue:

sine-sweep.png

I think I've found the issue you're referencing, but I don't think it's a theoretical one. It really seems like an off-by-one error or bad windowing or similar, not an inherent part of MinBLEP or any other BLEP schemes.

[Gordonjcp]

crazy tricks with pointers and precalculated lookup tables.

Hello Gordon, would you mind sharing how did you implement a lookup table for this algorithm? The most successful approach for me was to compute a 2D LUT for t & dT (which is both inefficient and not quite accurate - I get about 78dB noise floor compared to 98db without LUT). The table is also not efficiently packed. The values are all stored in the upper triangle, and it requires me to use some non-trivial mapping to not waste space on a microcontroller's flash.

I'll do you one better, almost 18 months after I posted, I'll just link to the code:

https://github.com/ErroneousBosh/slttblep

You'll notice that the table is computed for a phase angle from 0 to 255 fractions of a cycle, and since the AVR cannot divide but can kind of multiply the division is just multiplication by the reciprocal, again scaled from 0 to 255.

[9b0]

I'm trying to build a hard-synced polyblep oscillator. The amount of resources for this is pretty limited online (most of them are on this forum), and I feel like, I'm pretty close to my goal, but there's a final thing I cannot solve.

I have an issue with overlapping BLEPS, and I can't get my head around a solution. The issue is very disturbing at higher master oscillator frequencies and is less audible on lower ones. I've recorded a small video that showcases the issue:


It can be seen on the scope called BLEP, that the effect appears when two BLEP's overlap (the blue and orange ones).

I'm programming the oscillator in LUA (I'm not a coder anyway), since Alpha Forever has a LuaJIT compiler, and this allows me for quick prototyping and measuring. I'm calculating two BLEP's. The size of the BLEPs is 4 samples (that's why I have to mix them with 2 samples delay).

Maybe I should also scale the second BLEP?

Code: Select all

    local F={F1,F2}
    local p2=phase[2]
    for i=1,2 do
        inc[i]=F[i]*sRR -- calculate the incremental
        inc[i]=min(inc[i],0.25) -- limit the incremental
        phase[i]=phase[i]+inc[i] -- update the phase
        flip[i]=trunc(phase[i]) -- if phase>=1 then flip=1
    end
    if phase[1]>1 then
        phase[1]=phase[1]-1 -- reset phase 1
        d[1]=phase[1]/inc[1] -- calculate the intersample position of the phase crossing 1
        phase[2]=d[1]*inc[2] -- reset phase 2 with respect to phase 1
        scale=p2-phase[2]+inc[2] -- calculate the scaling factor for the blep based on the new value of phase 2
        polyBlep(blep[1],d[1],blepIndex,scale) -- calculate the blep
    elseif phase[2]>1 then
        phase[2]=phase[2]-1 -- reset phase 2
        d[2]=phase[2]/inc[2] -- calculate the intersample position of the phase crossing 1
        polyBlep(blep[2],d[2],blepIndex,1) -- calculate the blep
    end
    y=z[2]-blep[1][blepIndex]-blep[2][blepIndex] -- calculate the output
    for i=1,2 do
        blep[i][blepIndex]=0 -- reset the blep
    end
    z[2]=z[1] -- sample delay
    z[1]=phase[2] -- another delay
    blepIndex=(blepIndex%4)+1 -- increment the blep index
    return y*2-1

[mystran]

Maybe I should also scale the second BLEP?

Every BLEP needs to be scaled, except when the scaling factor would "happen to be" unity anyway, in which case obviously the scaling would do nothing. Getting the scaling right is very important.

The way to compute the scaling factor is that whenever you have a discontinuity, you compute the "limit" on both sides: take the waveform value at the (solved) time-offset where the jump happens before the actual jump has happened and then take the waveform value at the same point after the jump has happened. In terms of the continuous waveform, these are basically the left and right side limits to the discontinuity. Then you compute (after-before) and you have the step magnitude.

Now, when you're doing a natural reset, you're jumping from the very end of the waveform to the very beginning of the waveform, so both of these values are always fixed and you can hardcode them into the algorithm. Further, if your saw goes from -0.5 to 0.5 or 0 to 1 or some such, then the step magnitude is just -1 so "scaling" is just a sign change with respect to regular heaviside (and you can arrange the BLEPs to be inverted by default to skip scaling, as long as you follow this consistently everywhere else too).

When your doing slave reset at arbitrary point, when you solve the slave phase at the master reset point (but just before the actual reset is done; after the reset the phase would be zero again), you need to solve the slave waveform value at that solved phase. The target phase value (assuming we reset to zero phase) is still known, but the BLEP magnitude is the difference.

There's a few other "gotchas" with hardsync: it is possible for the slave to reset just before the master during the same sample in which case you need to solve the slave reset first, then it advances by some fraction of a sample and gets reset by the master which needs another BLEP. It's also possible that the slave gets reset by the master just before it would have reset naturally, in which case you need to only process the master reset and not the slave reset afterwards.

The easy way to figure out which resets to do, assuming that neither the master or slave frequency is so high as to reset several times per sample (this requirement can be relaxed too, but the code logic gets quite complicated), is to simply solve for the master reset first, but not actually process this reset for the slave yet, rather just store the time-offset into a variable. Then you check the slave's natural reset normally and if there's a natural reset, you compare it to the master reset time (if any) and process the slave reset only if it happens before the master reset. Then when you're done with the slave's own processing, whether or not it reset, if there was a sync reset you process it for the slave.

[9b0]

Thanks for the help Mystran, I'll try to work myself through your guides, and make use of them.

Something I'm not sure I'm understanding (linguistically):

take the waveform value at the (solved)

What do you mean by solved?

[mystran]

Thanks for the help Mystran, I'll try to work myself through your guides, and make use of them.

Something I'm not sure I'm understanding (linguistically):

take the waveform value at the (solved)

What do you mean by solved?

BLEP oscillators run in effectively in continuous time. We normally increment ("integrate") the phase for the whole sampling period, but then when we look at the phase at the end of said sampling period and it's above the reset threshold, we know that the continuous-oscillator "should have" reset somewhere in the middle of said sampling period.

So what we do is solve for the exact time-offset where it crossed the reset threshold. Under the assumption that frequency is constant over a single sampling intervals (which is the standard choice to make), the phase (the integral of frequency) is a linear function, which means we can solve the time-offset by solving a linear equation (which is easy enough). The actual reset does not happen "at the sample" but rather at the sub-sample time offset (in the "past" with respect to the current sampling instant we're trying to compute) that we get by solving said linear equation.

This time-offset is what we use to offset the BLEP in time, but when performing hard-sync it's also the time-offset where we need to compute the waveform value of the slave, so that we know exactly how large the jump back to the beginning is, so that we can scale the BLEP magnitude. The slave does not reset "at the sample" but rather at some previous time somewhere between the two sampling instants which we obtained as a solution to the linear equation where the master phase is exactly equal to the reset threshold.

Think of BLEPs as continuous waveforms, that we update optimistically and then look at the phase and then if we see that it's large enough that "something should have happened" we solve exactly where said "something" should have happened and then patch backwards.

[9b0]

Thanks! I understand it now, and this is what I do anyways. I've started this whole project a bit naively. I got polyBlep oscillators working anyways already for a long time, but I was somehow never ineterested in hard-sync. Lately I had an idea I wanted to realize, but it needed sync. I thought, I'd put it together quickly, since I had the idea of running two phase accumulators in parallel, and to apply the Bleps only to the slave phase to obtain a hard-synced saw. I run into lots of issues, and had to figure out how to do the phase reset. Weeks went by since. My current algorithm kind of works, I think, my current issue is going to be this one:

"When your doing slave reset at arbitrary point, when you solve the slave phase at the master reset point (but just before the actual reset is done; after the reset the phase would be zero again), you need to solve the slave waveform value at that solved phase. The target phase value (assuming we reset to zero phase) is still known, but the BLEP magnitude is the difference."

It's strange how underdocumented polyBleps are anyways. The latest paper I could find on hard-sync was on minBleps.

[9b0]

Thanks a lot for Your (Mystran) and Signalsmith's help (he helped on Reddit), I was able to get it done now. I think, this is the most I can ask for from a 4-point polyblep, but it should not be a problem to extend the same algorithm to 8, or more points. I think I've seen 8-point residuals in a paper.

Anyway, here's my result:

On very high master oscillator frequencies it has aliasing, but it sounds really clean on lower notes.

Calculating the slave first, and resetting the master without a blep if the slave and the master flipped in the same sample solved everything. The rest was a bug in my BLEP ringbuffer, that was a leftover of my branchless method from earlier.

Code: Select all

local F={F1,F2}
local p2=phase[2]
for i=1,2 do
    inc[i]=F[i]*sRR -- calculate the incremental
    phase[i]=phase[i]+inc[i] -- update the phase
    flip[i]=trunc(phase[i]) -- if phase>=1 then flip=1
end
if flip[2]>0 then
    phase[2]=phase[2]-1 -- reset phase 2
    d[2]=phase[2]/inc[2] -- calculate the intersample position of the phase crossing 1
    polyBlep(blep[2],d[2],blepIndex,1) -- calculate the blep
    if flip[1]>0 then -- if the master resets during the reset of the slave
        phase[1]=phase[1]-1 -- set the master phase without calculating a new blep
        flip[1]=0 -- set the flip to 0, so we do not care about reseting the master in this cycle anymore
    end
end
if flip[1]>0 then
    phase[1]=phase[1]-1 -- reset phase 1
    d[1]=phase[1]/inc[1] -- calculate the intersample position of the phase crossing 1
    phase[2]=d[1]*inc[2] -- reset phasse 2 based on the new value of phase 1
    scale=p2-phase[2]+inc[2] -- calculate the scaling factor for the blep based on the new value of phase 2
    polyBlep(blep[1],d[1],blepIndex,scale) -- calculate the blep
end
y=z[2]-blep[1][blepIndex]-blep[2][blepIndex] -- calculate the output (subtract both bleps from the naive saw delayed by 2 samples)
for i=1,2 do
    blep[i][blepIndex]=0 -- reset the blep that we already used
    flip[i]=0 -- set flips to 0 (who knows)
end
z[2]=z[1] -- sample delay
z[1]=phase[2] -- another delay
blepIndex=(blepIndex%4)+1 -- increment the blep index
return y*2-1

Thanks for your help again!

[mystran]

Thanks a lot for Your (Mystran) and Signalsmith's help (he helped on Reddit), I was able to get it done now. I think, this is the most I can ask for from a 4-point polyblep, but it should not be a problem to extend the same algorithm to 8, or more points. I think I've seen 8-point residuals in a paper.

Few notes about kernels:

First, the computation cost of polynomial kernels (whether BLEPs or normal interpolators) grows sort of quickly with the order (especially as you'd rather not have them in monomial form for numerical reasons), so beyond some point it's typically faster to just use a precomputed table that you (usually) interpolate linearly.

Second, as the BLEPs get longer, you probably want to consider output into some sort of buffer, rather than shuffling individual delays every sample. I usually tend to use a strategy where I have an output buffer that's long enough to fit a maximum external block size plus a BLEP and then just copy "overflow" from the end of the buffer back to the beginning (while zeroing rest of it) for every block... but a ring-buffer would work too.

Third, as the kernels get longer, you absolutely want to use a strategy where whenever you figure out what BLEP to insert, you insert the entire BLEP into the output buffer "on the spot" rather than tracking what BLEPs might be currently active. I'm not entirely sure if the code you posted already does this, but it's basically the only sane way to allow for arbitrary overlaps (which is always important, but becomes more and more obvious with longer kernels as the probability that they need to overlap grows).

[9b0]

Thanks! I definitely will try this, when I'll give longer Bleps a try. The workflow you write down sounds like convolution, but with an IR that is dependent on the size of the phase resets.