Alias-suppressed waveforms in hardware: DPW, PolyBLEP, other ideas?

[jrmoserbaltimore]

but one nice property of trapezoidal is A-stability which means that it's ROC contains the entire Laplace negative half-plane, so as long as the original system is stable, so is our trapezoidal integration too (at least in exact arithmetics).

This immediately makes sense. Now my problem is not knowing how to place the poles and zeroes (the mathematics of designing a filter are beyond my education level, give me 2 years).

Now.. historically this filter is very influential and arguably quite nice if your alternatives are things like biquads... but with modern understanding that actually doing proper trapezoidal integration of an SVF is essentially just as cheap, there's very few good reasons to bother with Chamberlin anymore.

So, solving problem X in exchange for not as bad problem Y, but we found another way to solve problem X without the downsides.

Not so quick. The details matter a ton. I still haven't properly read the paper, but it looks like it might be fine and essentially just a reordering of computations(?) when compared to the usual trapezoidal SVF.. and really the order of computations is NOT meaningful except as far as different methods of solving Ax=b take different number of clock-cycles.

Under contrasting approaches on page 8:

"While it is true that we placed the three outputs in phase alignment, and as such any delays between them have been removed, it is impossible to completely remove feedback delays, and have instantaneous signals everywhere in the filter. Instead, what we can do is move the 1-sample delays around, which gives the equations better resilience to errors."

Sounds pretty close to reordering the computations. The filter itself involves additional feed-forward paths and so forth, but that affects the transfer function and moves the poles and zeroes around…in ways I don't understand. It might just reorder the computations by moving the delays around.

With direct forms you essentially ALWAYS want to use series or parallel biquad decomposition, because the numerical behaviour gets truly horrid as the order grows. In fact, even in the 2nd order case the performance at low frequencies (where poles are generally pretty close to unit circle) tends to be wonky, to the point where if makes sense to use double precision in a digital EQ if you insist on direct forms.. so good luck trying to keep your bitwidths reasonable with these.

So avoid direct form.

The Chamberlin SVF and the LT SVF are supposed to be good with modifying the frequency, and one of my designs I wanted to attach the cutoff to a low-frequency oscillator. This resulted in me spending a lot of effort trying to find a good way to approximate cos(); I eventually settled on a 6th order polynomial approximation, meaning 3 multiplies, 3 adds. (My original approach was CORDIC, which used only addition and bit shifting to calculate an exact value, plus a final multiplication, but was slower because it calculates one bit at a time) The LT SVF requires calculating tangent, which is hard (division); Chamberlin didn't. I didn't know you could use trapezoidal integration though…I should consider that next time.

For non-modulated filters with a fixed frequency cutoff I'll look into the ladder-lattice topology. Best I've seen so far.

[mystran]

but one nice property of trapezoidal is A-stability which means that it's ROC contains the entire Laplace negative half-plane, so as long as the original system is stable, so is our trapezoidal integration too (at least in exact arithmetics).

This immediately makes sense. Now my problem is not knowing how to place the poles and zeroes (the mathematics of designing a filter are beyond my education level, give me 2 years).

You don't need an "education level" for this stuff. You just need to learn how to do it. In fact the true goal of "education" really is to teach you how to learn stuff, but it is not really a prerequisite, learning stuff is not a licensed profession.

Now... for standard types like Butterworth/Chebychev just grab the formulas from Wikipedia. It's handy to understand that if we want roots at p and q, we can write this as (z-p)*(z-q) or (1-p*z^-1)*(1-q*z^-1) which allows you to readily transform a bunch of poles and zeroes (that is roots of denominator and numerator respectively) into a rational.

[mystran]

