Call for help regarding anti-aliasing

[MikeNakis]

@S0lo

1. Imagine that we can do near infinite sampling rate. We will not have aliasing at all !! Why do you think is that ?

Because it would not be "using samples at certain discrete intervals" anymore. It would be taking samples continuously. So the error is minimized, and in the case of an infinite sampling rate, eliminated. It is exactly the same as drawing a line on a graphics display with infinitesimally small pixels. Of course, there will be no need to anti-alias that line.

2. If the naive waveform generates any frequencies above half of the sampling rate, we will have aliasing. Why do think is that?

It is unclear to me what you mean by "the naive waveform generates frequencies". In any case, when sampling naively, aliasing will occur at any ratio of signal frequency versus sampling frequency, it just might not be perceptible when the ratio is very low.

For example, when sampling a 440 Hz square wave at a 44100 Hz sampling rate, each cycle of the wave will take about 100 samples to describe, but not exactly 100. So, the naive sampling approach will usually produce 50 samples at 1.0 followed by 50 samples at -1.0, but every once in a while it will produce 51 samples at 1.0, or 51 samples at -1.0. This inconsistency of sometimes 50 and sometimes 51 samples is aliasing. It is just not perceptible to our ears.

The higher we go up the frequency scale, the more perceptible aliasing becomes, and when we reach half the sampling rate the result is complete and utter garbage, but I am afraid it is wrong to say that aliasing only starts when we reach that frequency. Unless your DSP textbook specifically defines aliasing as such, which means that my definition of aliasing differs from your DSP textbook's definition, and my definition is useful, while the DSP textbook's definition is not.

@imrae
My definition says that aliasing is inevitable when sampling naively. (Barring other circumstantial factors, such as sampling a signal that happens to only contain slope changes at the precise intervals that the samples are taken.) That's distinctly different from "competently implemented". It is the competent implementation that I am after.

[MikeNakis]

@imrae furthermore, on the subject of competent implementation of sampling, I had an idea, which I thought should solve the problem, but it did not.

The idea was that instead of reading the wave function (y) at specific, infinitesimally narrow points in time (t), which is the naive approach, I would compute the surface area of the signal between the wave function and y=0, from the previous sample until the current sample. Then, I would convert the surface area to a scalar value by dividing by the time between samples.

I did a proof of concept with a sine wave, where the surface area can easily be computed using integration, which involves the anti-derivative of sin(t), which is known to be -cos(t). The results came out perfect, but of course this does not say much, because a sine wave does not suffer from aliasing anyway. The proof of concept just showed that there was nothing fundamentally wrong with the approach.

Then I tried with waveforms described using straight line segments, where the surface area between two sampling points can again be trivially computed since it is the surface area of a trapezoid in which 3 of the 4 sides are as in a rectangle, and only the top side is slanted.

The results sounded much better than the naive approach, but nowhere near as good as the results of using polyblep: there were still plenty of undesirable artifacts. I do not know why.

It is issues like this that I want to discuss.

[signalsmith]

The way that aliasing is usually considered here (which I do think is helpful) is to think about the spectrum of the continuous-time waveform before sampling. For something like a square wave, the hard edges mean it contains a long tail of high frequencies - higher than a given sample-rate can contain.

When we sample such a signal, those too-high frequencies cause problems (technically: appearing at other lower frequencies, without any harmonic relationship). So the goal is to somehow produce a continuous signal which doesn't have those too-high frequencies.

The integration approach you describe is equivalent to blurring with a sample-wide box-filter (rectangular moving average) before sampling. This reduces frequencies that are too high for the sampling rate, but not completely, and will slightly dampen the desired high frequencies as well.

While I get that you don't want to learn the entirety of all DSP knowledge just to solve this one problem, having enough understanding to see why that solution is imperfect will help design a better one, either by yourself or in discussion with others using a shared vocabulary.

[MikeNakis]

@signalsmith thanks for your reply.

You state that the integration approach is equivalent to blurring with a sample-wide box-filter, but is this in fact so?

As far as I understand, a box filter replaces values with the average of adjacent values. I suppose that a box filter with a width of one sampling period would produce each sample as the average of the previous naive sample and the current naive sample, as if the original (ideal) signal between those two values was a straight line of which we are taking the midpoint.

