# Is there such a thing as band-limited non-linear distortion?

So if you generate a square wave by just switching a signal between two values, at sample boundaries, it produces an infinite series of harmonics, which alias and produce tones below your fundamental, which is very audible. The solution is Band-Limited Synthesis, either using additive synthesis or band-limited steps to produce waveforms that are the same as if you had band-limited the ideal mathematical square wave before sampling it:

http://flic.kr/p/83JMjT

But I just realized that if you apply large amplification to a digital sine wave and then clip it digitally, it will produce the same square wave shape, without the Gibbs phenomenon ripples. So it's then also producing aliased distortion products, right? So any non-linear distortion in the digital domain that produces harmonics outside of the Nyquist limits will produce aliased distortion products? (Edit: I've done some tests and confirmed that this part is true.)

Is there such a thing as band-limited distortion, to simulate (in the digital domain) the effects of distorting (in the analog domain) before band-limiting and sampling? If so, how do you do it? If I search for "bandlimited distortion" I find some references to Chebyshev polynomials, but I don't know how to use them or if they only work for sine waves or what:

This instrument does not attempt to generate band-limited distortion. Those interested in band-limited distortion should investigate the use of Chebyshev polynomials to generate the effect. Hyperbolic Tangent Distortion

"Chebyshev polynomial" -- shaping functions with the important property that they are intrinsically band-limited i.e. they do not introduce spurious spectral harmonics due to overlapping etc. Wave Shaper

• I'm not sure what you're asking. If you perform an operation that causes generation of frequency content outside the Nyquist region that you're operating in, then you will see aliasing, regardless of how you generated said content. What kind of analog distortion are you trying to simulate? One approach could be to upsample the signal to a sufficiently high sample rate first, then using the wider Nyquist region to perform your signal processing. You could then downsample back to the original rate after you're done. – Jason R Jun 11 '12 at 15:53
• @JasonR: Yes, for generating square waves, you can either do a truly band-limited method like additive synthesis, or you can approximate by upsampling first, generating the square wave in the naive way, and then downsampling (but there will still be some aliasing, just at a lower level). Likewise, you can approximate the distortion as you've said by upsampling first, but is there a way to generate it directly, with zero aliasing, analogous to the additive synthesis method for square wave generation? – endolith Jun 11 '12 at 15:57
• @JasonR: I'm asking about any non-linear distortion, in general, but something like emulating an analog guitar amp's distortion circuitry would be a good example. If I understand correctly, doing it naively in digital domain would produce distortion products that don't exist when distorted in the analog domain, some of which might be clearly audible at frequencies lower than the fundamental, etc. – endolith Jun 11 '12 at 16:00
• @endolith Chebyshev polynomials may be what you want. – datageist Jun 11 '12 at 19:36
• @datageist: Can you write an answer on how to use Chebyshev polynomials? Even a short one. – endolith Apr 15 '13 at 16:32

Applying a non-linear function will always introduce harmonics, and mixing non-linear functions with sampled versions of continuous signals does add the wrinkle you note above (where high-frequency harmonics are aliased to low frequencies.)

I can think of a few ways to proceed:

1. You can use an oversampling factor high enough to capture the extra harmonics (up to some arbitrary precision, e.g. your noise floor),
2. You can use a "softer" clipping function (see, for example, here) that has harmonics that die out sooner than the hard clipper. This is easier to model, but introduces its own distortion at low frequencies.
3. Building on the approach you suggested above, interpolate your sampled signal (e.g. using a Lagrange or Chebyshev interpolator) to construct a continuous-time model. Then, apply the hard clipper and low-pass in a simulated continuous-time domain. Sample the result.

You can combine (1) and (2). The third approach is complex, but gives you the best control over how much distortion to admit, and will probably scale better to very high fidelity requirements.

For non-linear functions that admit a series expansion (e.g. Taylor/Maclaurin), you can get a decent intuition for how fast the harmonics decay. The Maclaurin expansion of a function $f(x)$ is:

$$f(x)=\sum_{n=0}^\infty \left[ \frac{f^{(n)}(0)}{n!} x^n\right]$$

In your case, $f(x)$ is the clipping function. (You can't do this with a hard clipper, at least not naively!) If you consider the substitution $x=g(t)$, where $g(t)$ is your input signal, the $x^n$ becomes $g(t)^n$, which you can consider to be the convolution of your input signal with itself $n$ times. Thus, for low-pass signals, the $n$th term of the infinite summation has a bandwidth $n$ times that of your signal. To complete the picture, you need to figure out the amplitude associated with each term and decide how many terms in the summation are relevant.

(With a little thought, you might also be able to use this form directly to approximate the filtered non-linearity. That would require a good series representation for the clipper.)

• To clarify, #3 is not just oversampling with interpolation, it's finding parameters of a continuous Chebyshev polynomial that fits the sampled points and then working with those parameters and a model of the polynomial? – endolith Jun 12 '12 at 18:49
• I'm imagining a sequence of polynomial interpolators that are each active on a short range of samples. So, as a batch of new samples comes in, you'd construct an interpolator that's only active in a defined interval. Your continuous-time approximation to the sampled signal consists of these polynomials. (I'm thinking Lagrange, but Chebyshev is probably the same thing. I don't recall if Chebyshev interpolators match the sample points exactly. If not, you'd get discontinuities when switching between interpolators.) – Graeme Jun 12 '12 at 19:36

