Fourier series

The last chapter ended in the hush before a concert. Sound is the air pressure at the eardrum changing over time, that is, a function of just one variable, $p(t)$. A violin and a flute play the same note A, and for both the pressure curve repeats 440 times a second. But the curves themselves are unlike: the flute’s wave is almost smooth, the violin’s angular, like the teeth of a saw. The ear doesn’t confuse them. It breaks sound down into pure tones, sine waves, and hears which of them are present and how loud they are.

In this chapter we’ll open a recording studio and do with a signal what a sound engineer does: listen, split it into tracks, turn the knobs and mix it back together. There are two questions. Can any periodic signal be built from sine waves? And how can you tell from a finished curve which ones it’s made of?

One note, different voices

The simplest sound is a pure tone: the pressure varies as a sine wave, $p(t) = A\sin(2\pi\nu t)$. The number $\nu$ is the frequency, how many oscillations fit into a second; $A$ is the amplitude, on which the loudness depends. A tuning fork gives an almost pure tone. For the A above middle C $\nu = 440$ Hz, and one oscillation lasts $\frac{1}{440}$ of a second, about $2.27$ milliseconds.

Add a sine wave with frequency $880$ Hz. During one oscillation of the fundamental it manages exactly two, so the sum again repeats after $\frac{1}{440}$ of a second. The pitch of the note is the same; only the waveform changes. The frequencies $3 \cdot 440$, $4 \cdot 440$ and so on behave the same way: they all fit into the period of the fundamental a whole number of times. Musicians call them overtones, mathematicians harmonics.

From now on it’s convenient to measure time so that the period of the fundamental becomes $2\pi$: instead of $t$ we write $x = 2\pi\nu t$. Then harmonic number $n$ is the sine wave $A\sin(nx + \varphi)$, and by the formula for the sine of a sum it equals $A\sin\varphi\cos nx + A\cos\varphi\sin nx$. Any harmonic can be written as a mix of a cosine and a sine of the same frequency, and the shift $\varphi$, the phase, is hidden in the proportions of the mix.

The constant component: the average value of the signal over a period. Why it’s convenient to write it with a half will become clear once the coefficient formulas appear. The cosine part of the $n$th harmonic. The sine part of the same harmonic. The harmonic’s amplitude is $\sqrt{a_n^2 + b_n^2}$, and its phase is determined by the ratio $a_n : b_n$. Example: $f(x) = \sin x + \frac12\sin 2x + \frac13\sin 3x$. Here $a_0 = 0$, all $a_n = 0$, $b_1 = 1$, $b_2 = \frac12$, $b_3 = \frac13$. These are the first three harmonics of the “sawtooth” we started with: if the series $\frac{\sin nx}{n}$ is continued, the teeth get sharper.

Shuffle the phases a few times. The picture is new every time, yet to the ear the sound hardly changes. It seems the ear perceives the set of amplitudes rather than the waveform itself. This guess was put forward in 1843 by Georg Ohm, the same Ohm who discovered the law for electric circuits, and twenty years later Hermann Helmholtz built a theory of hearing on it. The rule isn’t absolute: in some sounds phases can be told apart after all. But it serves the sound engineer well, and the timbre of an instrument is described by the set of amplitudes of its harmonics.

A signal is the sum of two pure tones, 440 Hz and 660 Hz. After how long does it repeat?

After $\frac{1}{220}$ of a second After $\frac{1}{440}$ of a second After $\frac{1}{660}$ of a second Never: the frequencies are different

$440 = 2 \cdot 220$ and $660 = 3 \cdot 220$, so both tones are the second and third harmonics of the frequency $220$ Hz. In $\frac{1}{220}$ of a second the first makes exactly two oscillations and the second three, and the sum repeats. The fundamental frequency of $220$ Hz isn’t in the signal at all.

Where it all began: heat in a ring

The first to break functions down into sine waves weren’t musicians but a physicist who studied heat. Joseph Fourier went to Egypt with Napoleon, then became prefect of the Isère department in Grenoble and, between administrative duties, worked on heat conduction. In December 1807 he presented a memoir on the propagation of heat to the Paris Academy. The committee, which included Lagrange and Laplace, received it coolly, and the memoir wasn’t printed. Later the Academy made the propagation of heat the subject of its 1811 prize competition, and Fourier won it, although the report again reproached him for lack of rigour. The full theory came out in 1822 as the book The Analytical Theory of Heat.