Rant mode: Education can be a bit silly at times. You might learn in highschool (or whatever it's called in any given part of the world) that (a^b)*(a^c) = a^(b+c) when a,b,c are natural numbers.. and then perhaps a few years later that (a^b)*(a^c) = a^(b+c) when a,b,c are integers .. and if you're really lucky they might even tell you that a,b,c can be real numbers.. and then perhaps in college you learn that in fact (a^b)*(a^c) = a^(b+c) when a,b,c are complex numbers. And every time they show you complicated proofs that mostly just distract one from the fact that (a^b)*(a^c)=a^(b+c) for most of the types of numbers we care about most of the time.

At some point you'll probably also learn that exp(i*w)=cos(w)+i*sin(w) .. but what they might forget to tell you that now you can use exp(i*b)*exp(i*c) = exp(i*(b+c)) = cos(b+c)+i*sin(b+c) to do trigonometry about as easily as adding numbers together... and that's why I picked up this particular example, because this happens to be quite handy when doing DSP. I think really the most important part of "learning math" (at least for an engineer) is to learn how all the different simple things are related, so that if you can't solve a problem using one method, you can transform it to something else that you can solve... and often when you don't know how to solve something, it's not because you need to learn a lot of math, more often you just need to find the missing piece that relates it to what you already know.

Hint: most of what you learned about polynomials and rationals in high-school works exactly the same when you throw in a few complex numbers (perhaps even better, 'cos now there's slightly less special cases).

[Z1202]

Under contrasting approaches on page 8:

"While it is true that we placed the three outputs in phase alignment, and as such any delays between them have been removed, it is impossible to completely remove feedback delays, and have instantaneous signals everywhere in the filter.

I personally don't appreciate this phrasing. While I understand what the authors are aiming to say, this particular choice of wording can be and often is enormously misleading. It is possible to have instantaneous signals and there are even some variations in how those can be interpreted (e.g. Härmä suggests to understand it in the iterative sense, but it's also possible to understand it in some kind of bandlimited continuous sense, which allows convergence under a wider, actually all practically relevant, range of conditions). In my personal opinion, the majority of attempts to improve Chamberlin's SVF come from the reluctance to simply deal with instantaneous paths. OTOH allowing instantaneous paths gives simple and straightforward answers.

As for the necessity to compute the tangent, I think it's unavoidable if you want to have a BLT response, no matter what the method details are.

[jrmoserbaltimore]

Now... for standard types like Butterworth/Chebychev just grab the formulas from Wikipedia. It's handy to understand that if we want roots at p and q, we can write this as (z-p)*(z-q) or (1-p*z^-1)*(1-q*z^-1) which allows you to readily transform a bunch of poles and zeroes (that is roots of denominator and numerator respectively) into a rational.

W3C publishes an audio EQ cookbook with the basic formulas.

It is possible to have instantaneous signals and there are even some variations in how those can be interpreted (e.g. Härmä suggests to understand it in the iterative sense, but it's also possible to understand it in some kind of bandlimited continuous sense, which allows convergence under a wider, actually all practically relevant, range of conditions).

When the output feeds back into the input, you have to calculate the output first before you can calculate the next sample. Without a delay somewhere, it just becomes an infinitely-long equation.

[mystran]

When the output feeds back into the input, you have to calculate the output first before you can calculate the next sample. Without a delay somewhere, it just becomes an infinitely-long equation.

But that's the thing.. it does NOT become an infinitely long equation. It becomes a system of equations that you can solve... and that solution is always finite for any stable system.

This is the "paradigm shift" I was talking about earlier: when you accept that you can solve systems of equations, you just don't need to care about any delay free loops anymore... and because you don't need to care about delay free loops, you don't need to design filters as carefully constructed signal flow graphs with no delay free loops, but rather you can just design them as systems of differential equations... and then the exact way in which you solve Ax=b is mostly matter of what results in the least number of arithmetic operations... but it's really the same filter whatever method you choose, 'cos the solver is just "implementation detail" now.

[Z1202]

When the output feeds back into the input, you have to calculate the output first before you can calculate the next sample. Without a delay somewhere, it just becomes an infinitely-long equation.

But that's the thing.. it does NOT become an infinitely long equation. It becomes a system of equations that you can solve... and that solution is always finite for any stable system.

And this solution even has causal interpretations, where the best one (that I so far know of) is to consider delay-free loops as some kind of bandlimited continuous-time loops. But that's if you want a deeper understanding, practically these loops converge to equation solutions in all typical situations.

[odibo]

I need to divide by f0/fs (thus multiply by fs*1/f0) to compute PolyBLEP.

Division is hard.

There's no really fast way to approximate 1/x to a low error.

Maybe your PolyBLEP oscillator can be implemented without division. From my math textbook:

A / B = exp ( log (A) - log (B) )

In your code:

fs / f0 = exp2 ( log2 (fs) - log2 (f0) )

exp2 would be a table lookup plus arithmetic shift, log2(fs) a precalculated constant, and log2(f0) should be already available because that's how frequency information usually gets passed around in synths.

[Tone2 Synthesizers]

Just had a (completely new?) idea that could work really well (extremely low CPU + great sound + no latency):

The idea is this: Basically the minimum phase PolyBLEP is replaced with a complex sine resonator that emits a 'ping' at nyquist frequency once the saw/square waveform jumps.

The complex IIR resonator is second order (12dB), set at very high resonance so that it decays rather slowly and resonance frequency is set to nyquist. Once the squarewave/sawwave jumps the fractional part of the sample offset is added/distributed to the complex resonator's current amplitude values.

[Z1202]

Just had a (completely new?) idea that could work really well (extremely low CPU + great sound + no latency):

The idea is this: Basically the minimum phase PolyBLEP is replaced with a complex sine resonator that emits a 'ping' at nyquist frequency once the saw/square waveform jumps.

The complex IIR resonator is second order (12dB), set at very high resonance so that it decays rather slowly and resonance frequency is set to nyquist. Once the squarewave/sawwave jumps the fractional part of the sample offset is added/distributed to the complex resonator's current amplitude values.

I understand the idea, and I had similar thoughts myself earlier. But I somewhat doubt it can work. It's just that BLEP is a bit too far from an exponentially decaying sine, even though it has a sine component at Nyquist asymptotically (but it doesn't decay exponentially). Would be happy to be proven other wise, tho.

