Can harmonics be created or simulated a la carte

Since the input can be any type of Audio from 20-20k, filtering is not likely the idea. I understand to some limited degree the methods of adding harmonic distortion, but have not been encouraged towards a la carte control. The end result is to custom tailor a harmonic distortion effect.

Thanks

Without re-hashing all of the maths from RBJs response, I made a crummy little plug in that did more/less what you are proposing. It had a stupid simple GUI with one slide for nominal amplitude and several more to control harmonic levels. It was OK, but not as useful as one would hope for audio work. The idea was that given a sinusoid with a nominal input amplitude, it would calculate the polynomial required for the given harmonics. The bummer is that audio doesn’t work that way, at all. The sound doesn’t have any more harmonic ‘richness’ if you can call it that, it just sounds distorted in different ways.

To wit, my peripheral understanding of products advertising harmonic distortion use crossover filters to divide the signal into different bands, on which distortion is applied. By narrowing the bandwidth, you get something more closely approximating a sinusoid. Again though, my impression was that they were fine, but not exactly what I had in mind.

My experience was that as simple as ‘adding harmonics’ seems, it’s a pretty tough thing to implement. My apologies for not providing a more rigorous analysis, but figured some first hand experience might be relevant.

• Thanks Dan. I appreciate the effort to test an idea. Jan 3 '21 at 5:36

The answer is "sorta". Lot's of assumptions, though. This is related to that other question where you've been.

This is about waveshaping of a single sinusoid (without harmonics) by performing a polynomial operation on a sinusoid having an amplitude of 1 and any frequency. It is memoryless, which will mean that the waveshaping operation is (at least to start with) does not include filtering and has no states like a filter does. There may be cases where you want to filter either the input and/or output, to modify the amplitude of the harmonics based on pitch of the note.

Leaving the filtering aside for another big topic, this polynomial mapping function is useful that you can control what the highest harmonic generated is and limit aliasing. Aliased waveshaping often sounds cheap and digital, not nice and analog. But here is the polynomial (having order $$N$$) functional form, and a good way to implement it with Horner's rule:

\begin{align} p_N(x) &= \sum\limits_{n=0}^{N} \ a_n \ x^n \\ \\ &= a_{\small{0}} + \Bigg(a_{\small{1}} + \bigg(a_{\small{2}} + \Big(a_{\small{3}} + \,... \big(a_{\small{N-2}} + (a_{\small{N-1}} + a_{\small{N}} \,x \,)x \, \big)x \ ...\Big)x \, \bigg)x \, \Bigg)x\\ \end{align}

So if the input $$x(t)$$ is a sinusoid of frequency $$f_0$$, the highest possible harmonic you will get outa this is at frequency $$N\,f_0$$. If it's an even-symmetry polynomial ($$a_n=0 \ \text{if}\ n\ \text{odd}$$), then only even-numbered harmonics will come out. Similarly for odd-symmetry polynomials ($$a_n=0 \ \text{if}\ n\ \text{even}$$).

Now the next thing you gotta learn about are Tchebyshev polynomials which can be defined as:

\begin{align} T_m(x) & = \begin{cases} \cos(m\arccos x) \qquad & \text{ for }~ |x| \le 1 \\ \\ \dfrac{1}{2} \bigg( \Big(x-\sqrt{x^2-1} \Big)^m + \Big(x+\sqrt{x^2-1} \Big)^m \bigg) \qquad & \text{ for }~ |x| \ge 1 \\ \end{cases} \\ \\ & = \begin{cases} \cos(m\arccos x) \qquad \quad & \text{ for }~ -1 \le x \le 1 \\ \\ \cosh(m \operatorname{arcosh}x) \qquad \quad & \text{ for }~ 1 \le x \\ \\ (-1)^m \cosh\big(m \operatorname{arcosh}(-x)\big) \qquad \quad & \text{ for }~ x \le -1 \\ \end{cases} \\ \end{align}

and then results in this recursion rule which proves that they are polynomials and one way to get the polynomial coefficients, $$a_n$$. The functions above surely don't look like polynomials, but they are. For orders $$m=0$$ and $$m=1$$, it is clear that

\begin{align} T_0(x) & = 1 \\ T_1(x) & = x \\ \end{align}

