Over a year ago, I wrote “Adding harmonic distortions with Arduino Teensy”. In that post, I happened upon a way to apply any arbitrary profile of harmonics using a Teensy-based waveshaper (just except that waveshapers categorically can’t vary the phase of each harmonic). I used that project to replicate the “musical” distortion of tube amplifiers, but I had forgotten about the other purpose of those things: sound synthesis. Guitar pedals, soft clippers, and bit-crushers–that is, pushing distortion to the point of musicality in its own right. That was where all the established literature was! Even in 1979, there was “A Tutorial on Non-Linear Distortion or Waveshaping Synthesis”.

I happened upon the same thing in an overcomplicated way, but it helped me learn the practical nuances of waveshapers. Let me show a mix of my own and the established method. Compared to what I did before, it’s easier to implement and much more concise.

How to generate a waveshaper function in four steps!

  1. Decide what amplitude ratios $\alpha_n$ each $n$-th harmonic should have with the fundemental frequency

  2. Build a preliminary function $f_0(x)$ as the weighted sum of the Chebyshev polynomials

    \[f_0(x) = T_1(x) + \sum_{n=2}^N \alpha_n T_n(x)\]

    where the first Chebyshev polynomials are

    $T_0(x) = 1$

    $T_1(x) = x$

    $T_2(x) = 2x^2-1$

    $T_3(x) = 4x^3-3x$

    $T_4(x) = 8x^4-8x^2+1$

    $T_5(x) = 16x^5-20x^3+5x$

    and more Chebyshev polynomials can be derived from the recursive formula

    \[T_{n+1}(x) = 2 x T_n(x)-T_{n-1}(x)\]

  3. Shift $f_0(x)$ so that it maps zero to zero (for preventing constant DC) by evaluating $f_0(x)$ at $x=0$ then subtracting that

    \[f_1(x) = f_0(x)-f_0(0)\]

  4. Normalize $f_1(x)$ by (brute-force) finding the maximum absolute value for $-1 < x < 1$ then dividing by that

    \[f_2(x) = \frac{f_1(x)}{f_{\text{1,absmax}}}\]

The above function, $f_2(x)$, is your final function. Evaluate it at as many points within $-1 < x < 1$ as can fit in your waveshaper’s LUT! If the input sine wave swings exactly within $-1 < x < 1$, then the ratios $\alpha_n$ will be realized. Otherwise, different and smaller ratios will occur.

Using this method, I can perfectly replicate my old post!

  1. In that old post, I chose to give the second harmonic a weight of 0.2 and no weight to the higher ones, so $\alpha_2 = 0.2$ and $\alpha_n = 0$ for $n > 2$.

  2. The sum reduces to a single Chebyshev polynomial term, so the preliminary function is

    \[f_0(x) = x + 0.2 (2x^2-1)\]

  3. We can calculate that $f_0(0)=-0.2$, so our new function must be

    \[f_1(x) = x+0.2 (2x^2-1)+0.2\]

  4. Plotting $f_1(x)$ reveals that it achieves a maximum absolute value of 1.4 at $x=1$, so our final function must be

    \[f_2(x) = \frac{x+0.2*(2x^2-1)+0.2}{1.4}\]

That function simplifies to $\frac{2}{7}x^2+\frac{5}{7}x$.

New and old plots

What’s new? What’s established?

Steps 1 and 2 are inspired by that 1979 paper, but steps 3 and 4 were my own. Step 4 was a practical move. It responds to the Teensy waveshaper only accepting points between -1 and 1. If that isn’t true for your waveshaper, then you’re at liberty to do something else–as long as you don’t clip your signal unintentionally! On the other hand, step 3 represents a constraint that I haven’t seen online. Specifically, we shouldn’t want to map zero to something non-zero because that would cause constant DC! This is why I subtracted $f_0(0)$.

However, I think it’s interesting to ask why we would have a constant DC problem in the first place. That requires digging into step 2. Why do we take a weighted sum of the Chebyshev polynomials? What are Chebyshev polynomials? And from there, how could they break?

Applying Chebyshev polynomials in waveshaping?

Chebyshev polynomials can be used in a rigorous approach to building waveshapers. In this context, their claim to fame is that they can twist a $\cos(x)$ wave into its $n$-th harmonic, $\cos(nx)$. Somehow, high school algebra and trigonometry are all you need to use them.

You already know one if you remember the double-angle identity $\cos(2x)=2 \cos^2(x)-1$. Now if you will, imagine un-substituting $\cos(x)$, and you get the Chebyshev polynomial $T_2(x) = 2x^2-1$. Then, imagine a double-double-angle formula, $\cos(4x)=2 \cos^2(2x)-1$, and expand that to $8 \cos^4(x)-8 \cos^2(x)+1$. Un-substituting $\cos(x)$ gives the Chebyshev polynomial $T_4(x)=8x^4-8x^2+1$. However, there is another way to jump from $T_2(x)$ to $T_4(x)$, and that is the recursion formula $T_{n+1}(x) = 2 x T_n(x)-T_{n-1}(x)$. Feel free to try it yourself, and notice that this alternative method involves stepping through $T_3(x)$ first. That would have corresponded to the triple-angle formula!

T_4(x) and cos(4x) plots

Okay, my only reason for bringing up $T_4(x)$ was this intuitive-looking plot. Phase shifts usually get in the way of the intuitive nature, but not for $n=4$. That aside!

At the same time, how we should use a Chebyshev polynomial becomes clearer. That is, they’re actually functions of $\cos(x)$ that spit out $\cos(n x)$! In symbols, $T_n(\cos(x)) = \cos(n x)$.

The final product of waveshapers should be a sum of the original wave plus some new harmonics. According to $T_n(\cos(x)) = \cos(n x)$, the Chebyshev polynomials stand in for the new harmonics, and so we simply take a weighted sum.

\[f_0(x) = T_1(x) + \sum_{n=2}^N \alpha_n T_n(x)\]

Here $T_1(x) = \cos(x)$ by defintion, so it is our original wave. Then, each $\alpha_n T_n(x)$ term adds the $n$-th harmonic, scaled by $\alpha_n$. That is what step 2 really means.

Breaking polynomial waveshapers with arbitrary inputs

That said, something is wrong with the result of step 2. Often, we’d want to input a wave of some varying amplitude $a(t) \leq 1$ (i.e. an ADSR envelope) or even an arbitrary input. I don’t know where the impacts end. At the very least, I can address one of them: DC shifts.

For $a(t) = 0$, a waveshaper will see nothing but zero, and it may decide to map that to something nonzero. This is because Chebyshev polynomials weren’t defined with that in mind. For example, $T_2(0)=-1$. If my headphones saw -1 volts at DC, they’d blow. From my old post, I only saw that happen when I added even harmonics, and I’ve seen that adding a constant equal to $\alpha_n$ or $-\alpha_n$ would correct that. Ultimately though, the easiest way to correct this effect is to just evaluate the waveshaper function at $x=0$, then subtract that value. That’s step 3. Finally, this addition technically could be considered as a part of the weighted sum in part 2. Specifically, it would be as $\alpha_0$ of another Chebyshev polynomial, $T_0(x)=1$, that represents a “zeroth” harmonic or 0 Hz. I just wouldn’t wish upon anyone the task of solving for $\alpha_0$ without evaluating the rest of the sum first.