Take a thin wire ring and heat part of it unevenly. We’ll give the position of a point on the ring by the angle $x$, and its temperature at time $t$ by the function $u(x, t)$. Heat flows from hot to cold, and a point warms up the faster, the warmer its neighbours are on average. This “difference from the neighbours” is measured by the second derivative in $x$. If units are chosen so that the coefficient is one, we get the heat equation

$$\frac{\partial u}{\partial t} = \frac{\partial^2 u}{\partial x^2}.$$

We know partial derivatives from the last chapter. Fourier noticed that in this equation harmonics live separately. The function $u = e^{-n^2 t}\sin nx$ satisfies it: the derivative in $t$ gives a factor $-n^2$, and the second derivative in $x$ gives the same factor $-n^2$. The same is true of $e^{-n^2 t}\cos nx$. Each harmonic stays a harmonic and simply fades at its own rate. The equation is linear, so a sum of solutions is a solution too. If the initial temperature can be broken down into harmonics, the solution is ready:

The damping factor. The higher the harmonic’s number, the faster it fades: by the time $t = 0.1$ the first keeps $e^{-0.1} \approx 90\,\%$ of its amplitude, the third $e^{-0.9} \approx 41\,\%$, the fifth $e^{-2.5} \approx 8\,\%$, and the tenth five thousandths of a per cent. Harmonic number $n$ of the initial temperature $u(x, 0)$. The constant $\frac{a_0}{2}$ is the average temperature of the ring. It doesn’t change: the ring gives no heat to the outside, heat only flows around it. The formula holds if the initial profile can be expanded in a series; for $t > 0$ the factors $e^{-n^2t}$ fall so fast that the series can be differentiated term by term.

The problem has come down to a single question: how do you break the initial temperature down into harmonics? And it can be anything at all, for example half the ring hot and half cold. So Fourier had to claim that even a function with a step can be broken down into sine waves.

A step from sine waves

That’s exactly what outraged Lagrange. Sine waves are smooth, and their sums are smooth too: where would a step come from? Let’s test it on the most awkward signal, the square wave: $\operatorname{sq}(x) = 1$ for $0 < x < \pi$ and $\operatorname{sq}(x) = -1$ for $-\pi < x < 0$, repeated with period $2\pi$ from there on. Fourier’s answer:

A common factor. Where it comes from, the coefficient formulas will show. The fundamental, with amplitude $\frac{4}{\pi} \approx 1.27$: it’s even taller than the step. The third harmonic. It presses down the hump in the middle and makes the edges steeper. All the harmonics are odd. For the square wave the second half of the period is the first turned upside down: $\operatorname{sq}(x + \pi) = -\operatorname{sq}(x)$. When shifted by $\pi$ the sine wave $\sin nx$ is multiplied by $(-1)^n$, so only the odd harmonics have this property, and the even ones have no business being in a square wave. Example: at $x = \frac{\pi}{2}$ the sines are $1, -1, 1, -1, \dots$, and the formula promises $\frac{4}{\pi}\bigl(1 - \frac13 + \frac15 - \frac17 + \dots\bigr) = 1$, that is, $1 - \frac13 + \frac15 - \dots = \frac{\pi}{4}$. That’s the slow Leibniz series from the chapter on series. The partial sums close in on one unhurriedly: $1.273$, $0.849$, $1.103$, …, and after ten harmonics $0.968$.

Each partial sum is a continuous function, and at the point $x = 0$ all the sines are zero. So the series converges to zero there, exactly the middle of the jump. The step appears only in the limit. Lagrange was wrong: a limit of continuous functions can be discontinuous. But he wasn’t worried for nothing: such a limit has to be handled with care. We’ll find out exactly where the series converges and to what. But first let’s work out where the coefficients $\frac{4}{\pi}$, $\frac{4}{3\pi}$, $\frac{4}{5\pi}$ come from.

Why the tracks don’t interfere

We’ve already solved a similar problem for arrows. In the chapter on vectors the projection of a vector onto a unit axis was a dot product. If the axes $\vec e_1$ and $\vec e_2$ are perpendicular and $\vec a = a_1\vec e_1 + a_2\vec e_2$, then $\vec a \cdot \vec e_1 = a_1$: the second term drops out because $\vec e_2 \cdot \vec e_1 = 0$. A coordinate can be found without solving any equations: you only need to take the dot product with the axis.

Let’s do the same with functions. A vector has a few coordinates; a function has as many “coordinates” as there are points $x$, and the sum of products of coordinates turns into the integral of a product.