A few approaches to alias-free nonlinear distortion (in increasing order of difficulty):

1. Subband distortion: Use a low pass filter to extract the lower end of the signal. If you choose a cutoff frequency of $\frac{f_s}{2N}$ you can apply any non-linear transfer function $f$ with derivatives starting at $f^{N+1}$ vanishing to avoid aliasing. Add just the distortion output of the low band back to the original signal. For this, modify your nonlinear transfer function to only produce the distortion, not the incoming signal. This approach is often totally sufficient for audio, because the distortion contribution of the higher frequencies would be inaudible in most cases.

2. Oversampling: This is similar to approach 1, but does require the extra steps of upsampling and downsampling to make sure that all information in the original signal band contributes to the distorted signal, if you think this is needed. It also allows you to distort to higher order, because parts of the back-folded spectrum will be removed by the downsampling filter. That means if you oversample by a factor $N$, the derivatives of the nonlinear transfer function starting at order $2N$ have to vanish. The procedure is a simple as upsampling, applying the nonlinear transfer function, downsampling.

3. Using a local analytic solution: Any smooth nonlinear transfer function can be applied to a signal that is described by a power series. If you let the transfer function act on a power series, you will get a power series back. If you restrict yourself to truncated power series of a certain order, and the transfer function to a certain degree of smoothness, you can write the transfer function as a map on the coefficients of the truncated series. That means instead of $f:\mathbb{R}\to\mathbb{R}$, you consider $f:\mathbb{R}^N\to\mathbb{R}^M$ taking an input series expansion of length $N$ to an output series expansion of length $M>N$. With this understanding, you can use a local finite approximation of the input signal in terms of a power series and map it to a power series approximation of the output. Then you can analytically integrate the output series to create a box-car antialiasing filter to get the output sample value. All these computations can be done symbolically, and because you have to include local features of the input signal, you will finally get a nonlinear filter that uses past values of your input signal to generate the current output.

4. Constrain based algebraic design: In the previous item, you have seen that antialiasing nonlinear distortion leads to nonlinear filters. Of course, not all nonlinear filters are alias free, but some may be. So the obvious question is for a criterium to make such a filter strictly alias free and how to design it. As it turns out, an equivalent statement to being free of aliasing is, that the non-linear filter commutes with sub-sample translations. So you must make sure that it does not make a difference if you translate first and then filter, or filter first and then translate. This condition leads to very strict design constraints for nonlinear filters, but depends on how you realise the signal translation. For example, the ideal translation would require infinitely many coefficients for the nonlinear filter. So you have to approximate signal translation to finite order to get a finite nonlinear filter. Alias-freeness scales with the approximation you use, but you have very good control over it. After you've worked through the mathematics of this approach, you can design any (not just smooth) nonlinear transfer function as a nearly ideal digital model in the form of a nonlinear filter. I cannot possibly sketch the details here, but maybe you can find some inspiration from this description.

• Do you know what the alias-free version of the non-hysteretic monomial transfer functions would be? In other words, the alias-free version of $y[t] = x[t]^k$ for some positive integer k. Or do you know of any published work on this topic? – Mike Battaglia Aug 22 '18 at 5:38
• Or, another (perhaps related) question - if you go with the local finite approach that you mentioned, you get a map from truncated power series to truncated power series. When you then attempt to lowpass the truncated power series, which would normally be convolving with a sinc function, do you get any simple expression for the result? Can the result be expressed again as a truncated power series, and if so, what does it look like? – Mike Battaglia Aug 22 '18 at 15:29
• @MikeBattaglia, maybe you can create a new question so that I can give a thorough answer there. To answer your second question, you don't use a SINC kernel but in the simplest case a boxcar kernel. Higher order kernels do work, but in order to get an analytic result you have to restrict yourself to certain kernels. – Jazzmaniac Aug 22 '18 at 15:50
• Thanks - created a new question here for the filter design question: dsp.stackexchange.com/q/51533/18276 – Mike Battaglia Aug 27 '18 at 3:31

One algorithm that qualifies as bandlimited distortion is the use of Chebyshev polynomials with sinusoidal inputs for waveshaping synthesis. The Chebyshev polynomials [of the first kind] can be defined as $$T_n(x) = \mathrm{cos}(n \, \mathrm{arccos}(x)).$$

