Husserl tutorial series(2). Resampling: when Average is Better

Ernest

Tutorial 2: Introduction

This is part 2 of a series on optimal programming in gen~ codebox. The first part showed how to make an optimal LFO, with special considerations of optimal programming of branches and functions. This part examines resampling, which is a copy of a 'Facebook note' I wrote some years ago, then Facebook most unhelpfully stopped supporting 'Facebook notes,' an annoyance which programmers frequently encounter when trying to use uinique and special features. I'd have to say the best advice I could give as an engineer is not to do anything which can't be easily repurposed, and I really should take my own advice better. There will be more code examples of resamplers on my blog soon, at https://yofiel.com. Other tutorials in this series:

Designing a good LFO in gen~ Codebox: https://cycling74.com/forums/gen~-codebox-tutorial-oscillators-part-one
Resampling: when Average is Better: https://cycling74.com/forums/gen~-codebox-tutorial-oscillators-part-2
Wavetables and Wavesets: https://cycling74.com/forums/gen~-codebox-tutorial-oscillators-part-3
Anti-Aliasing Oscillators: https://cycling74.com/forums/husserl-tutorial-series-part-4-anti-aliasing-oscillators
Implementing Multiphony in Max: https://cycling74.com/forums/implementing-multiphony-in-max
Envelope Followers, Limiters, and Compressors: https://cycling74.com/forums/husserl-tutorial-series-part-6-envelope-followers-limiting-and-compression
Repeating ADSR Envelope in gen~: https://cycling74.com/forums/husserl-tutorials-part-7-repeating-adsr-envelope-in-gen~
JavaScript: the Oddest Programming Language: https://cycling74.com/forums/husserl-tutorial-series-javascript-part-one
JavaScript for the UI, and JSUI:<a href="https://cycling74.com/forums/husserl-tutorial-9-javascript-for-the-ui-and-jsui"> https://cycling74.com/forums/husserl-tutorial-9-javascript-for-the-ui-and-jsui
Programming pattrstorage with JavaScript: https://cycling74.com/forums/husserl-tutorial-series-programming-pattrstorage-with-javascript
Applying gen to MIDI and real-world cases. https://cycling74.com/forums/husserl-tutorial-series-11-applying-gen-to-midi-and-real-world-cases
Custom Voice Allocation. https://cycling74.com/forums/husserl-tutorial-series-12-custom-voice-allocation

Optimal Resampling: when Average is Better

Before the really good stuff, some oddities in the way almost everyone does the simplest form of resampling, ‘2x linear interpolation,’ are well worth considering. It’s invariably written as adding a fraction of the difference between two samples to the earlier one, to estimate a point half between::

    interp  = z1 +(z0 -z1) /2

I’ve never heard anyone remark how even this simple equation contains an unnecessary subtraction. Some optimize by replacing the divide by an multiply, but stop there. Actually, you don’t need to make the fractional component and add it to the base. All you’re actually doing is averaging the numbers:

     interp  = (a +b) *.5

On modern CPUs, floating-point divides take ~40 clock cycles, multiplies take 4 clocks, and add/subtracts take 2 clocks, assuming the data and instructions are already cached. So this lower form is ~8 times faster. Also, one can also extend the linear technique and get even more optimization. Typically, people extend the first equation above for 4x linear upsampling like this:

  x1  = z0 + (z1 -z0) *.25;

  x2  = z0 + (z1 -z0) *.5;

  x3  = z0 + (z1 -z0) *.75;

  x4  = z1;

That’s the obvious thing to do. So you might think my little trick here only optimizes the middle upsample. But there’s another way to get similar quality with only 2x upsampling! Formally called 4-way piecewise linear interpolation, it can improve the accuracy for the small price of a few adds. It also smooths effectively near Nyquist, handles reverses in wave direction well, and, using my nifty first trick above, optimizes to almost nothing!

  x  = (z0 + z1 +z2 +z3) *.25

Unless you want to add weighting factors, though, you can’t use odd numbers of points, or extend it much beyond four points. And with the optimization below, you can do MUCH better, with for just about the same amount of CPU as 4x-linear upsampling.

The optimal Resampler? The distortion of linear interpolation is ~-50dB. In some cases the noise from resampling doesnt matter, but depending on what you’re resampling, it can be quite audible as white noise, and for some content it can even reach -45dB.

For something better, there’s much argument about the best compromise between performance and quality. “sinc windowed filtering” is definitely the best if you make the window large enough. But larger sinc windows are CPU intensive. The Catmulli Rom Spline, which is a four-point, third-order (‘cubic’) polynomial, is generally better than alternatives of the same order, sometimes by about an order of magnitude (you’d need a 32-point sinc to do better). Here’s a function for spline resampling in MAX .genexpr format:

splineInterp4(in1, x){                  // x = offset between z1 and z2.     History z0, z1, z2, z3;    
     x  = x * x * x * (-0.5*z0 + 1.5*z1 - 1.5*z2 + 0.5*z3)
         + x * x * (z0 - 2.5*z1 + 2*z2 - 0.5*z3)
         + x  * (-0.5*z0 + 0.5*z2)
         + z1;
     z3 = z2; z2 = z1; z1 = z0; z0 = in1;
     return  x;
 }

But it takes 18 multiplies and 11 add/subtracts! Compared to other methods, its low on your cpu, but you can still improve it. Suppose you want 2x upsampling. Then you know the interval between the two middle points is always 0.5, hence, you can ‘solve the equation’ for x=.5 and get this (I cheat using Wolfram Alpha to solve it for me):

splineInterp45(in1){
     History z0, z1, z2, z3;    
     x = - 0.0625 * z0 
         + 0.5625 * z1 
         + 0.5625 * z2 
         - 0.0625 * z3;
     z3 = z2; z2 = z1; z1 = z0; z0 = in1;
     return  x;
 }

So your 29 floating-point ops have become 8, just like that! A side note: solving the equation this way makes it the same as ‘weighting factors’ on a 4-way piecewise linear approximation, but the generalized form for determining points at any arbitrary interval on a curve has more applications, so that’s how it’s usually implemented, even though it’s less efficient if you have regular, fixed intervals, like in resampling.

If you want 4x resampling, you can just do the same thing for .25 and .75, in the same way as for 4x linear oversampling. If you do, you might be astonished to find you can do 4x resampling for only 17 multiplies and 16 add/subtracts, because you can do even more term reduction on the downsampling. You can also stack splines in layers to take in more points and further increase quality.

These optimizations can make a huge difference on, for example, a 32-voice synthesizer with multiple upsamplers in each voice for oscillators and filters. The functions in commercial DAW software are more flexible and contain more range checking, and other options perhaps. But if you are doing alot of resampling, it can easily take up an entire core thread on your CPU, even with the faster computers of 2019.

But I dare say, you won’t find many other DSP engineers out there who would go through the effort. Most just use built-in libraries and have no idea how they work, which is one reason why high-voice-count software synthesizers of any quality are still quite rare. And the experts are too busy arguing Blackman versus Hamming windows for sinc. It takes an average mathematician like me to bother, I guess :)

enrico wiltsch

Thx!