The first function, for example a signal. The second function, for example the harmonic we’re projecting onto. Example: $\langle \sin x, \sin x\rangle = \int_{-\pi}^{\pi}\sin^2 x\,dx = \pi$, $\langle 1, 1\rangle = 2\pi$, $\langle x, \sin x\rangle = \bigl[\sin x - x\cos x\bigr]_{-\pi}^{\pi} = 2\pi$. This product has the same properties as the dot product of vectors: it’s symmetric, linear in each argument, and $\langle f, f\rangle = \int f^2\,dx \ge 0$ plays the role of the squared length.

The trick is to turn a product of sine waves into a sum: the areas under the terms can be computed separately.

Take two sine waves, $\p1{\sin mx}$ and $\p2{\sin nx}$; in the picture $m = 3$ and $n = 2$, and the crests can be dragged to change the frequencies. We need the signed area under their product on the interval $[-\pi, \pi]$. Subtract the formulas for the cosine of a difference and the cosine of a sum: $\cos(\alpha - \beta) - \cos(\alpha + \beta) = 2\sin\alpha\sin\beta$. With $\alpha = mx$, $\beta = nx$ we get $\p3{\sin mx\sin nx} = \p4{\tfrac12\cos(m - n)x} - \p5{\tfrac12\cos(m + n)x}$. In the picture the difference of the two cosine waves coincides exactly with the product. The area under $\cos kx$ for an integer $k \ne 0$ on an interval of length $2\pi$ is zero. The antiderivative $\frac{\sin kx}{k}$ is zero at both ends, because $\sin(\pm k\pi) = 0$. In the picture you can see it like this: exactly $|k|$ periods fit on the interval, and each hump above the axis is balanced by an identical dip. If $m \ne n$, both cosine waves have non-zero frequencies $m - n$ and $m + n$, and the whole integral is zero. If $m = n \ge 1$, the first cosine wave degenerates into the constant $\frac12$, the area under it is $\frac12 \cdot 2\pi = \pi$, and the second still gives zero. Cosines are dealt with in the same way through $\cos mx\cos nx = \frac12\cos(m - n)x + \frac12\cos(m + n)x$. Finally, $\sin mx\cos nx = \frac12\sin(m + n)x + \frac12\sin(m - n)x$, and the area under a sine on an interval symmetric about zero is zero: the sine is odd, and its hump on the right is balanced by a dip on the left.

The same argument about areas gives a useful rule about symmetry. A signal is called even if $f(-x) = f(x)$, like $|x|$, and odd if $f(-x) = -f(x)$, like the square wave.

The idea: the product of an odd function and an even one is odd, and for an odd function the areas to the left and right of zero cancel out.

That $\p1{f}$ is odd means its graph goes to itself under a $180^\circ$ rotation about the origin: the left half is the right half turned through half a turn. The cosine $\p2{\cos nx}$ is even: its left half is the mirror image of the right half in the vertical axis, because $\cos(-nx) = \cos nx$. The product $\p3{f(x)\cos nx}$ is odd: $f(-x)\cos(-nx) = -f(x)\cos nx$. Each of its points on the right has a matching point on the left with the opposite sign. So the signed area under the product on $[-\pi, 0]$ equals the area on $[0, \pi]$ with the opposite sign, and together they give zero: $\langle f, \cos nx\rangle = 0$. For $n = 0$ this means $\int f\,dx = 0$. For an even $f$ it’s all the same: the product $f(x)\sin nx$ is odd, and $\langle f, \sin nx\rangle = 0$.

Splitting into tracks

Now the coefficients are found the same way as the coordinates of a vector. Let $f = \frac{a_0}{2} + \sum_k (a_k\cos kx + b_k\sin kx)$. Take the dot product with $\cos nx$: all the terms but one are orthogonal to this “probe” and drop out, and the remaining one gives $a_n\langle \cos nx, \cos nx\rangle = \pi a_n$. For $n = 0$ the probe is one: $\langle f, 1\rangle = \frac{a_0}{2}\cdot 2\pi = \pi a_0$. That’s why the constant is written with a half: the formula for $a_0$ becomes the same as for the others.

A question of honesty remains: can an infinite sum be multiplied by $\cos nx$ and integrated term by term, like a finite one? Not always. This needs uniform convergence: the remainder of the series after the $N$th term must be small on the whole interval at once, not at each point separately.

We’ll find a coordinate as for a vector: take the dot product with the axis. Uniformity will do the work of swapping the sum and the integral.