I won't go into details here, but there's a relatively simple derivation using De Moivre's Theorem which shows how the above formula defines a set of polynomials (cf. Sec. 7.2 in Digital Filters by Hamming). The trigonometric definition is useful because it makes it easy to see how each $T_n(x)$ maps a unit amplitude cosine into its nth harmonic, i.e.

$$T_n(\mathrm{cos}(kx)) = \mathrm{cos}(n\,\mathrm{arccos}(\mathrm{cos}(kx))) = \mathrm{cos}(nkx). \tag{1}$$

The polynomials themselves can be easily generated by using the following recurrence relation:

$$T_0(x) = 1 \\ T_1(x) = x \\ T_n(x) = 2x\,T_{n-1}(x) - T_{n-2}(x).$$

Here are the first few:

\eqalign{ T_0(x) &= 1 \\ T_1(x) &= x \\ T_2(x) = 2x(x - 1) &= 2x^2 - 1 \\ T_3(x) = 2x(2x^2 - 1) - x &= 4x^3 - 3x \\ T_4(x) = 2x(4x^3 - 3x) - (2x^2 - 1) &= 8x^4 - 8x^2 + 1 \\ \ldots }

Another way to gain intuition for property $(1)$ is to consider, for example, the operation of $T_2$ on the complex exponential form of $\mathrm{cos}(x)$

\eqalign{2\,\mathrm{cos}^2(x) - 1 &= 2\left(\frac{e^{ix} + e^{-ix}}{2}\right)^2 - 1 \\ &= \frac{2}{4}(e^{i2x} + 2e^{ix}e^{-ix} + e^{-i2x}) - 1 \\ &= \left(\frac{e^{i2x} + e^{-i2x}}{2}\right) + \frac{2}{2} - 1 \\ &= \mathrm{cos}(2x). }

By computing a Chebyshev Series

$$f(x) = \sum_{n = 0}^{\infty}{a_n\,T_n(x)}$$

truncated to an appropriate value of $n$, you get a waveshaper (i.e. $f(x)$) that, when applied to a unit amplitude cosine, will generate an arbitrary [integer] set of bandlimited harmonics.

• Thanks! For waveforms other than a single sinusoid, what happens? Bandlimited intermodulation or not? – endolith Apr 21 '13 at 17:10
• Yeah, you get intermodulation, but it is bandlimited. For most full-bandwidth audio, that just means you know the oversampling factor you need to use. And it will be true of all polynomial waveshapers up to $x^n$. – datageist Apr 21 '13 at 21:14

@robert-bristow-johnson explains this very clearly on comp.dsp:

you have to oversample to a finite extent. if you represent the (memoryless, i assume) non-linearity as a finite order polynomial (that approximates whatever curve you are trying to implement), then whatever is the order of polynomial is the same factor of oversampling needed and no aliases will occur. then low-pass filter (at that oversampled rate) to get rid of all frequency components higher than your original Nyquist, then downsample and you will not have aliasing.

In other words, if your nonlinearity is a polynomial, the highest frequency that can be produced by the distortion will be the highest frequency in your signal times the order N of the polynomial. (The polynomial nonlinearity is multiplying the signal by itself N times, so its spectrum gets convolved with itself and spreads out by the same ratio.)

So then you know the maximum frequency (whether Nyquist or some lower limit for your application), and you know the order of the polynomial, so you can oversample enough to prevent aliasing, do the distortion, and then low-pass filter and downsample.

In fact, you can reduce the oversampling rate by letting some aliasing happen, as long as it's contained in the band that will be removed before downsampling:

another little trick is that you need not care about aliasing that folds down to the area that you'll LPF out. so a 5th order polynomial needs only to have an oversampling ratio of 3. those top 2 harmonics might alias, but won't get back into the baseband. when downsampling, you filter those aliased harmonics out. so i think the hard and fast rule is

oversampling ratio = (polynomial order + 1)/2

• i seem to have overstated the oversampling requirement is a little in that post. you can allow some foldover of harmonics as long as the aliased harmonics don't make it back down to the original baseband. that means, if the order of the polynomial representing the memoryless non-linearity is $N$, all you need to oversample by is a factor of $\frac{N+1}{2}$. e.g. if you oversample by a factor of 4 (like upsample to 192 kHz), you can apply a 7th order polynomial, filter out (at $f_s$=192 kHz) all the crap generated above 20 kHz, and then downsample back to 48 kHz. – robert bristow-johnson Feb 10 '16 at 6:37
• @robert I saw and added that part as you were commenting. If you rewrite this as your own answer I'll delete it and accept yours – endolith Feb 10 '16 at 6:39
• oh dear. to dress up a nice answer requires work. – robert bristow-johnson Feb 10 '16 at 6:43
• @robertbristow-johnson Well, the imaginary internet points are yours if you want them – endolith Feb 10 '16 at 13:42