[jrmoserbaltimore]

The idea is this: Basically the minimum phase PolyBLEP is replaced with a complex sine resonator that emits a 'ping' at nyquist frequency once the saw/square waveform jumps.

The complex IIR resonator is second order (12dB), set at very high resonance so that it decays rather slowly and resonance frequency is set to nyquist. Once the squarewave/sawwave jumps the fractional part of the sample offset is added/distributed to the complex resonator's current amplitude values.

Test it out. If it works, write a research paper and submit it to a journal.

[Tone2 Synthesizers]

Nope, since I don't get paid for this

[2DaT]

exp2 would be a table lookup plus arithmetic shift, log2(fs) a precalculated constant, and log2(f0) should be already available because that's how frequency information usually gets passed around in synths.

Division is the same, you can approximate reciprocals by lookup table and doing a NR iteration on that to improve upon. If I remember correctly, the iteration is 2 multiplies and 2 adds and even then, these don't have to be full precision if we don't need a precise result.

The incorrect assumption that division is impossibly hard and expensive comes from the fact that IEEE-754 compliant division has to be correctly rounded in the last bit, which implies more than 40 computed bits for a 23-bit single precision mantissa (exact amount and why that is you can find in the specialized literature).

If precise result isn't needed, 12-14 bits of precision are reasonable (Lookup+NR).

For polyblep in particular you can even try to get away without doing a NR iteration by doing a lookup without correction. It doesn't get easier than that.

[jrmoserbaltimore]

Storing a half-waveform for a band-limited saw is wild man.

Let's say you have an 8 sample period, indexes [0, 1, 2, 3, 4, 5, 6, 7] (sample 8 = sample 0)

When you reach sample 4, binary 3'b100, the value is 0. This is the same as sample 0. In fact, for each sample number:

s[4] = -s[0]
s[5] = -s[3]
s[6] = -s[2]
s[7] = -s[1]

That means you play back the 8 samples in the order [0, 1, 2, 3, 0, 3, 2, 1], or in binary:

3'b000
3'b001
3'b010
3'b011
3'b000 (phase: 3'b100)
3'b011 (phase: 3'b101)
3'b010 (phase: 3'b110)
3'b001 (phase: 3'b111)

Let's say we count up to 8 though, to track phase normalized to tau. That gives the numbers as noted.

Let's say if phase & 3'b100 is non-zero (single bit check), we use 4 minus the phase:

3'b000
3'b001
3'b010
3'b011
3'b000 (4 - 4)
3'b111 (4 - 5 = -1)
3'b110 (4 - 6 = -2)
3'b101 (4 - 7 = -3)

Notice the last two bits: they match up with the order in which we need to select samples. Since we overflow at 8 in this example, the next sample will just be sample 0 again.

This requires an extra adder, but that's a lot smaller than 1024 more 24-bit entries. In hardware this can be done in a roundabout way by only using the top (here, 4) bit of the phase as an input and wiring it to the corresponding input bit on the adder and the subtract signal. Basically, in this example, when the 3rd bit is on, the input is 4 - [phase], and when the 3rd bit is off, the input is 0 + [phase]. (In software you'd just use an if statement; flowing through the adder in both cases is faster in hardware.)

This was not as simple as I wanted it to be dammit. I've solved it now.

I also assume I'll need to find the largest sample value from all of my tables and divide all samples in all tables by that value to get all the band-limited saws to the same volume (the individual sine components will be the same amplitude, but there will be different numbers of them, so the waveform's peaks will be different). Mind you I failed algebra, check my transcript, you can see Elementary Algebra [D] and Calculus 1 And Analytical Geometry [A], so you might want to check my math….