The function $f$ is a uniform limit of continuous partial sums, so it’s continuous (this is proved in the chapter on rigorous analysis) and the integrals in the formulas exist. Fix $n$ and take $\varepsilon > 0$. For $N \ge \max(n, N_0)$ we have $|f(x) - S_N(x)| < \varepsilon$ on the whole interval, and $|\cos nx| \le 1$. So $\left|\int_{-\pi}^{\pi}f\cos nx\,dx - \int_{-\pi}^{\pi}S_N\cos nx\,dx\right| \le \int_{-\pi}^{\pi}\varepsilon\,dx = 2\pi\varepsilon$. The sum $S_N$ is finite and can be integrated term by term. By the orthogonality theorem one term survives: $\int S_N\cos nx\,dx = a_n \cdot \pi$ for $n \ge 1$ and $\frac{a_0}{2}\cdot 2\pi = \pi a_0$ for $n = 0$. We’ve got $\left|\int_{-\pi}^{\pi}f\cos nx\,dx - \pi a_n\right| \le 2\pi\varepsilon$ for every $\varepsilon > 0$. A non-negative number that’s smaller than every positive one is zero, so $\pi a_n = \int f\cos nx\,dx$. For $b_n$ the argument is the same with the factor $\sin nx$.

The formulas can be applied to any integrable function, even without knowing in advance whether anything converges. That’s how the main object of the chapter is defined.

The signal over one period. The probe for the cosine part of the $n$th harmonic; at $n = 0$ it’s one, and $a_0$ is twice the average value of the signal. The probe for the sine part. We divide by $\pi$ because $\langle \sin nx, \sin nx\rangle = \pi$: that’s the squared “length” of the axis. Example: for the square wave $a_n = 0$, because it’s odd, and $b_n = \frac{1}{\pi}\Bigl(\int_{-\pi}^{0}(-\sin nx)\,dx + \int_0^{\pi}\sin nx\,dx\Bigr) = \frac{2}{\pi}\int_0^{\pi}\sin nx\,dx = \frac{2\bigl(1 - (-1)^n\bigr)}{\pi n}$. For odd $n$ this is $\frac{4}{\pi n}$, for even $n$ zero. We’ve got exactly the Fourier series of the square wave.

Find the coefficient $b_2$ of the function $f(x) = x$ on $(-\pi, \pi)$ (its periodic extension is the “sawtooth”).

By parts: $\int_{-\pi}^{\pi}x\sin nx\,dx = \Bigl[-\frac{x\cos nx}{n}\Bigr]_{-\pi}^{\pi} + \frac1n\int_{-\pi}^{\pi}\cos nx\,dx = -\frac{2\pi\cos n\pi}{n} + 0$. So $b_n = -\frac{2\cos n\pi}{n} = \frac{2(-1)^{n+1}}{n}$, and $b_2 = -1$. The whole series: $x \sim 2\bigl(\sin x - \frac{\sin 2x}{2} + \frac{\sin 3x}{3} - \dots\bigr)$.

Find $a_1$ for $f(x) = |x|$ on $[-\pi, \pi]$. The answer can be written in terms of $\pi$, for example 4/pi.

The function is even, so $a_n = \frac{2}{\pi}\int_0^{\pi}x\cos nx\,dx = \frac{2}{\pi}\Bigl[\frac{x\sin nx}{n} + \frac{\cos nx}{n^2}\Bigr]_0^{\pi} = \frac{2\bigl((-1)^n - 1\bigr)}{\pi n^2}$. At $n = 1$ we get $a_1 = -\frac{4}{\pi} \approx -1.273$. For even $n$ the coefficients are zero, and the odd ones fall off as $\frac{1}{n^2}$, faster than for the square wave: the “triangle” $|x|$ has no jumps, only corners.

The last observation is worth remembering: the smoother the signal, the faster its coefficients fall off. A jump gives $\frac1n$, a corner $\frac{1}{n^2}$. The reason is integration by parts: each time the signal’s smoothness lets you move a derivative onto it, an extra factor $\frac1n$ appears in the coefficient. That’s why in the ring from the earlier section the sharp corners vanished first: they’re held up by high harmonics, and those fade as $e^{-n^2t}$.

You can compute coefficients for piecewise-defined functions, look at partial sums and at the rotating vectors of the harmonics in the Fourier series visualiser, which also has standard textbook problems.

When the series converges to the signal