However, computing the *surface area* between the two sampling points on the horizontal and between the wave function and y=0 on the vertical (by means of integration or whatever) is intended to guarantee that we are accounting for all of the intricate variation of the ideal signal within that interval, and not just a straight line. At least that was my intention when I was coming up with this approach.

In any case, yes, it is that discussion with others using a shared vocabulary that I am seeking.

[signalsmith]

You state that the integration approach is equivalent to blurring with a sample-wide box-filter, but is this in fact so?

Just to be clear: I'm not talking about sampling the signal and then filtering. I'm talking about filtering/blurring the continuous-time signal, and then sampling. So it would include (but smooth out) all the sub-sample detail.

The goal is to get a (theoretical) continuous-time signal without the problematic too-high frequencies, so that when we sample that it doesn't cause problems. This approach doesn't completely do that, but it's a step in that direction which is why it sounds better.

This 1-sample-wide continuous-time box-filter/average is equivalent to your "surface area" approach. (It's also similar to some other useful techniques like DPW, which is usually explained with calculus.)

[imrae]

In a competent implentation of sampling the ADC includes a steep lowpass filter below the Nyquist frequency. So the square wave that is recorded does not actually contain vertical edges, and has some ripples near the steps. The reconstruction filter in the DAC fills in the corresponding edge from this information; it does not matter how close the mid-point lies to the sampling grid, the output signal can still match the (band-limited) input signal.

[kryptonaut]

An 'intuitive hand-wavy' description of Blep antialiasing (assuming a little prior knowledge):

A pure analog sawtooth or squarewave contains harmonics which theoretically extend to infinite frequency, decreasing in amplitude as the harmonic frequency increases.

These high harmonics cannot be represented in a sampled version of the wave, because they are higher than the Nyquist frequency (half the sample rate) and would be aliased to different frequencies.

If you could subtract all these above-Nyquist harmonics from the analog signal before sampling, then you'd be ok to sample without aliasing. The sum of these removed harmonics turns out to look like 'little packets of wiggliness' surrounding each step in the waveform - so if you remove one packet of wiggliness at each jump then you have antialiased that jump, by removing the high harmonics.

What Blep antialiasing does is to remove something that approximates to a packet of wiggliness, positioned with sub-sample accuracy, at each step in the waveform. This synthetic packet can be several samples long (generally using a lookup table) giving a fairly good approximation, or it can be a 2-sample polynomial PolyBlep giving a less good but acceptable approximation.

Integrating the waveform is worse than PolyBlep in that it doesn't really approximate the wiggles very well, but it may be better than nothing.

[Music Engineer]

I repeat, the problem with it is that it only works with a specific set of predefined waveforms, whereas what I want is a general purpose solution that will work with any waveform.

i only skimmed the thread but so far, but it seems to me that it is way too much focused on bleps whereas you probably should rather look into mip-mapping techniques in this case (which i did not yet see being mentioned - but i may have overlooked it). take your prototype waveform, create a bunch of lowpass-filtered versions of it with progressively less bandwidth (i.e. lower cutoff) and then during playback, select an appropriate version of the waveform depending on the pitch (higher pitches will select more heavily bandlimited versions of the waveform). i personally like to use a waveform length of 2048 samples, doing the filtering via FFT/iFFT, cutoffs going down in octaves, no downsampling of the low-bandwidth waves, 2x oversampling during playback, using linear interpolation. that's the setup, i have settled to in my own wavetable oscillators - used, for example, in straightliner. if i remember correctly, this should suppress the aliasing by amounts in the 100 dB ballpark

[S0lo]

@S0lo

It is unclear to me what you mean by "the naive waveform generates frequencies". In any case, when sampling naively, aliasing will occur at any ratio of signal frequency versus sampling frequency, it just might not be perceptible when the ratio is very low.

Yes, it happens with square, saw and other complex waveforms when constructed naively (ie. without regard to anti-aliasing) because they contain harmonics (sine waves) that go beyond half the sampling rate.

By the way, you can construct a saw tooth by adding many many sine waves with amplitudes that follow the 1/k formula (https://en.wikipedia.org/wiki/Sawtooth_wave)

If you stop before half sampling rate, you'll have no aliasing. It works, but it's very bad for CPU usage.

Edit: you can equivalently construct any waveform you want by adding arbitrary sine waves at desired harmonically related or in-harmonic frequencies as long as all frequencies don't exceed half sampling rate (nyquist). But again it's not practical for CPU. I think your better off using inverse FFT in this case.

[rafa1981]

i only skimmed the thread but so far, but it seems to me that it is way too much focused on bleps...

I was going to write that also, wavetable synthesis probably has to be considered too.

There are many implementations online to learn from if you learn better by reading code than papers (e.g. Surge, Vital?).

[JustinJ]

i only skimmed the thread but so far, but it seems to me that it is way too much focused on bleps whereas you probably should rather look into mip-mapping techniques in this case (which i did not yet see being mentioned - but i may have overlooked it). take your prototype waveform, create a bunch of lowpass-filtered versions of it with progressively less bandwidth (i.e. lower cutoff) and then during playback, select an appropriate version of the waveform depending on the pitch (higher pitches will select more heavily bandlimited versions of the waveform). i personally like to use a waveform length of 2048 samples, doing the filtering via FFT/iFFT, cutoffs going down in octaves, no downsampling of the low-bandwidth waves, 2x oversampling during playback, using linear interpolation. that's the setup, i have settled to in my own wavetable oscillators - used, for example, in straightliner. if i remember correctly, this should suppress the aliasing by amounts in the 100 dB ballpark

If you mip-map don't you end up with gaps before Nyquist or overflows that reflect back and alias? Currently I'm doing the FFT/iFFT calculation every 256 samples to take into account pitch mods etc. Seems to perform Ok and I'm sure I read somewhere that Urs of U-He was doing a similar thing with Zebra2 (which is why I tried it). The nice thing about this scheme is that you can monkey about with the waveform shape algorithmically and you kind of get it for almost nothing. You've also got the option of messing about with the partials before you iFFT back to a band-limited wave.

Goes down to well below 100dB aliasing wise and if you use Lagrange interpolation it goes even further. Very, very clean.

[mystran]

If you mip-map don't you end up with gaps before Nyquist or overflows that reflect back and alias? Currently I'm doing the FFT/iFFT calculation every 256 samples to take into account pitch mods etc.

In general, yeah, you do end up with gap and/or overflow, but if you have at least one mip per octave and oversample by 2x (which is what Music Engineer said he's doing) then you can pretty much keep all the gaps and overflows in the oversampled region and it's mostly a non-issue.

Dynamic FFT is another possibility and slightly more flexible as it allows you to also apply spectral filters before the IFFT, though the downside is that it's quite inefficient if you want a reasonable modulation rate. If you think about it, you're computing some table (let's say 2k samples which is a typical wavetable size choice) and only really using 256 (interpolated) samples out of it, right?

As it turns out, there's another option too, which in some cases actually turns out to be more efficient than FFT when you increase the moduleration rate: just do "naive" IDFT (where "naive" is in quotes, because it certainly needs to be heavily optimized with SIMD and recurrences).

That might sound crazy, until you realize that straight DFT doesn't have any dependency between different samples and you can compute it only for the points you need (which means we don't need interpolation as we can just compute the desired "fractional" positions directly, so quality tends to be better too).

This is effectively now additive synthesis, except you can do the trigonometric recurrences over frequency (rather than over time) which means you get to reset for every new sample (so no error accumulation over time; waveforms keep their form as harmonics won't drift).

Ofcourse, like everything there's a downside to this method as well: if you want unison, you need to run the computation separately for every unison voice. With an FFT, you can reuse the same table which is cheap (even though it can sound a bit cheesy with spectral filters shifting with detune). With straight DFT on the other hand 8 voice unison means 8 times the CPU load; this isn't bad for high frequencies where you don't have the many harmonics anyway (since you can skip the DFT computation once past Nyquist), but it's arguably a bit ugly when you have the combination of lots of unison and low frequencies.

The bottom line though is that there are quite a few techniques and a lot depends on both the set of features you want and the quality vs. CPU cost tradeoff you're willing to make.