Harmonic distortion of sine signal. plotting resulting spectrum

Question

I'm trying to learn about how the spectrum of a sine signal changes when fed into a waveshaper and i've learnt a great deal from you guys here on stackexchange. Specifically how a transfer function such as $y = x + x^7$ adds the fundamental frequency of $x$ and then odd overtones up to the $7^{th}$ overtone.

Since i'm no mathematician and just now learning some dsp for an effect plugin i'm developing i need to be able to plot the spectrum of $y$ to learn more about the behavior how ploynomials affect this spectrum.

I could use Fourier tranform in matlab on $y$, but it feels like an unnecessary step.

I would like to be able to plot it like in this stackexchang question: Add odd/even harmonics to signal?

Thanks in advance.

Answer 1

Here is a way that you can approach this mathematically. It's not the only way. Another way to approach this is with Tchebyshev polynomials and I will leave that to someone else.

This approach requires understanding two ingredients.

1. an identity that is a consequence of Euler's formula:

$$ \cos(\omega t) = \tfrac12(e^{j\omega t}+e^{-j\omega t}) $$ or $$ \sin(\omega t) = \tfrac1{2j}(e^{j\omega t}-e^{-j\omega t}) $$

2. the binomial theorem:

$$\begin{align} (a+b)^N &= \sum_{n=0}^N {N \choose n} \, a^{N-n} \, b^n \\ &= \sum_{n=0}^N {N \choose n} \, a^n \, b^{N-n} \\ \end{align}$$ or $$\begin{align} (a+b)^N &= \sum_{n=0}^N \frac{N!}{(N-n)! \, n!} \, a^{N-n} \, b^n \\ &= \sum_{n=0}^N \frac{N!}{(N-n)! \, n!} \, a^n \, b^{N-n} \\ \end{align}$$

Now you might see where this is going. So using the cosine identity (so i don't have to mess around with powers of $j$ in the sine identity):

$$\begin{align} \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\ \\ &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, (e^{j\omega t})^{k-n} \, (e^{-j\omega t})^n \\ \\ &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-n)\omega t} \, e^{-jn\omega t} \\ \\ &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\ \end{align}$$

It appears that we will have to deal with the two cases of even $k$ vs. odd $k$ differently.

---

For even $k \ge 2$:

$$\begin{align} \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\ \\ &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\ \\ &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \frac{k!}{\tfrac{k}{2}! \, \tfrac{k}{2}!} + \sum_{n=\tfrac{k}{2}+1}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \right) \\ \\ &= \frac{1}{2^k} \left(\frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{n! \, (k-n)!} \, e^{j(k-2(k-n))\omega t} \right) \\ \\ &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{1}{2^k} \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, \Big( e^{j(k-2n)\omega t} + e^{-j(k-2n)\omega t} \Big) \\ \\ &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{1}{2^k} \sum_{n=0}^{\tfrac{k}{2}-1} \frac{k!}{(k-n)! \, n!} \, 2 \cos\big((k-2n)\omega t\big) \\ \\ &= \frac{1}{2^k} \frac{k!}{\left(\tfrac{k}{2}!\right)^2} + \frac{2 k!}{2^k} \sum_{n=1}^{\tfrac{k}{2}} \frac{1}{\left(\tfrac{k}{2}+n\right)! \, \left(\tfrac{k}{2}-n\right)!} \, \cos(2n\omega t) \\ \end{align}$$

Note that only DC and even-numbered harmonics come outa the even powers of $\cos^k(\omega t)$.

---

For odd $k \ge 1$:

$$\begin{align} \cos^k(\omega t) &= \frac{1}{2^k}(e^{j\omega t}+e^{-j\omega t})^k \\ \\ &= \frac{1}{2^k} \sum_{n=0}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \\ \\ &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=\tfrac{k+1}{2}}^k \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} \right) \\ \\ &= \frac{1}{2^k} \left( \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, e^{j(k-2n)\omega t} + \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{n! \, (k-n)!} \, e^{j(k-2(k-n))\omega t} \right) \\ \\ &= \frac{1}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, \Big( e^{j(k-2n)\omega t} + e^{-j(k-2n)\omega t} \Big) \\ \\ &= \frac{1}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{k!}{(k-n)! \, n!} \, 2 \cos\big((k-2n)\omega t\big) \\ \\ &= \frac{2 k!}{2^k} \sum_{n=0}^{\tfrac{k-1}{2}} \frac{1}{\left(\tfrac{k+1}{2}+n\right)! \, \left(\tfrac{k-1}{2}-n\right)!} \, \cos\big((2n+1)\omega t\big) \\ \end{align}$$

Note that only odd-numbered harmonics come outa the odd powers of $\cos^k(\omega t)$.

---