Any integrable function has coefficients, but that doesn’t yet mean the series converges, let alone to the function itself. Fourier was sure it always did; Lagrange had doubts. The first to honestly prove convergence for a wide class of functions was Peter Gustav Lejeune Dirichlet, in 1829.

The square wave, the “sawtooth”, the “triangle” and any signal drawn with line segments and arcs are piecewise smooth. The proof will need a lemma about fast oscillations: if you multiply a function by a very rapid sine wave, the area under the product almost disappears.

The trick is first to prove the lemma for a step function, where everything is visible, and then to approximate any function by steps.

A single step: $\p1{g = c}$ on the interval $[\alpha, \beta]$. Then $\int_\alpha^\beta c\sin\lambda t\,dt = \frac{c}{\lambda}(\cos\lambda\alpha - \cos\lambda\beta)$, and the absolute value of this number is at most $\frac{2|c|}{\lambda}$. In the picture neighbouring humps and dips of the rapid sine wave cancel each other, and only a sliver at the edges is left. A step function $\p2{s}$ of $m$ steps with heights $c_1, \dots, c_m$: the integrals add, and $\left|\int_a^b s(t)\sin\lambda t\,dt\right| \le \frac{2(|c_1| + \dots + |c_m|)}{\lambda}$. As $\lambda$ grows this tends to zero. Move the diamond on the $\lambda$ scale to the right: the more rapid the sine wave, the smaller the leftover. Now any piecewise continuous $\p3{g}$. On each closed piece it’s continuous, and so uniformly continuous (Cantor’s theorem, chapter on rigorous analysis). Take $\varepsilon > 0$ and split the interval so finely that on each part $g$ changes by less than $\varepsilon$. The step function $s$ equal to $g$ at the left end of each part differs from $g$ by less than $\varepsilon$ everywhere except finitely many points, which don’t affect the integral. Estimate: $\left|\int g\sin\lambda t\right| \le \int|g - s|\cdot|\sin\lambda t| + \left|\int s\sin\lambda t\right| \le \varepsilon(b - a) + \frac{2\sum|c_i|}{\lambda}$. For large $\lambda$ the second term is less than $\varepsilon$, so the integral is less than $\varepsilon(b - a + 1)$ in absolute value. The number $\varepsilon$ was arbitrary, so the limit is zero. With cosine it’s word for word the same.

The idea: write the partial sum as an average of the signal with a weight that gathers near zero as $N$ grows. The contribution of everything away from zero is killed by the Riemann–Lebesgue lemma.

Substitute into the partial sum $S_N(x)$ the coefficient formulas with integration variable $u$. By the formula for the cosine of a difference $a_k\cos kx + b_k\sin kx = \frac{1}{\pi}\int_{-\pi}^{\pi}f(u)\cos k(u - x)\,du$, and after the substitution $u = x + t$ (for a periodic function the integral over any interval of length $2\pi$ is the same) $S_N(x) = \frac{1}{\pi}\int_{-\pi}^{\pi}f(x + t)\,\p1{D_N(t)}\,dt$, where $D_N(t) = \frac12 + \cos t + \cos 2t + \dots + \cos Nt$ is the Dirichlet kernel. The kernel folds up into a fraction. From the product formula $2\sin\frac t2\cos kt = \sin\bigl(k + \frac12\bigr)t - \sin\bigl(k - \frac12\bigr)t$, and summing over $k$ from $1$ to $N$ together with $2\sin\frac t2\cdot\frac12$, everything cancels except $\sin\bigl(N + \frac12\bigr)t$. So $D_N(t) = \frac{\sin(N + \frac12)t}{2\sin\frac t2}$. In the picture it’s a tall narrow peak of height $N + \frac12$ at zero, with fading waves on either side. The area under the kernel: the cosines give zero over a period, leaving $\int_{-\pi}^{\pi}\frac12\,dt = \pi$. The kernel is even, so $\p2{\frac{1}{\pi}\int_0^{\pi}D_N\,dt = \frac12}$, and the same on the left. Write $\p3{h(t) = f(x + t) - f(x + 0)}$. Multiply $f(x + 0)$ by the right half of the area and $f(x - 0)$ by the left half, and subtract from $S_N(x)$. The difference $S_N(x) - \frac{f(x + 0) + f(x - 0)}{2}$ becomes the sum of two integrals: the right one, $\frac{1}{\pi}\int_0^{\pi}h(t)D_N(t)\,dt$, and a similar left one over the interval $[-\pi, 0]$, with $f(x - 0)$ in place of $f(x + 0)$. In the right integral put $\p4{g(t) = \frac{h(t)}{2\sin(t/2)}}$; then it equals $\frac{1}{\pi}\int_0^{\pi}g(t)\sin\bigl(N + \frac12\bigr)t\,dt$. As $t \to 0^+$ the function $g$ doesn’t blow up: $g(t) = \frac{h(t)}{t}\cdot\frac{t}{2\sin(t/2)}$. By Lagrange’s theorem the first factor equals $f'(\xi)$ at a point $\xi$ between $x$ and $x + t$ and tends to the limit of the derivative from the right; the second tends to $1$. So $g$ is piecewise continuous on $[0, \pi]$. By the Riemann–Lebesgue lemma with $\lambda = N + \frac12$ the right integral tends to zero; the left one does the same. So $S_N(x) \to \frac{f(x + 0) + f(x - 0)}{2}$.

