KVR Audio

Amusesmile · Post by **Amusesmile** » Thu Nov 01, 2012 9:24 pm

Alright, I've spent the past day or two researching BLEPS and minBLEPS and various other methods for generating bandlimited waveforms. There's a general consensus that hard sync is awesome and you can only really pull it off with the minBLEP method as described by Eli Brandt in this paper:

http://www.cs.cmu.edu/~eli/papers/icmc01-hardsync.pdf

So I'm going to go with minBLEPs. There's various, somewhat conflicting but universally dizzying example code floating around on this forum as well as musicdsp.org. The most promising that I've come across is Daniel Werner's method (http://www.experimentalscene.com/articles/minbleps.php) which I have running now. Basically does all the difficult work of actually initializing the minBLEP table, which ends up as an array of floats that looks something like this (depending on the number of zero crossings and oversampling rate):

0.0037742
0.00912649
0.0168908
0.0302262
0.0543585
0.0956485
0.1601
0.251603
0.370363
0.512015
0.667792
0.825851
0.973478
1.09959
1.19687
1.26284
1.29977
1.31327
1.31035
1.29735
...

Ok awesome, so I should be able to use this generated minBLEP with my naive sawtooth class to get a nice, bandlimited, hard-syncable sawtooth. That's where I'm stuck. I understand that you're supposed to insert these minBLEPS wherever you have a zero crossing in the naive signal which is easy enough to figure out, but then you also need to keep track of how many are in progress, which also isn't that difficult but I'm a little shaky on how everything is added together and output. I found this in a thread from way back:

"You mix in the minBLEP step any time there is a discontinuity in the signal.

Assuming you have the discontinuity at position 3.4, for example, the procedure goes like this:
minblep startpos = 0.4*oversampling
sample[x] = trivialsaw[x] + minblep[startpos+x*16];
once x > minblep length/oversampling you stop mixing in that minblep."

So by 3.4 I assume he means phase offset 3.4, or 3.4 times through the waveform, which is the same as phase 0.4. Alright that's all well and good. So if I understand this correctly, you just detect zero crossings in the naive signal, and whenever that happens you take the naive signal's current phase and use that offset to reach into the minBLEP table and output the signal value + the corresponding minBLEP value. You keep doing that until the end of the minBLEP table is reached, at which point the BLEP essentially dies. If you're at a high frequency so that many minBLEPs are overlapping, you still just do the same thing with all of them at once (you could have like ~10 or more). Is this correct? That seems too simple or otherwise somehow fundamentally flawed.

Any advice is appreciated and (on my honor), when I get this working I'll post an open source example with sweeping saw, sweeping square, and hard synced versions of both. If this is the best antialiasing solution it shouldn't be so difficult to track down!

AdmiralQuality · Post by **AdmiralQuality** » Thu Nov 01, 2012 11:40 pm

I know, funny how there's still nothing to just "copy" it from.

Still, I'm kind of glad it's a bit hard, it forces us to understand what's happening.

I use another method I designed with some math-help from a well qualified friend. Same idea, just without the lookup table, and the function shape looks different.

My issue is actually in the overlapping them, that's the final bit I'm working on now in a new anti-aliasing implementation. Overlapping them is easy enough, but the combining of them is giving me some grief. In my case the function wave is all positive, so I multiply it by the naive waveform. But multiplying two (or more) on top of each other (by recursion) is NOT sounding nice. (This just gave me some ideas of a few more things to look at. Maybe something needs inversion...)

Anyway, following along... Curious and yes, that would be a good resource.

mystran · Post by **mystran** » Fri Nov 02, 2012 12:06 am

The basic algorithm (for a saw) is to increment oscillator phase, then while (and yes I mean a while-loop) the phase is larger than 1 (or whatever reset point you want) you (1) solve the exact time within the sample that the phase crosses 1 (2) pick a branch from the BLEP (which is essentially a large polyphase FIR) and mix/scale it into a buffer (starting from current sample position). Then substract 1 (or other reset value) from the phase, and repeat (possibly mixing more BLEPs).

Once you're done, add the trivial value into the same buffer, with current sample offset + the latency of the BLEP (for linear phase this is half the length and minBLEPs are a fail). Once the trivial value has been added, you advance to the next sample and repeat the whole thing. Note that this algorithm also solves the whole latency thing trivially, will work gracefully above Nyquist, and can easily be extended to any piece-wise polynomial waves (with some additional CPU cost ofcourse).

Even if you decide to go with the minBLEPs, I would STILL suggest rendering it all into a buffer, because that's the only sensible way to do it.

For more complicated waves I suggest replacing the loop with a while(true) and branching out with a break once you know that you don't need any more BLEPs. For pulse-waves for example, you'd keep track of the state (high or low) and assuming low-state first, you'd compare "low and phase > pw" (solve for crossing PW to get the blep branch) and "high and phase > 1" (reset back to the beginning). Note that for PWM you want to treat the PW as a piece-wise linear (rather than constant during the sample) or you'll run into some artifacts.

Sync is a bit trickier if you want to do it right, but basically you have to find out whether slave has any events with (sub-sample) offsets preceding a master reset, then process those events, then do the reset. Easiest if you combine both oscillators into a single piece of code (which also means you can reuse the same buffer, as long as you do volume-scaling at the BLEP mixing stage).

Amusesmile · Post by **Amusesmile** » Fri Nov 02, 2012 6:38 am

AdmiralQuality wrote: Still, I'm kind of glad it's a bit hard, it forces us to understand what's happening.

Yeah, I feel you. The scavenger hunt is a good portion of the fun. Still, I'm making a project of finding, cleaning and archiving a bunch of the older stuff. It's starting to be hard to find. On that subject, this is probably the most useful link I've found on BLEPs, and you can't even access it online anymore. I was using badly cached google versions until I found it on the wayback machine:

http://web.archive.org/web/201106110423 ... /bl-synth/

Cool. So mystran, thanks so much for the walkthrough. I'm going to try to go through this step by step. I think I almost have it conceptually.

Just so I know I have this straight, in your example we're working with 3 buffers. The first is a naive table saw, the second is a pre-computed BLEP, and the last is a circular mix/sum buffer that we actually read from at the current sample offset to get an ouput value. Into the mix/sum buffer we're always putting in the current naive saw delayed by half the size of the BLEP and adding any BLEPs that happen at zero crossings. Alright if this is correct I'll have it soon (knock on wood). Thanks.

AdmiralQuality · Post by **AdmiralQuality** » Fri Nov 02, 2012 7:54 am

kryptonaut · Post by **kryptonaut** » Fri Nov 02, 2012 12:47 pm

Having used a MinBLEP-like algorithm for creating hard-synced waveforms, here's my take on it.

The MinBLEP table as described is essentially an oversampled table giving the response of a low-pass filter to a 0->1 step input. Typically the response will be a slow-ish rise to somewhere above 1, then a few decaying oscillations about 1.

Basically the idea is to use this precomputed result to simulate what would have happened if the whole naive waveform had been oversampled and sent through the low-pass filter to bandlimit it at (or just below) Nyquist, and then resampled back down to the original rate.

So any time there is a step of size S in the naive waveform it needs to be replaced with S times the pre-computed filter response, downsampled appropriately. Since the step will generally occur 'between' samples, we start the lookup in the table at an appropriate offset. For example, assume the table is 16x oversampled. If the step occurs 0.25 of the way through a 'normal' sample interval then we want to start reading the table at the (0.25x16)th location, stepping by 16 for each subsequent sample. You can interpolate the table for additional accuracy if you like.

Now, how do we use this to actually generate a bandlimited sawtooth? A naive sawtooth might be generated by an accumulator which is decreased by a fixed amount each sample step, and whenever the accumulator falls below -0.5 it is increased by 1. We want to replace that abrupt step with the filtered version taken from our MinBLEP table. We could set a pointer to the start of the table and read it in step with our sample output rate, adding the result to the accumulator (which would now be falling below -0.5), but when we reached the end of the MinBLEP table we'd have some fiddling to do because the table would contain values of around 1 and the accumulator would be around -1.5 - a better approach is to offset all the table entries by -1 so that they represent a filtered unit step from -1 up to 0. Then when we reach the end of the table we can simply stop adding values since they will be 0 anyway. To compensate for this offset we also offset the accumulator in the opposite direction whenever we start reading the MinBLEP table.

So now when the accumulator reaches -0.5 (or less) we add 1 to it as before, but also start adding entries from the (offset by -1) MinBLEP table. Those values will initially be around -1, so the net result is that the accumulator only rises slowly towards +0.5, oscillates a bit (due to the ripples in the MinBLEP table), and then by the time we reach the end of the MinBLEP the accumulator is back to where the naive version would have been.

However we also have to consider the fact that another step may arrive before the response from the previous one has decayed away, and also that there may be more than one step per sample (e.g. if we are simulating hard-sync). The way to do this is to maintain a circular buffer whose length is at least the number of samples it takes for the MinBLEP table to settle to a constant value (or near enough), divided by the oversampling rate of the table. Maintain a pointer into this circular buffer which increments in step with the sample output buffer. When a step occurs in the naive waveform we read a complete subsampled (and offset in time by the appropriate fraction of the supersampling rate) copy of the MinBLEP table and add it into the circular buffer starting at the current position. The values to be added also need to be scaled by the amplitude of the step in the naive wave. Then in addition to taking the naive step, we add the content of the current entry in the circular buffer (which may now include contributions from previous MinBLEPs). The entry in the circular buffer can then be zeroed ready for the next time it is used.

It's a complicated description, but to my mind the best way to think of it is just as a pre-computed low-pass filter response, which only gets added at abrupt steps in the waveform. Using an oscilloscope to visualise the result is invaluable when developing this kind of algorithm as it can be tricky to ensure all the replacement steps have the right amplitude and sign - but when it finally works it works well! Getting this going is partly what motivated me to produce J-Scope: http://jaggedplanet.com/VST/jscope.html

mystran · Post by **mystran** » Fri Nov 02, 2012 12:57 pm

Amusesmile wrote: Just so I know I have this straight, in your example we're working with 3 buffers. The first is a naive table saw, the second is a pre-computed BLEP, and the last is a circular mix/sum buffer that we actually read from at the current sample offset to get an ouput value. Into the mix/sum buffer we're always putting in the current naive saw delayed by half the size of the BLEP and adding any BLEPs that happen at zero crossings. Alright if this is correct I'll have it soon (knock on wood). Thanks.

When I say "naive saw" (or other "naive waves") I simply mean analytically calculated values for the waves without any anti-aliasing considerations. In case of saw, this is simply the phase-value (or 2*phase - 1 for bipolar). For pulse-waves the "naive wave" would be something like (state == low ? -1 : 1). You don't need any buffer for this.

So for saw/pulse you need 1 buffer (the mix buf) and 1 read-only table (the BLEP data).

Amusesmile · Post by **Amusesmile** » Fri Nov 02, 2012 8:55 pm

Blarg. Alright I've got something resembling a bandlimited saw at low frequencies but it's still aliasing like hell as you go higher. Thanks so much for both of those descriptions. Yeah mystram sorry- I was computing a saw shape into a table and reading from that just like I do with my sine oscillator class. The phase method is more efficient for sure. Thanks.

So now I'm creating a unipolar naive wave from 0.0 - 1.0 and while the phase is above 1.0, we subtract 1.0 from the phase, compute the exact sub-sample moment at which we crossed 1.0 (exactCrossTime = (phaseIncrement-phase)/phaseIncrement), and then insert our BLEP into the circular buffer offset by the exactCrossTime.

Here's my update loop right now:

Code: Select all

double update()
{
    phase += phaseIncrement;
    index = (index+1)%4;
    while(phase >= 1.0)
    {
        phase -= 1.0;
        double exactCrossTime = (phaseIncrement-phase)/phaseIncrement;
        for(int i = 0; i < 3; i++)
        {
            double tempIndex = (exactCrossTime*32.0)+(i*32.0);
            double tempFraction = tempIndex-floor(tempIndex);
            circularBuffer[(index+i)%4] += (1.0-LERP(tempFraction, BLEP[(int)floor(tempIndex)], BLEP[(int)ceil(tempIndex)]));
        }
    }
    
    circularBuffer[index] += phase;
    double output = circularBuffer[index];
    circularBuffer[index] = 0.0;
    return output;
}

My BLEP[] table is 128 long, with 32 times oversampling, so my circular buffer is only 4 samples long (128/32). Not sure what's wrong exactly because conceptually I think I'm doing what both of you walked me through. I don't think this is the issue but just so we're on the same page, the BLEP[] looks like this:

Code: Select all

Amusesmile · Post by **Amusesmile** » Fri Nov 02, 2012 9:01 pm

Or here's a simpler version of the code without the linear interpolation on the BLEP samples:

Code: Select all

double update() 
{ 
    phase += phaseIncrement; 
    index = (index+1)%4; 
    while(phase >= 1.0) 
    { 
        phase -= 1.0; 
        double exactCrossTime = (phaseIncrement-phase)/phaseIncrement; 
        for(int i = 0; i < 3; i++) 
        { 
            int tempIndex = (int)(exactCrossTime*32.0)+(i*32.0); 
            circularBuffer[(index+i)%4] += (1.0-BLEP[tempIndex]); 
        } 
    } 
    
    circularBuffer[index] += phase; 
    double output = circularBuffer[index]; 
    circularBuffer[index] = 0.0; 
    return output; 
}

kryptonaut · Post by **kryptonaut** » Fri Nov 02, 2012 10:06 pm

Well I can't see anything obviously wrong with the code. Does the aliasing gradually appear as the frequency rises, or does it suddenly start? How bad is it? One thing that I wonder is whether 4 samples is enough to provide good bandlimiting, but I don't know if that's the issue. Might be worth increasing it as an experiment, though - looking at the numbers in the BLEP table I'm not convinced that it's finished oscillating at the end of the table.

I don't think it's the problem here but if phase ever did get bigger than 2, you'd end up with a negative value for exactCrossTime, which might cause trouble.

I think the best option is to pipe the output into an oscilloscope and a spectrum analyser to get a clearer picture of what's going on as the frequency is gradually increased.

Amusesmile · Post by **Amusesmile** » Fri Nov 02, 2012 10:42 pm

Aight, well I'm glad to hear that there are no obvious errors. While still frustrating it's much less disheartening to know that the concept is right.

I built a scope in the environment I'm testing in. It's obvious that there's a DC offset introduced especially now that I've increased the size of the BLEP, taking your suggestion. Also the waveform gets smaller and smaller as the frequency increases. So that's probably most of the problem right there. Any idea what's causing the offset? Some of the earlier threads talk about this so I'm going to try to figure it out there. Here are photos of the scope output:

http://amusesmile.com/temp/saw1.png
http://amusesmile.com/temp/saw2.png
http://amusesmile.com/temp/saw3.png

Amusesmile · Post by **Amusesmile** » Fri Nov 02, 2012 11:11 pm

Wow, ok so subtracting the phaseIncrement at every step seems to have completely fixed the DC offset, however the aliasing is still horrible. This kind of raises more questions than it answers because it means that the DC offset is being introduced by the phaseIncrement in the first place somehow. Or else it's intrinsically linked. Arg:

http://amusesmile.com/temp/saw4.png

Code: Select all

    
circularBuffer[index] += (phase-phaseIncrement);
double output = 1.0*circularBuffer[index];
circularBuffer[index] = 0.0;
return output;

kryptonaut · Post by **kryptonaut** » Fri Nov 02, 2012 11:40 pm

The offset is in part because your raw wave goes between 0 and 1 rather than -0.5 to +0.5 But the BLEP introduces a further DC offset which increases with frequency because it effectively slows down and delays the downward step, so (in this case) the waveform spends less time in the negative part of the cycle than the naive wave would, so ends up with a net positive bias.

As to why the aliasing problem is persisting - I'm surprised not to see any ringing effects at the steps. Are you sure you have the BLEP oversampling ratio correctly specified in the calculation of tempIndex? I would expect the values that get put into the circular buffer to start near 1, fall towards zero for a sample or two, then oscillate around zero in a decaying manner. When these are added to the naive phase I'd expect to see a bit of overshoot then some ringing before settling back to the straight naive saw, especially at low to mid frequencies. If you can get in with a debugger, check the values being stored in the circular buffer to see what they look like.

Amusesmile · Post by **Amusesmile** » Sat Nov 03, 2012 12:23 am

kryptonaut wrote:The offset is in part because your raw wave goes between 0 and 1 rather than -0.5 to +0.5 But the BLEP introduces a further DC offset which increases with frequency because it effectively slows down and delays the downward step, so (in this case) the waveform spends less time in the negative part of the cycle than the naive wave would, so ends up with a net positive bias.

No I totally understand the 0.0 to 1.0 offset. We can subtrack 0.5 from the time to go from -0.5 to 0.5 instead. No worries there.

The extra DC offset introduced makes sense because the wave is taking longer to snap down like you're saying. If the aliasing is fixed I can empirically find an amount to subtract based on the phaseIncrement and the oversampling amount and zero crossings.

So the aliasing is a big problem. I tried outing the values every time we add a blep. Won't the ringing only last ~8 samples (depending on the blep buffer size) no matter what frequency the wave is at though? Or is that a fundamental problem with how I'm approaching this? The BLEP always lasts the same amount of time right?:

add blep with 4 zero crossings and 16 oversample
tempIndex: 13.7233 blep add: 0.962785 circularBuffer[3]: 0.962785
tempIndex: 29.7233 blep add: 0.417643 circularBuffer[4]: 0.417643
tempIndex: 45.7233 blep add: -0.124307 circularBuffer[5]: -0.124307
tempIndex: 61.7233 blep add: 0.0248652 circularBuffer[6]: 0.0248652
tempIndex: 77.7233 blep add: -0.0525997 circularBuffer[7]: -0.0525997
tempIndex: 93.7233 blep add: -0.00420365 circularBuffer[0]: -0.00420365
tempIndex: 109.723 blep add: -0.0138244 circularBuffer[1]: -0.0138244
tempIndex: 125.723 blep add: -0.00579836 circularBuffer[2]: -0.00579836

add blep with 4 zero crossings and 16 oversample
tempIndex: 15.1778 blep add: 0.931389 circularBuffer[3]: 0.931389
tempIndex: 31.1778 blep add: 0.273993 circularBuffer[4]: 0.273993
tempIndex: 47.1778 blep add: -0.105348 circularBuffer[5]: -0.105348
tempIndex: 63.1778 blep add: 0.0157142 circularBuffer[6]: 0.0157142
tempIndex: 79.1778 blep add: -0.050905 circularBuffer[7]: -0.050905
tempIndex: 95.1778 blep add: 0.0018565 circularBuffer[0]: 0.0018565
tempIndex: 111.178 blep add: -0.0187399 circularBuffer[1]: -0.0187399
tempIndex: 127.178 blep add: -0.000365654 circularBuffer[2]: -0.000365654

add blep with 4 zero crossings and 16 oversample
tempIndex: 0.632341 blep add: 1.00019 circularBuffer[4]: 1.00019
tempIndex: 16.6323 blep add: 0.922327 circularBuffer[5]: 0.922327
tempIndex: 32.6323 blep add: 0.243952 circularBuffer[6]: 0.243952
tempIndex: 48.6323 blep add: -0.0988503 circularBuffer[7]: -0.0988503
tempIndex: 64.6323 blep add: 0.0126158 circularBuffer[0]: 0.0126158
tempIndex: 80.6323 blep add: -0.0499097 circularBuffer[1]: -0.0499097
tempIndex: 96.6323 blep add: 0.00271958 circularBuffer[2]: 0.00271958
tempIndex: 112.632 blep add: -0.0194804 circularBuffer[3]: -0.0194804

mystran · Post by **mystran** » Sat Nov 03, 2012 1:27 pm

The BLEP data you posted looks wrong. The BLEP data should look like a filtered step, with the step removed. What it looks like is more like a low-pass "bump" without the integration or the step-removal.

Secondly, it doesn't look like "minimum phase" (which I would argue you don't want anyway) so you should delay the trivial wave by half the "BLEP" length. What this means is that you should add the trivial wave with an offset.

Thirdly, forget linear interpolation, and just use enough branches. Cheaper. I find "enough" branches to be somewhere around 4k or so. As for the length of a single branch, I'd test with 16 or 32 points or something, because at that point you can still expect "clean" results.. so you know any aliasing is a bug, and not just a limitation of a cheap kernel.

Fourthly: just make your ringbuffer a few times larger. You can try to minimize the size once the basic stuff works. It doesn't hurt to have some extra buffer space while debugging though, so you don't need to worry about 1 sample offsets.

minBLEPS once and for all