So now what if we pass this sinusoid:

$$ x(t) = A \, \cos(\omega t) $$

through this polynomial mapping function:

$$\begin{align} y(t) &= a_0 + a_1 x(t) + a_2 x^2(t) + a_3 x^3(t) + ... + a_{K-1} x^{K-1}(t) + a_K x^K(t) \\ &= \sum_{k=0}^{K} a_k \big( x(t) \big)^k \\ \end{align} $$

?

I believe we have to split the even-order terms from the odd-order terms in the summation

$$\begin{align} y(t) &= \sum_{k=0}^{K} a_k \big( x(t) \big)^k \\ \\ &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \big( x(t) \big)^{2k} + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \big( x(t) \big)^{2k+1} \\ \\ &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \big( A \cos(\omega t) \big)^{2k} + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \big( A \cos(\omega t) \big)^{2k+1} \\ \\ &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} A^{2k} \cos^{2k}(\omega t) + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} A^{2k+1} \cos^{2k+1}(\omega t) \\ \\ &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} A^{2k} \left( \frac{1}{2^{2k}} \frac{(2k)!}{(k!)^2} + \frac{2 (2k)!}{2^{2k}} \sum_{n=1}^{k} \frac{1}{\left(k+n\right)! \, \left(k-n\right)!} \, \cos(2n\omega t) \right) + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} A^{2k+1} \left( \frac{2 (2k+1)!}{2^{2k+1}} \sum_{n=0}^{k} \frac{1}{(k+1+n)! \, (k-n)!} \, \cos\big((2n+1)\omega t\big) \right) \\ \\ &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{A^{2k}}{2^{2k}} \frac{(2k)!}{(k!)^2} + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} \sum_{n=1}^{k} a_{2k} \frac{2 A^{2k} (2k)!}{2^{2k}\left(k+n\right)! \left(k-n\right)!} \, \cos(2n\omega t) \qquad + \sum_{k=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} \sum_{n=0}^{k} a_{2k+1} \frac{2 A^{2k+1} (2k+1)!}{2^{2k+1}(k+1+n)!(k-n)!} \, \cos\big((2n+1)\omega t\big) \\ \\ &= a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{A^{2k}}{2^{2k}} \frac{(2k)!}{(k!)^2} + \sum_{n=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} \sum_{k=n}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{2 A^{2k} (2k)!}{2^{2k}\left(k+n\right)! \left(k-n\right)!} \, \cos(2n\omega t) \qquad + \sum_{n=0}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} \sum_{k=n}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \frac{2 A^{2k+1} (2k+1)!}{2^{2k+1}(k+1+n)!(k-n)!} \, \cos\big((2n+1)\omega t\big) \\ \end{align} $$

$\lfloor u \rfloor$ is the floor() function, which is the largest integer no larger than the argument $u$.

$$ \lfloor u \rfloor \le u < \lfloor u \rfloor + 1 $$ or $$ u-1 < \lfloor u \rfloor \le u $$

Also note that if $K$ is even $$ \left\lfloor \tfrac{K-1}{2} \right\rfloor < \left\lfloor \tfrac{K}{2} \right\rfloor = \tfrac{K}{2} \quad \text{and} \quad \left\lfloor \tfrac{K-1}{2} \right\rfloor < \tfrac{K-1}{2} $$

If $K$ is odd $$ \left\lfloor \tfrac{K}{2} \right\rfloor = \left\lfloor \tfrac{K-1}{2} \right\rfloor = \tfrac{K-1}{2} \quad \text{and} \quad \left\lfloor \tfrac{K}{2} \right\rfloor < \tfrac{K}{2} $$

---

$$\begin{align} y(t) &= \sum_{k=0}^{K} a_k \big( A \cos(\omega t) \big)^k \\ \\ &= \sum_{n=0}^{K} b_n \cos(n \omega t) \\ \end{align}$$

where, apparently the DC component of the Fourier cosine series is:

$$ b_0 = a_0 + \sum_{k=1}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{(k!)^2} \left(\frac{A}{2}\right)^{2k} $$

The coefficient for the even harmonic terms:

$$ b_{2n} = 2 \sum_{k=n}^{\left\lfloor \tfrac{K}{2} \right\rfloor} a_{2k} \frac{(2k)!}{\left(k+n\right)! \left(k-n\right)!} \left(\frac{A}{2}\right)^{2k} $$

And the coefficient for the odd harmonic terms:

$$ b_{2n+1} = 2 \sum_{k=n}^{\left\lfloor \tfrac{K-1}{2} \right\rfloor} a_{2k+1} \frac{(2k+1)!}{(k+1+n)!(k-n)!} \left(\frac{A}{2}\right)^{2k+1} $$

Answer 2