Dirichlet’s theorem pays off straight away. Expand $f(x) = x^2$ on $[-\pi, \pi]$. The function is even, so there are no sines; integrating by parts twice, we get $a_0 = \frac{2\pi^2}{3}$ and $a_n = \frac{4(-1)^n}{n^2}$:

$$x^2 = \frac{\pi^2}{3} + 4\sum_{n=1}^{\infty}\frac{(-1)^n\cos nx}{n^2}.$$

The periodic extension of the parabola is continuous, including at the point $x = \pi$, where neighbouring arcs meet at height $\pi^2$. There $\cos n\pi = (-1)^n$, and by Dirichlet’s theorem $\pi^2 = \frac{\pi^2}{3} + 4\sum\frac{1}{n^2}$. Hence $\sum_{n=1}^{\infty}\frac{1}{n^2} = \frac{\pi^2}{6}$: the answer to the Basel problem, which Euler struggled with in the chapter on series, in two lines.

What does the Fourier series of the square wave converge to at the point $x = 0$?

To zero To one The series diverges there To $\frac{4}{\pi}$

All the sines are zero at zero, so each partial sum equals $0$. This agrees with Dirichlet’s theorem: $\frac{f(0 - 0) + f(0 + 0)}{2} = \frac{-1 + 1}{2} = 0$, the middle of the jump.

Ringing at the edge of a step

Look at the partial sums of the square wave near the jump. Each one overshoots the level one, then dives below it and oscillates like a string after it’s plucked. You might think the overshoot disappears as $N$ grows. It doesn’t: it gets ever narrower and hugs the jump, but its height doesn’t drop below $1.179$.

The trick is to stretch the picture at the jump by a factor of $2M$. On the new scale the height of the first peak turns out to be a sum of rectangles under the curve $\frac{\sin t}{t}$.

Put $t = 2Mx$. The picture shows the partial sum on this scale, $\p1{S_M\bigl(\frac{t}{2M}\bigr)}$: for any $M$ its first peak comes at $t = \pi$. Let’s prove it. The derivative is $S_M'(x) = \frac{4}{\pi}\sum_{j=1}^{M}\cos(2j - 1)x$. From $2\sin x\cos(2j - 1)x = \sin 2jx - \sin(2j - 2)x$ the sum telescopes: $\p2{S_M'(x) = \frac{2}{\pi}\cdot\frac{\sin 2Mx}{\sin x}}$. On $\bigl(0, \frac{\pi}{2M}\bigr)$ the derivative is positive, and just after $\frac{\pi}{2M}$ negative, so at $x = \frac{\pi}{2M}$, that is, $t = \pi$, the sum reaches a maximum. The value at the peak: $S_M\bigl(\frac{\pi}{2M}\bigr) = \frac{4}{\pi}\sum_{j=1}^{M}\frac{\sin t_j}{2j - 1}$, where $t_j = \frac{(2j - 1)\pi}{2M}$. Since $\frac{1}{2j - 1} = \frac{\pi}{2Mt_j}$, we get $\frac{2}{\pi}\sum_{j=1}^{M}\p3{\frac{\sin t_j}{t_j}\cdot\frac{\pi}{M}}$. This is the sum of the areas of $M$ rectangles of width $\frac{\pi}{M}$ with heights at the midpoints: a Riemann sum for $\frac{2}{\pi}\int_0^{\pi}\frac{\sin t}{t}\,dt$. The function $\frac{\sin t}{t}$ is continuous on $[0, \pi]$ (at zero it equals $1$), so its Riemann sums converge to the integral, and the height of the peak tends to $\p4{\frac{2}{\pi}\int_0^{\pi}\frac{\sin t}{t}\,dt \approx 1.17898}$. This number doesn’t depend on $M$: the peak gets narrower along with the scale $x = \frac{t}{2M}$, but it doesn’t come down. The overshoot $0.17898$ is $8.949\,\%$ of a jump of height $2$.