And from analysis of the definition of $$T_m(x)$$, this recursion formula is can be shown to apply:

$$T_{m+1}(x) = 2 x\,T_m(x) - T_{m-1}(x)$$

which shows how the polynomial order increases as $$m$$ does and how to build the Tchebyshev polynomial.

Now the first thing to notice about these Tchebyshev polynomials is what they do with a sinusoid having an amplitude of 1 being input to the polynomial function.

$$x(t) = \cos(2 \pi f_0 t)$$

What comes out of the $$m$$-th order Tchebyshev polynomial is:

\begin{align} y_m(t) &= T_m\big( x(t) \big) \\ &= T_m\big( \cos(2 \pi f_0 t) \big) \\ &= \cos\Big(m \arccos\big( x(t) \big) \Big) \\ &= \cos\Big(m \arccos\big( \cos(2 \pi f_0 t) \big) \Big) \\ &= \cos\big(m (2 \pi f_0 t) \big) \\ &= \cos\big(2 \pi (m \,f_0) t \big) \\ \end{align}

So what goes in is a unit-amplitude sinusoid of frequency $$f_0$$ and what comes out is a unit-amplitude sinusoid of frequency $$m\, f_0$$, the $$m$$-th harmonic. The amplitude of that $$m$$-th harmonic can be scaled with its own coefficient, $$b_m$$.

$$b_m \cos\big(2 \pi (m \,f_0) t \big) = b_m T_m\big( \cos(2 \pi f_0 t) \big)$$

To create a variety of up to $$N$$ harmonics each with their independent amplitude, you can just sum it up.

\begin{align} y(t) &= \sum\limits_{m=0}^{N} b_m \cos\big(2 \pi (m \,f_0) t \big) \\ &= \sum\limits_{m=0}^{N} b_m T_m\big( \cos(2 \pi f_0 t) \big) \\ &= \sum\limits_{m=0}^{N} b_m T_m\big( x(t) \big) \\ &= \sum\limits_{n=0}^{N} a_n \ \big( x(t) \big)^n \\ & = p_N \big( x(t) \big) \end{align}

So you simply pass your input sinusoid (of unit amplitude) through this simple polynomial function $$p_N(x)$$ having coefficients $$a_n$$ and what comes out is a sum of sinusoids all having their own amplitudes of $$b_m$$.

