Can harmonics be created or simulated a la carte

Question

I am trying to learn about adding harmonics (audio) with a goal to make a plugin that can allow adjustment of the amplitude of each harmonic being added. I would not have the ability to write this but wanted to get an opinion about whether it is possible. Example: to add a 2nd harmonic and adjust its level separately, add a 3rd and adjust its level etc. Ideally have a GUI that allows the ability to add harmonics and turn a knob or slider to adjust the amplitude of each harmonic up to around the 10th. The goal is to find out from some someone if it is possible and if so I’d pursue someone to write this for me.

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

Answer 1

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.

Answer 2

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