There’s no contradiction with Dirichlet’s theorem. It speaks about each point separately: for a fixed $x > 0$ the peak will sooner or later slip to the left of $x$, and after that $S_M(x)$ will calmly approach one. But the convergence isn’t uniform: near the jump there’s always a point where the partial sum is off by almost $0.18$. It couldn’t be otherwise: a uniform limit of continuous functions is continuous, and the square wave isn’t.

The overshoot can be removed if the series isn’t cut off abruptly. In 1900 Lipót Fejér proposed taking the arithmetic mean of the partial sums, $\sigma_N = \frac{S_0 + S_1 + \dots + S_{N-1}}{N}$. In it the $k$th harmonic enters with weight $1 - \frac{k}{N}$, that is, high harmonics are muted smoothly. Switch on the Fejér means in the magnifier: there are no ripples, and the step is slightly blurred. Why there’s no overshoot at all is shown by the kernel.

Gibbs ringing is familiar to anyone who has seen over-compressed JPEG images: ripples appear around black letters on a white background. These are the discarded high harmonics, and we’ll come back to them.

Circles on circles

The formulas get shorter if we recall Euler’s formula from the chapter on the number $e$: $e^{i\varphi} = \cos\varphi + i\sin\varphi$. Then $\cos nx = \frac{e^{inx} + e^{-inx}}{2}$ and $\sin nx = \frac{e^{inx} - e^{-inx}}{2i}$, and the pair “cosine plus sine” turns into a pair of rotations with frequencies $n$ and $-n$.

The complex amplitude. Its modulus is the radius of a circle, its argument the starting angle. The link with the earlier coefficients: $c_n = \frac{a_n - ib_n}{2}$, $c_{-n} = \frac{a_n + ib_n}{2}$, $c_0 = \frac{a_0}{2}$. Rotation round the unit circle: over a period the point $e^{inx}$ makes $|n|$ turns, clockwise when $n < 0$. Example: $\cos x = \frac12 e^{ix} + \frac12 e^{-ix}$, two circles of radius $\frac12$ rotating in opposite directions. The imaginary parts of their points cancel all the time, and the sum swings along the real axis from $-1$ to $1$.

The formula also works for a closed curve in the plane: a point $z(x) = X(x) + iY(x)$ that goes round the curve in time $2\pi$ is a complex-valued function with period $2\pi$. Its Fourier series is a chain of rotating circles: the centre of each sits on the rim of the one before, and each spins at its own whole-number frequency. Astronomers knew such chains a millennium and a half before Fourier. In Ptolemy’s Almagest (around 150 AD) a planet moves on a small circle, the epicycle, whose centre moves on a large circle, the deferent. This explained retrograde motion: from time to time Mars stops in the sky and loops backwards.

Look at the “Mars from Earth” sample. Mars goes round the Sun in an almost circular orbit of radius $1.52$ astronomical units, the Earth in one of radius $1$, and in 15 years the Earth makes 15 turns and Mars almost exactly 8. The position of Mars relative to the Earth is the difference of two rotations, $z(t) = 1.52\,e^{8it} - e^{15it}$: a Fourier series with two terms. Ptolemy’s deferent is the orbit of Mars itself, and the epicycle the Earth’s orbit turned inside out. That’s why, in the simplest version of the model, the epicycle of an outer planet goes round in exactly a year. Ptolemy’s real model is more complicated: the centre of the deferent is offset, and the uniformity of the motion is measured from a special point, the equant. But in one respect the widget is honest: circles on circles can draw anything at all, even the profile of an elephant, so a good match with observations explains nothing by itself. The explanation came with Kepler’s ellipses (1609) and Newton’s law of gravitation.

The digital mixing desk

In a computer, sound is stored not as a formula but as samples: on a CD, 44,100 values a second. The integral in the formula for $c_n$ is then replaced by a sum over the samples, and we get the discrete transform.

