Investigating the math of waveshapers: Chebyshev polynomials
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 Teensybased 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 bitcrushers–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 NonLinear 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!

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

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^21$
$T_3(x) = 4x^33x$
$T_4(x) = 8x^48x^2+1$
$T_5(x) = 16x^520x^3+5x$
and more Chebyshev polynomials can be derived from the recursive formula
\[T_{n+1}(x) = 2 x T_n(x)T_{n1}(x)\] 
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)\] 
Normalize $f_1(x)$ by (bruteforce) 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!

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$.

The sum reduces to a single Chebyshev polynomial term, so the preliminary function is
\[f_0(x) = x + 0.2 (2x^21)\] 
We can calculate that $f_0(0)=0.2$, so our new function must be
\[f_1(x) = x+0.2 (2x^21)+0.2\] 
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^21)+0.2}{1.4}\]
That function simplifies to $\frac{2}{7}x^2+\frac{5}{7}x$.
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 nonzero 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 doubleangle identity $\cos(2x)=2 \cos^2(x)1$. Now if you will, imagine unsubstituting $\cos(x)$, and you get the Chebyshev polynomial $T_2(x) = 2x^21$. Then, imagine a doubledoubleangle formula, $\cos(4x)=2 \cos^2(2x)1$, and expand that to $8 \cos^4(x)8 \cos^2(x)+1$. Unsubstituting $\cos(x)$ gives the Chebyshev polynomial $T_4(x)=8x^48x^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_{n1}(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 tripleangle formula!
Okay, my only reason for bringing up $T_4(x)$ was this intuitivelooking 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.