FFT does exactly what you are after. For example, in Python:

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as fft
from __future__ import division

# Number of samples
N = 500
# Sampling frequency
fs = 1000
# Fundamental frequency
f0 = 50

t = np.linspace(0, N/fs, N)
tone = np.sin(f0 * 2*np.pi*t)

# The signal after waveshaping
y = tone + 3 * tone**3 + 2 * tone**7

yf = fft.fft(y)
xf = np.linspace(0.0, 1/(2/fs), N/2)

fig, ax = plt.subplots()
ax.plot(xf, np.abs(yf[:N//2]))
plt.show()

That script is using the polynomial $x+3x^3+2x^7$ (arbitrarily).

Answer 3

I'ld like to provide an alternative frequency domain derivation for the discrete-time spectrum of memoryless polynomial nonlinearities for inputs of the form $x[n] = A \cos(\omega_0 n)$, which would complement the time domain derivation provided by RBJ. As revelaed, it's longer to describe but shorter to compute.

Assume a sufficiently oversampled sinusoidal input signal : $$x[n] = A \cos(\omega_0 n) ~~~,~~~\text{ for } -\infty < n < \infty \tag{1}$$

We are interested in the spectral representation of the signal $y_m[n]$ defined with the following nonlinearity : $$ \boxed{ y_m[n] =\left( x[n] \right)^m } \tag{2}$$

Since, multiplication in the time domain transforms into convolution in the frequency domain

$$f[n] g[n] \longleftrightarrow \frac{1}{2\pi} F(\omega) \star G(\omega)\tag{3} $$ we have the following relationship between the DTFT of $y_k[n]$ and the DTFT of $x[n]$ , which is induced for $k=1,2,...,m$

\begin{align} y_1[n] = x[n] \implies Y_1(\omega) &= X(\omega) \\ y_2[n] = x^2[n] \implies Y_2(\omega) &= \left( \frac{1}{2\pi} \right) X(\omega) \star X(\omega) \\ y_3[n] = x^3[n] \implies Y_3(\omega) &= \left( \frac{1}{2\pi} \right)^2 X(\omega) \star X(\omega) \star X(\omega) \\ ... \\ y_m[n] = x^m[n] \implies Y_m(\omega) &= \left( \frac{1}{2\pi} \right)^{m-1} \underbrace{ X(\omega) \star ... \star X(\omega)}_{m-1 \text{ times}} \tag{4} \\ \end{align}

where $Y_k(\omega) = DTFT\{y_k[n]\}$ as expected.

Above derivation simply suggests the following recursive relation $$ \boxed{ Y_{m+1}(\omega) = \frac{1}{2\pi} Y_m(\omega) \star X(\omega) } \tag{5}$$

Furthermore, recognising the DTFT of $x[n]$ as $$X(\omega) = A \pi \left[~~~ \delta(\omega + \omega_0) ~~ + ~~~ \delta(\omega - \omega_0) ~~~\right] \tag{6}$$

The recursion for $Y_{m+1}(\omega)$ becomes like this: first left shift $Y_m(\omega)$ by $\omega_0$, then right shift $Y_m(\omega)$ by $\omega_0$, and add the two (and multiply by $ {A\pi}/{2\pi} =A/2$).

Based on this procedure, the following figure-1 is obtained for the spectrum of $y_m[n] = x^m[n] = \left(A \cos(\omega_0 n) \right)^m$ for the values $m=1,2,3$.

In figure-1, the impulse weights are shown as integers, omitting the $(A/2)^m$ weight for clarity, so that the binomial nature of them are revealed. By a visual inspection of this figure, one can arrive at the following closed form expression for $Y_m(\omega)$ :

$$\boxed{ Y_m(\omega) = 2\pi ( \frac{A}{2} )^m \sum_{k=0}^{m} {m \choose k} \delta(\omega + m \omega_0 - 2\omega_0 k) } \tag{7}$$ where $k=0$ indicates the leftmost impulse at $\omega = -m \omega_0$ , and $k=m$ indicates the rightmost one at $\omega = m \omega_0$, and between two consecutive impulses is a distance of $2 \omega_0$.

This closes the first part of the problem, which shows the closed form spectrum expression for a single term of the form $y[n] = x^m[n]$ when $x[n]$ is $A \cos(\omega_0 n)$.

Extending it to include the general polynomial sum of order $m$, as did RBJ, we have : $$ y[n] = P_K(x[n]) = a_0 x^0[n] + a_1 x^1[n] + a_2 x^2[n] + ... + a_K x^K[n] = \sum_{r=0}^{K} a_r x^r[n] \tag{8}$$

Indeed, $Y(\omega)$ can trivially be extended by utilizing the linearity of the DTFT. Thus we can immediately conclude the result as:

$$\boxed{ Y(\omega) = \sum_{r=0}^{K} a_r Y_r(\omega) = \sum_{r=0}^{K} a_r \left( 2\pi ( \frac{A}{2} )^r \sum_{k=0}^{r} {r \choose k} \delta(\omega + r \omega_0 - 2\omega_0 k) \right) }\tag{9}$$

Yet, however, one also wishes to explicitly compute the impulse weights $B_k$ for each impulse $B_k \delta(\omega - k\omega_0)$ located at the frequency $\omega = k \omega_0$ so that the spectrum $Y(\omega)$ is expressed as

$$ Y(\omega) = \sum_{k=-K}^{K} B_k \delta(\omega - k \omega_0) \tag{10} $$

To do so requires to divide the indices $k$ into even and odd values and to perform two independent summations as did RBJ in his answer. However, I would like to make a simple trick and compute the weights in a single formula.

For this purpose first, I have modified the spectrum $Y_m(\omega)$ of $y_m[n]$ as in the figure-2:

In the figure-2, the black-or-white? square impulses are actually nonexistant, hence their weights are $0$. Yet with this trick we assume that the spectrum $Y_m(\omega)$ includes impulses at every $k\omega_0$ stop, rather than at $2\omega_0 k$ as previosuly was assumed. We have to change the formula of $Y_m(\omega)$ a bit to reflect this modification:

$$\boxed{ Y_m(\omega) = 2\pi (\frac{A}{2} )^m \sum_{k=0}^{2m+1} \begin{bmatrix} m \\ k/2 \end{bmatrix} \delta(\omega + m \omega_0 - k\omega_0) }\tag{11}$$

where $$ \begin{bmatrix} m \\ k/2 \end{bmatrix} = \begin{cases} {m \choose k/2} ~~~&, k \text{ even } \\ 0 ~~~&, k \text{ odd} \end{cases} \tag{12}$$

can be called as the conditional combinatorial operator.

Now, since the coefficients $B_k$ are conjugate-symmetric for real inputs; i.e., $$ B_{-k} = B_{k}^* $$ we shall evaluate them only for one half of the figure-2, which I chose to be the positive half; i.e., for $k = 0,1,..,K$.

The formula for $Y_m(\omega)$ which includes only the positive frequency impulses is defined as (not to be confused with analytic signal)

$$ \begin{align} Y_m^+(\omega) &= 2\pi (\frac{A}{2} )^m \sum_{k=m}^{2m+1} \begin{bmatrix} m \\ k/2 \end{bmatrix} \delta(\omega + m \omega_0 - k\omega_0) \\ &= 2\pi (\frac{A}{2} )^m \sum_{k=0}^{m} \begin{bmatrix} m \\ (k+m)/2 \end{bmatrix} \delta(\omega - k\omega_0) \tag{13}\\ \end{align} $$

At this point, the following figure-3 reveals the range of summations along sub-coefficients of each $Y_m^+(\omega)$ for each final coefficient $B_k$

From figure-3 (note that green line index $L$ is replaced with $r$ in the below formulas) it can be seen that for computing each final coefficient $B_k$ for the impulse located at $\delta(\omega - k \omega_0)$ the following sum is needed:

$$\boxed{ B_k = 2\pi \sum_{r=k}^K a_r (A/2)^r \begin{bmatrix} r \\ (k+r)/2 \end{bmatrix} }\tag{14}$$

Then, finally the spectrum $Y(\omega)$ for the polynomial memoryless nonlinear mapping $$y[n] = \sum_{r=0}^{K} a_r x^r[n] \tag{15}$$

in terms of final coefficients $B_k$ for each impulse, is found to be:

$$ Y(\omega) = \sum_{k=-K}^{K} B_k \delta(\omega - k \omega_0) \tag{16} $$

which becomes

$$ \boxed{ Y(\omega) = 2\pi \sum_{k=-K}^{K} \left( \sum_{r=|k|}^K a_r (A/2)^r \begin{bmatrix} r \\ \frac{|k|+r}{2} \end{bmatrix} \right) \delta(\omega - k \omega_0) } \tag{17} $$

This closes the second part of the problem, assuming there are no errors in the derivations.

NOTE: The answer provided by RBJ computes the weights of the cosine terms, denoted as $b_{2n}$ and $b_{2n+1}$, for even and odd terms separation. Here, I provided frequency domain solution corresponding to the coefficients of frequency domain impulses, located at $\omega = k \omega_0$, which is denoted as $B_k$ (for both even and odd terms), and which has the following relation to those time domain coefficients $b_n$ as:

$$ \boxed{ \frac{ B_{k} }{\pi} = b_k ~~~,~~~ \text{ for } -K \leq k \leq K } \tag{18}$$