The samples of the signal over one period (or in one window of a recording). The probe: harmonic number $k$, taken at the same $N$ points. If the $N$ samples are taken to cover the period $2\pi$, then $x_n = f\bigl(\frac{2\pi n}{N}\bigr)$, and $\frac{X_k}{N}$ is the formula for $c_k$ with the integral replaced by a sum of rectangles. Example: $N = 4$, samples $1, 0, -1, 0$ (a cosine at four points). The only non-zero samples are $x_0 = 1$ and $x_2 = -1$, so $X_1 = 1 + (-1)\cdot e^{-2\pi i \cdot 2/4} = 1 + (-1)(-1) = 2$; likewise $X_3 = 2$, while $X_0 = X_2 = 0$. The cosine is two equal “circles” numbered $1$ and $-1 \equiv 3$.

Remarkably, the discrete transform can be inverted exactly, with no questions about convergence at all. At its heart are the roots of unity from the chapter on complex numbers. Write $\omega = e^{2\pi i/N}$.

The idea: add up the terms as arrows, placing each at the end of the one before, and see that the chain closes up.

The $N$th roots of unity are the vertices of a regular $N$-gon on the unit circle; in the picture $N = 8$. The numbers $\p1{\omega^{kd}}$, $k = 0, \dots, N - 1$, are roots of unity too: each is obtained from the one before by a rotation through the angle $\frac{2\pi d}{N}$. The point $\omega^d$ can be moved around the circle. Place the arrows end to end. Each is rotated relative to the one before by the same angle $\frac{2\pi d}{N}$, and after $N$ arrows the total rotation is $2\pi d$, a whole number of turns. The chain goes round a regular polygon (perhaps several times) and returns to the start: the sum is zero. Algebraically: for $q = \omega^d$ the sum $S = 1 + q + \dots + q^{N-1}$ satisfies $qS - S = q^N - 1 = \omega^{dN} - 1 = 0$. If $q \ne 1$, then $S = 0$. If $q = 1$, which happens exactly when $d$ is a multiple of $N$, all the arrows point the same way, and $S = N$.

Computing the DFT from the definition is expensive: $N$ numbers, each with $N$ terms, $N^2$ multiplications in all. For $N = 1024$ that’s more than a million. In 1965 James Cooley and John Tukey published a method that came to be called the fast Fourier transform.

One DFT of length $N$ has been reduced to two of half the length plus $\frac N2$ multiplications by $\omega^{-k}$. Repeating until the length is one, we get $\log_2 N$ levels of $\frac N2$ multiplications each: for $N = 1024$ that’s $5120$ instead of $1{,}048{,}576$, two hundred times fewer. For $2^{20}$ samples, about 24 seconds of CD sound, the gain is more than a hundred-thousandfold. It later turned out that Gauss had devised the same trick around 1805, when he was calculating the orbits of the asteroids Pallas and Juno. He didn’t publish the note in his lifetime; it was printed in his collected works in 1866 and appreciated only in the 1980s.

Almost all digital signal processing rests on the fast Fourier transform and its relatives. The equaliser in a music player multiplies each frequency component by its own factor: turning up the bass means boosting the low frequencies. The JPEG format (a 1992 standard) divides a picture into $8 \times 8$ squares of pixels and expands each in cosines. This is the discrete cosine transform, proposed by Nasir Ahmed, T. Natarajan and K. R. Rao in 1974. The eye is poor at noticing high harmonics, so they’re stored coarsely or thrown away: hence the ringing at sharp edges. The MP3 format (a 1993 standard) does the same with sound: it breaks it down by frequency and throws away what’s masked by louder neighbouring frequencies and so can’t be heard anyway.

The mixdown

One last job remains in the studio: getting your hand in. On the first level of the trainer the coefficients are read straight off a trigonometric polynomial or guessed from the function’s parity; on the second, the polynomial first has to be obtained from products and powers of sines; on the third, the coefficients are computed by the formulas for the “sawtooth”, the “triangle”, the parabola and the square wave.

Where next

The complex form has made everything simpler. Instead of two formulas for $a_n$ and $b_n$ there’s one for $c_n$, instead of pairs of “cosine and sine” there are rotating circles, and differentiating a harmonic has become very easy: $\bigl(e^{inx}\bigr)' = in\,e^{inx}$, the derivative just multiplies it by a number. What happens if the whole of calculus is carried over into the complex numbers, looking at functions whose argument and value are both complex, differentiating them and integrating them along curves in the plane? There we’ll also find the answer to an old riddle: why the Taylor series of the function $\frac{1}{1 + x^2}$, smooth on the whole line, converges only for $|x| < 1$. That’s the chapter on complex analysis.