Since $$T_m(x)$$ is a polynomial of order $$m$$, adding $$N$$ of these polynomials having orders of $$0$$ up to $$N$$, will result in an aggregate polynomial of order $$N$$ that will have coefficients of $$a_n$$. The next problem to solve is how to determine the coefficients of $$p_N(x)$$ (which we're calling $$a_n$$) from the known harmonic coefficients $$b_m$$.

So here we will express the power-series terms of the Tchebyshev polynomial explicitly as:

$$T_m(x) = \sum\limits_{n=0}^{m} \ a_{m,n} \ x^n$$

$$a_{m,n}$$ is the coefficient that multiplies the $$n$$-th power term, $$x^n$$, in the $$m$$-th Tchebyshev polynomial $$T_m(x)$$. $$a_{m,n}$$ is defined only for $$0 \le n \le m$$. We can assume that for all $$n>m$$, then $$a_{m,n}=0$$. For even $$m$$, then all odd-order ($$n$$ odd) coefficients will be zero, and likewise for $$m$$ odd. Because

\begin{align} T_0(x) & = 1 \\ T_1(x) & = x \\ \end{align}

then we know that

\begin{align} a_{0,0} & = 1 \\ a_{1,0} & = 0 \\ a_{1,1} & = 1 \ . \\ \end{align}

Then we can use the salient Tchebyshev recursion formula above to determine the polynomial coefficients for Tchebyshev polynomial coefficients.

\begin{align} T_{m+1}(x) &= 2 x\,T_m(x) \ - \ T_{m-1}(x) \\ \\ \sum\limits_{n=0}^{m+1} a_{m+1,n} \ x^n &= 2 x\,\left(\sum\limits_{n=0}^{m} a_{m,n} x^n \right) \ - \ \sum\limits_{n=0}^{m-1} a_{m-1,n} x^n \\ &= \sum\limits_{n=0}^{m} 2 a_{m,n} x^{n+1} \ - \ \sum\limits_{n=0}^{m-1} a_{m-1,n} x^n \\ &= \sum\limits_{n=1}^{m+1} 2 a_{m,n-1} \ x^n \ - \ \sum\limits_{n=0}^{m-1} a_{m-1,n} \ x^n \\ &= \sum\limits_{n=1}^{m+1} 2 a_{m,n-1} \ x^n - \sum\limits_{n=1}^{m-1} a_{m-1,n} \ x^n - a_{m-1,0} \\ &= \sum\limits_{n=1}^{m-1} 2 a_{m,n-1} x^n + 2 a_{m,m-1} \ x^m + 2 a_{m,m} x^{m+1} - \sum\limits_{n=1}^{m-1} a_{m-1,n} x^n - a_{m-1,0} \\ &= \sum\limits_{n=1}^{m-1} (2 a_{m,n-1}-a_{m-1,n}) x^n \ + 2 a_{m,m-1} x^m + 2 a_{m,m} x^{m+1} - a_{m-1,0} \\ \end{align}

Comparing terms having the same power of $$x$$, we can see

\begin{align} a_{m+1,0} &= -a_{m-1,0} \\ a_{m+1,n} x^n &= (2 a_{m,n-1}-a_{m-1,n}) x^n \qquad & \text{for} \ 1 \le n \le m-1 \\ a_{m+1,m} x^m &= 2 a_{m,m-1} x^m \\ a_{m+1,m+1} x^{m+1} &= 2 a_{m,m} x^{m+1} \\ \end{align}

or

$$a_{m+1,n} = \begin{cases} -a_{m-1,n} \qquad & n=0 \\ 2 a_{m,n-1}-a_{m-1,n} \qquad & 1 \le n \le m-1 \\ 2 a_{m,n-1} \qquad & m \le n \le m+1 \\ \end{cases}$$

Now, in the end, your net polynomial coefficients, $$a_n$$, with the weightings of harmonic amplitudes, $$b_m$$, are

\begin{align} y & = p_N(x) \\ &= \sum\limits_{n=0}^{N} a_n \ x^n \\ &= \sum\limits_{m=0}^{N} b_m T_m(x) \\ &= \sum\limits_{m=0}^{N} b_m \ \sum\limits_{n=0}^{m} a_{m,n} \ x^n \\ &= \sum\limits_{m=0}^{N} \sum\limits_{n=0}^{m} b_m \ a_{m,n} \ x^n \\ &= \sum\limits_{n=0}^{N} \sum\limits_{m=n}^{N} b_m \ a_{m,n} \ x^n \\ \end{align}

Then the net coefficient that multiplies $$x^n$$ is the sum of all the Tchebyshev coefficients for the same power, $$x^n$$, in the Tchebyshev polynomials each weighted by the corresponding $$b_m$$

$$a_n = \sum\limits_{m=n}^{N} b_m a_{m,n} \qquad 0 \le n \le N$$

• Thanks Robert for your time to reply with this great detail. It’s over my head but a quick question is does the concept apply to a range of inputs ie vocal or full mix or instruments? I notice your reference to “single sinusoid” and the application would rarely ever be a pure sine but rather complex input. Jan 1 '21 at 5:05
• listen, i'll continue and complete an answer anyway. you actually can tune a polynomial to take a sinusoid as input and return a set of harmonics with prescribed amplitudes. for a sine wave input, that can be done. now how it works for some other waveform going in (like a monophonic guitar note) that is less clear how you would tune the harmonics for it. Jan 1 '21 at 8:20
• If your input is going to be any arbitrary sound (even any arbitrary music), please clarify that in your question. Jan 1 '21 at 19:05
• I wanted to lay out the math explicitly. I can maybe show you a short snippet of C code that will compute the Tchebyshev polynomial coefficients, $a_{m,n}$ in a small 2-dimensional table. You will never need to get to polynomial orders of more than 9 or 10. So on initialization, you can compute the table of $a_{m,n}$. From that table and from your set of harmonic amplitudes $b_m$, you can compute the net polynomial coefficients $a_n$. Jan 2 '21 at 9:07
• Robert just FYI I sent an email. Jan 3 '21 at 18:57