# Harmonic distortion of sine signal. plotting resulting spectrum

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?

• Welcome to DSP.SE! Why would the Fourier transform be unnecessary? It provides exactly what you are looking for. – Tendero Sep 5 '18 at 15:23
• Just some terminology FYI, the input-output relationship you have described $y = f(x)$ is technically a transfer characteristic of a memoryless nonlinear (MNL) system, whereas the term "transfer function" in the context of signal processing usually refers to the ratio of the Laplace transforms of the output divided by the input, i.e. $H(s) = Y(s)/X(s)$, of a linear time invariant (LTI) system. – Robert L. Sep 5 '18 at 22:05
• Hi. Thanks for clearing that up. There's a lot to learn indeed. – antwoord Sep 8 '18 at 12:13

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

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

• note that i am "borrowing" two graphics from that other answer referred to in the question. – robert bristow-johnson Sep 6 '18 at 0:02
• and will someone please check my math. i don't wanna play too fast-and-loose with the limits in the summations. – robert bristow-johnson Sep 6 '18 at 0:03
• ok. The thundercloud and below seems consistent. Assuming the formula for $cos^k(\omega t)$ from the first part was also right, then it seems ok... – Fat32 Sep 7 '18 at 10:02
• I will look into this. Seems to be a good idea to try to learn this. I will look into Euler's formula and the Binomial theorem. It's a lot to delve in to but now i know where to start. Thanks. – antwoord Sep 8 '18 at 12:19
• Binomial Theorem is college-level algebra (or advanced high-school algebra). Euler's formula requires a little bit of calculus and (of course) complex variables math (what you also get in calculus) and is fundamental for electrical engineering (or mechanical engineering) to understand and model sinusoidal analysis. it's where we originally get the notion of negative frequencies from. – robert bristow-johnson Sep 9 '18 at 2:52

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

• I tried it in Matlab and i got the correct result. Thanks – antwoord Sep 8 '18 at 13:37
• @antwoord No problem! If this solved your problem, please mark it as accepted so it can help future readers. – Tendero Sep 8 '18 at 13:39

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

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

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)$$ 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,3,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}} \\ \end{align}

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

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

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

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 the figure-1, the impulse weights are shown as integers, omitting the $(A/2)^m$ weighting, so that the binomial nature of them are revealed. By the 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) }$$ 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]$$

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

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

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 previsouly did. 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) }$$

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

is the conditional combinatorial operator (as I named it so!).

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)\\ \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 the 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} }$$

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

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

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

And 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 by $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 }$$

• the derivations were s'posed to be straight-forward. but they became a mess. i think my only screwup is in the swapping of summations in the final polynomial summation. i believe my expressions for even and odd power monomials are still correct. – robert bristow-johnson Sep 6 '18 at 23:07
• i think i fixed the mistake. i realize that the expressions in the thundercloud look formidable, but i tried not to skip any steps. it starts out simple enough and each step should be able to be clearly justified from the previous step. – robert bristow-johnson Sep 7 '18 at 2:02
• @robertbristow-johnson Hi! I finally harmonized them as much as I could by replacing [m , b, c] of mine with [ K , a , b] of yours. But the final b of yours is B in my derivation, for the sake of differentiating frequency from time domains. They yield the same results, yet there is a difficulty in running those $b_{2n}$ and $b_{2n+1}$, what's the range of $n$ there ? – Fat32 Sep 13 '18 at 0:22
• maybe i should modify mine. i gotta go through yours and pick out the corresponding $e^{j n \omega t}$ in my answer to the $\delta(f - \tfrac{n \omega}{2\pi})$ in your answer. they should be the same. – robert bristow-johnson Sep 13 '18 at 0:27
• @robertbristow-johnson you will only have to multiply each $b_n$ by $\pi$ so that it will be my $B_n$... It's up to you. They are quite ok, seems to me. But again up to you...I have to go now, may look at it later. – Fat32 Sep 13 '18 at 0:30