# Building functions out of sines and cosines Regression, as we have been doing it, is a question with three parts: a set of basis functions, a rule for combining them, and a criterion for choosing the combination. Polynomial regression uses $\{1, t, t^2, \dots\}$, combines them linearly, and picks the coefficients that minimize the sum of squared residuals. Fourier series is the same question with a different basis. The basis is the **fundamental set** $ S = \{\,1,\ \cos t,\ \sin t,\ \cos 2t,\ \sin 2t,\ \cos 3t,\ \sin 3t,\ \dots\,\}, $ and every element of it repeats every $2\pi$: $\sin 2t$ has period $\pi$, but it also repeats every $2\pi$, and so on down the list. A linear combination of $2\pi$-periodic functions is $2\pi$-periodic, so anything we build from this set, $ f(t) = a_0 + \sum_{n=1}^{\infty}\bigl[a_n \cos(nt) + b_n \sin(nt)\bigr], $ is a $2\pi$-periodic function. The term with $n = 1$ is the **fundamental**; the terms with $n = 2, 3, \dots$ are its **harmonics**, and the whole thing is a **Fourier series**. The converse is the remarkable part: essentially any $2\pi$-periodic function, including ones with corners and jumps, can be written this way. Which coefficients? That is the regression question. ## The coefficients, by orthogonality Multiply $f$ by one member of the set and integrate over a period. Every cross term vanishes, $ \int_{-\pi}^{\pi} \cos(nt)\cos(mt)\,dt = \begin{cases} 0, & n \ne m \\ \pi, & n = m \end{cases}, \qquad \int_{-\pi}^{\pi} \sin(nt)\sin(mt)\,dt = \begin{cases} 0, & n \ne m \\ \pi, & n = m \end{cases}, \qquad \int_{-\pi}^{\pi} \cos(mt)\sin(nt)\,dt = 0, $ so only the matching coefficient survives: $ a_0 = \frac{1}{2\pi}\int_{-\pi}^{\pi} f(t)\,dt, \qquad a_n = \frac{1}{\pi}\int_{-\pi}^{\pi} f(t)\cos(nt)\,dt, \qquad b_n = \frac{1}{\pi}\int_{-\pi}^{\pi} f(t)\sin(nt)\,dt . $ This is the normal equations $X^{\mathsf T}X\,\beta = X^{\mathsf T}y$ with a design matrix whose columns are *orthogonal*: $X^{\mathsf T}X$ is diagonal, so the system decouples and each coefficient is a single inner product. Two consequences follow immediately, and both are visible in the applets below. > [!note] Two things orthogonality buys you > **Adding a term never changes the ones you already have.** With polynomials, fitting a quadratic changes the linear coefficient you found when you fit a line. With the trigonometric basis, $a_1$ and $b_1$ are what they are whether you stop at $N = 1$ or go on to $N = 100$. The partial sums > $S_N(t) = a_0 + \sum_{n=1}^{N}\bigl[a_n \cos(nt) + b_n \sin(nt)\bigr]$ > are *nested*: $S_{N+1}$ is $S_N$ plus one more term. > > **The squared error is bookkeeping.** Since the basis is orthogonal, the mean square of $f$ splits into a sum over the terms, > $\frac{1}{2\pi}\int_{-\pi}^{\pi} f(t)^2\,dt = a_0^2 + \frac12\sum_{n=1}^{\infty}\bigl(a_n^2 + b_n^2\bigr),$ > which is **Parseval's identity**. The share of the left-hand side that the first $N$ terms account for is the (uncentered) $R^2$ of the regression truncated at $N$ — the draw applet below reports the centered version, with $a_0$ taken out — and the RMS residual $\|f - S_N\|$ is the square root of what is left over. ## Before Fourier: the truncated Taylor series You have built functions out of simpler pieces before. The Taylor series of $f$ about a centre $a$, $ f(t) = \sum_{n=0}^{\infty} c_n (t - a)^n, \qquad c_n = \frac{f^{(n)}(a)}{n!}, $ uses the basis $\{1, (t-a), (t-a)^2, \dots\}$, and its truncation $T_N$ is the polynomial that agrees with $f$ in value, slope, curvature and so on through the $N$th derivative *at the one point* $a$. That is where its strengths and weaknesses both come from. Near $a$ it is superb, and the error there is exactly zero. Away from $a$ the terms $c_n (t-a)^n$ grow without bound, so the polynomial eventually leaves $f$ no matter how many terms you take, and for functions with a singularity, like $\ln(1+t)$ or $1/(1-t)$, it fails everywhere outside a fixed **radius of convergence** — the distance from $a$ to the nearest bad point. A polynomial cannot be periodic, either: whatever $f$ does, $T_N$ runs off to $\pm\infty$. <div class="applet" data-applet="taylor-series"></div> Things worth trying: - Take $\sin t$ about $a = 0$ and press **build up**. The even terms vanish (sine is odd), each odd term reaches a little farther out, and it takes $N = 11$ or so before $T_N$ tracks $\sin t$ across $[-\pi, \pi]$. The Fourier series of $\sin t$ is *one term*. - Slide the centre $a$. The whole approximation moves with it: good near $a$, bad far away. A Fourier partial sum has no centre; it is equally good (or bad) everywhere on the period. - Switch to $1/(1-t)$ or $\ln(1+t)$. The orange band is $|t - a| < R$. Inside it, more terms help; outside it, $T_N$ diverges from $f$ *faster* as $N$ grows — watch the term stack, where the top strip's scale explodes. Move $a$ toward the singularity and the band shrinks to nothing. - Look at the coefficient bars. For $e^t$, $\sin t$ and $\cos t$ they die like $1/n!$, faster than anything; for $\ln(1+t)$ and $1/(1-t)$ about $a = 0$ they decay like $1/n$ or not at all. Fourier coefficients have a decay story too, and it is about smoothness on the whole period rather than at a point: that is the next section. The contrast is the point of the page. Taylor is **local**: derivatives at a point, a polynomial basis, a radius. Fourier is **global**: integrals over a period, a basis of bounded waves that already have the right periodicity, and coefficients chosen by orthogonality — least squares over the whole interval — rather than by matching derivatives at one spot. ## Four functions with closed-form coefficients The first applet builds four standard $2\pi$-periodic functions from their first $N$ terms. The integrals above can be done by hand for each; the results are the readout line under the controls. | $f(t)$ on $[-\pi, \pi]$ | $a_0$ | $a_n$ | $b_n$ | |---|---|---|---| | square, $\operatorname{sgn}(t)$ | $0$ | $0$ | $\dfrac{4}{n\pi}$ ($n$ odd), $0$ ($n$ even) | | triangle, $1 - 2\lvert t\rvert/\pi$ | $0$ | $\dfrac{8}{n^2\pi^2}$ ($n$ odd), $0$ ($n$ even) | $0$ | | sawtooth, $t/\pi$ | $0$ | $0$ | $\dfrac{2(-1)^{n+1}}{n\pi}$ | | pulse, $1$ for $\lvert t\rvert < \pi/4$, else $0$ | $\dfrac14$ | $\dfrac{2\sin(n\pi/4)}{n\pi}$ | $0$ | Three patterns in the table are worth reading before pressing anything. An **odd** function ($f(-t) = -f(t)$: square, sawtooth) has only sines; an **even** one (triangle, pulse) has only cosines, and the integrals say so before you compute them. **Half-wave symmetry** ($f(t + \pi) = -f(t)$: square, triangle) kills every even harmonic. And the **rate of decay** of the coefficients tracks the smoothness of $f$: a jump gives $1/n$ (square, sawtooth, pulse), a corner gives $1/n^2$ (triangle), and each extra derivative of smoothness buys another power of $n$. That decay rate decides how many terms you need. <div class="applet" data-applet="fourier-waves"></div> Things worth trying: - Start with the **square** wave at $N = 1$ and press **build up**. The first term is a sine of amplitude $4/\pi \approx 1.27$, already taller than the square wave it is approximating; every subsequent odd term is a smaller, faster sine that flattens the top and sharpens the sides. The even strips in the stack are flat lines: those coefficients are zero. - Watch the readout as $N$ grows. The RMS error falls, the Parseval share climbs toward $100\%$, but **max $S_N$** settles at about $1.179$ and stays there. That is the **Gibbs phenomenon**: near a jump, the partial sums overshoot by about $9\%$ of the jump no matter how many terms you take. The overshoot gets *narrower* but not *shorter*, so the pointwise error at the jump never goes away even though the RMS error goes to zero. - Switch to the **triangle** and build it up again. With $1/n^2$ decay, $N = 3$ is already visually indistinguishable from the target, and there is no overshoot because there is no jump. - The **pulse** has a nonzero mean, so $a_0$ appears at the top of the stack as a flat magenta line, and its $a_n = 2\sin(n\pi/4)/(n\pi)$ vanishes whenever $n$ is a multiple of $4$: look at the spectrum. - **Listen.** The **fundamental only** button plays $220\,\text{Hz}$; **play $S_N$** adds the first $N$ harmonics (at $220n\,\text{Hz}$) and **play $f$** adds two hundred. The square wave and the sawtooth have the same fundamental and the same $1/n$ decay but different harmonic *content* (odd-only versus all), which is why a clarinet and a violin playing the same note do not sound alike. Press **terms above $N$** to hear only what $S_N$ is leaving out; as $N$ grows, that residual gets quieter and higher. ## Draw your own The second applet takes any function you draw on one period, computes the integrals numerically (a trapezoid rule on $512$ samples, which for a periodic integrand is the natural quadrature: on a full period it is spectrally accurate, and for the harmonics it is exactly the discrete Fourier transform), and rebuilds the function from its first $N$ terms. The design matrix here has $2N + 1$ columns and $512$ rows; because the columns are orthogonal on the sample grid, the normal equations are still diagonal, and each coefficient is still one dot product. <div class="applet" data-applet="fourier-draw"></div> Things worth trying: - Draw something smooth and wavy with the mouse. Then drag $N$ down. Notice how few terms it takes: for most hand-drawn shapes, $N = 8$ captures well over $95\%$ of the variance. Most of the information in a smooth curve lives in its lowest few frequencies; that is the whole premise of compression. - Press **a step** and repeat. Now $N = 8$ still rings visibly at the two jumps, and the spectrum decays like $1/n$ instead of dropping off a cliff. Jumps are expensive in this basis. - Press **noise**. The spectrum is now flat out to $n = 64$, and no $N$ short of $64$ reconstructs the curve. That is the flip side: high-frequency content that is *not* structure costs exactly as much as content that is. When we fit a Fourier model to a time series later in the term, the coefficients we keep are the signal and the ones we drop are, we hope, the noise; the applet shows what that bet looks like when it pays off and when it does not. - Draw a curve whose right end sits well above its left end. Because the function is $2\pi$-periodic, the applet joins the ends, and the joint is a jump you did not draw. Its $1/n$ tail shows up in the spectrum. This is why, in the sea level analysis, detrending *before* the Fourier fit matters: a trend is a sawtooth in disguise. ## The complex form, for later Euler's identity $e^{i\theta} = \cos\theta + i\sin\theta$ packs each cosine–sine pair into one exponential. With $c_n = (a_n - i b_n)/2$ and $c_{-n} = \overline{c_n}$, $ f(t) = \sum_{n=-\infty}^{\infty} c_n\, e^{int}, \qquad c_n = \frac{1}{2\pi}\int_{-\pi}^{\pi} f(t)\,e^{-int}\,dt, $ and $|c_n| = \frac12\sqrt{a_n^2 + b_n^2}$ is half the amplitude of the $n$th harmonic, which is what the spectrum panels above are plotting. The complex form is the one the FFT computes, and it is the form that turns into the Fourier *transform* when the period is sent to infinity. The applets on this page live in the finite-period world; the next step is what happens to the picture when it is not. ## The applets in code Each applet reduced to the few lines that compute what it shows, in the four languages of the course. Fully working versions of these — with the checks against closed forms, the figures, and the sound written to WAV files — are kept in the `fourier-codes` repository (`python/`, `matlab/`, `r/`, `mathematica/`; one file per applet, named after it). > [!example]- `taylor-series` — the truncated Taylor polynomial $T_N$ of $\sin t$ about $a$, against $\sin t$ > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > from math import factorial > > N, a = 5, 0.0 > > t = np.linspace(-4, 4, 801) > > c = [np.sin(a + n * np.pi / 2) / factorial(n) for n in range(N + 1)] # c_n = f^(n)(a)/n! for sin > > T = sum(c[n] * (t - a)**n for n in range(N + 1)) # T_N(t) > > plt.plot(t, np.sin(t), "0.6", lw=3); plt.plot(t, T, "r"); plt.ylim(-4, 4); plt.show() > > plt.bar(range(N + 1), np.abs(c)); plt.show() # the coefficients > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > N = 5; a = 0; t = linspace(-4, 4, 801); > > c = sin(a + (0:N) * pi / 2) ./ factorial(0:N); % c_n = f^(n)(a)/n! for sin > > T = c * (t - a).^((0:N)'); % T_N(t) > > plot(t, sin(t), 'k', t, T, 'r'); ylim([-4 4]); figure; bar(0:N, abs(c)) > > ``` > > > [!info]- R > > ```r > > N <- 5; a <- 0; t <- seq(-4, 4, length.out = 801) > > c <- sin(a + (0:N) * pi / 2) / factorial(0:N) # c_n = f^(n)(a)/n! for sin > > T <- colSums(c * outer(0:N, t, function(n, t) (t - a)^n)) # T_N(t) > > plot(t, sin(t), type = "l", col = "grey", lwd = 3, ylim = c(-4, 4)); lines(t, T, col = "red") > > barplot(abs(c), names.arg = 0:N) # the coefficients > > ``` > > > [!info]- Mathematica > > ```mathematica > > NN = 5; a = 0; > > T[t_] = Normal @ Series[Sin[t], {t, a, NN}] (* T_N about a; Normal drops the O-term *) > > Plot[{Sin[t], T[t]}, {t, -4, 4}, PlotRange -> {-4, 4}, PlotStyle -> {Gray, Red}] > > ListPlot[Abs @ CoefficientList[T[t], t], Filling -> Axis] (* the coefficients (about a = 0) *) > > ``` > > [!example]- `fourier-waves` — the square wave from its first $N$ terms, its spectrum, and the sound of $S_N$ > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > N = 5 > > t = np.linspace(-np.pi, np.pi, 1001); n = np.arange(1, N + 1) > > b = np.where(n % 2, 4 / (n * np.pi), 0) # square wave: a_n = 0, b_n = 4/(n pi) for odd n > > S = b @ np.sin(np.outer(n, t)) # S_N(t) = sum b_n sin(n t) > > plt.plot(t, np.sign(t), "0.6", lw=3); plt.plot(t, S, "r"); plt.show() > > plt.stem(n, b); plt.show() # the spectrum > > ts = np.arange(44100 * 1.5) / 44100 # hear it: harmonic n at 220 n Hz with amplitude b_n > > y = b @ np.sin(2 * np.pi * 220 * np.outer(n, ts)) > > from scipy.io import wavfile; wavfile.write("S_N.wav", 44100, (0.5 * y / abs(y).max() * 32767).astype(np.int16)) > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > N = 5; t = linspace(-pi, pi, 1001); n = 1:N; > > b = 4 ./ (n * pi) .* mod(n, 2); % square wave: a_n = 0, b_n = 4/(n pi) for odd n > > S = b * sin(n' * t); % S_N(t) = sum b_n sin(n t) > > plot(t, sign(t), 'k', t, S, 'r'); figure; stem(n, b) % the partial sum, then the spectrum > > ts = (0:44100*1.5 - 1) / 44100; % hear it: harmonic n at 220 n Hz with amplitude b_n > > sound(0.2 * (b * sin(2 * pi * 220 * n' * ts)), 44100) > > ``` > > > [!info]- R > > ```r > > N <- 5; t <- seq(-pi, pi, length.out = 1001); n <- 1:N > > b <- ifelse(n %% 2 == 1, 4 / (n * pi), 0) # square wave: a_n = 0, b_n = 4/(n pi) for odd n > > S <- colSums(b * sin(outer(n, t))) # S_N(t) = sum b_n sin(n t) > > plot(t, sign(t), type = "l", col = "grey", lwd = 3); lines(t, S, col = "red") > > plot(n, b, type = "h") # the spectrum > > ts <- (0:(44100 * 1.5 - 1)) / 44100 # hear it: harmonic n at 220 n Hz with amplitude b_n > > y <- colSums(b * sin(outer(n, 2 * pi * 220 * ts))) > > tuneR::play(tuneR::Wave(round(32767 * 0.2 * y), samp.rate = 44100, bit = 16)) # or writeWave(...) > > ``` > > > [!info]- Mathematica > > ```mathematica > > NN = 5; b[n_] := If[OddQ[n], 4/(n Pi), 0]; (* square wave: a_n = 0, b_n = 4/(n Pi) for odd n *) > > S[t_] := Sum[b[n] Sin[n t], {n, NN}]; (* the partial sum *) > > Plot[{Sign[t], S[t]}, {t, -Pi, Pi}, PlotStyle -> {Gray, Red}] > > ListPlot[Table[b[n], {n, NN}], Filling -> Axis] (* the spectrum *) > > EmitSound @ Play[S[2 Pi 220 t], {t, 0, 1.5}] (* hear it: harmonic n at 220 n Hz *) > > ``` > > [!example]- `fourier-draw` — numerical coefficients of any sampled period, rebuilt with $N$ terms > > > [!info]- Python > > ```python > > import numpy as np, matplotlib.pyplot as plt > > M, N = 512, 8 > > t = -np.pi + 2 * np.pi * np.arange(M) / M # one period, endpoint not repeated > > y = np.where(t < -1, -0.6, np.where(t < 1, 0.9, -0.2)) # "the drawing": any samples you like > > n = np.arange(1, N + 1); C, Sn = np.cos(np.outer(n, t)), np.sin(np.outer(n, t)) > > a0 = y.mean(); a = 2 * C @ y / M; b = 2 * Sn @ y / M # a_n = (1/pi) int f cos(nt) etc., by the trapezoid rule > > S = a0 + a @ C + b @ Sn > > print("R^2 =", 1 - np.sum((y - S)**2) / np.sum((y - a0)**2)) > > plt.plot(t, y, "0.5", lw=3); plt.plot(t, S, "r"); plt.show() > > ``` > > > [!info]- MATLAB / Octave > > ```matlab > > M = 512; N = 8; t = -pi + 2*pi*(0:M-1)/M; % one period, endpoint not repeated > > y = -0.6*(t < -1) + 0.9*(t >= -1 & t < 1) - 0.2*(t >= 1); % "the drawing": any samples you like > > n = (1:N)'; C = cos(n * t); Sn = sin(n * t); > > a0 = mean(y); a = 2 * (y * C') / M; b = 2 * (y * Sn') / M; % a_n = (1/pi) int f cos(nt) etc., trapezoid rule > > S = a0 + a * C + b * Sn; > > R2 = 1 - sum((y - S).^2) / sum((y - a0).^2) > > plot(t, y, 'k', t, S, 'r') > > ``` > > > [!info]- R > > ```r > > M <- 512; N <- 8; t <- -pi + 2 * pi * (0:(M - 1)) / M # one period, endpoint not repeated > > y <- ifelse(t < -1, -0.6, ifelse(t < 1, 0.9, -0.2)) # "the drawing": any samples you like > > n <- 1:N; C <- cos(outer(n, t)); Sn <- sin(outer(n, t)) > > a0 <- mean(y); a <- 2 * as.vector(C %*% y) / M; b <- 2 * as.vector(Sn %*% y) / M # trapezoid rule > > S <- a0 + colSums(a * C + b * Sn) > > 1 - sum((y - S)^2) / sum((y - a0)^2) # R^2 > > plot(t, y, type = "l", col = "grey", lwd = 3); lines(t, S, col = "red") > > ``` > > > [!info]- Mathematica > > ```mathematica > > M = 512; NN = 8; t = N[-Pi + 2 Pi Range[0, M - 1]/M]; > > y = Which[# < -1, -0.6, # < 1, 0.9, True, -0.2] & /@ t; (* "the drawing": any samples you like *) > > a0 = Mean[y]; a = Table[2 Mean[y Cos[n t]], {n, NN}]; (* a_n = (1/Pi) Integral f Cos[n t], trapezoid rule *) > > b = Table[2 Mean[y Sin[n t]], {n, NN}]; > > S = a0 + Sum[a[[n]] Cos[n t] + b[[n]] Sin[n t], {n, NN}]; > > 1 - Total[(y - S)^2]/Total[(y - a0)^2] (* R^2 *) > > ListLinePlot[{Transpose[{t, y}], Transpose[{t, S}]}, PlotStyle -> {Gray, Red}] > > ``` > > [!quote] Attribution > The idea of hearing the partial sums while watching them build — and the wave-stack and draw-your-own applets this page grew from — comes from Jez Swanson's *[An Interactive Introduction to Fourier Transforms](https://www.jezzamon.com/fourier/)* (source on [GitHub](https://github.com/Jezzamonn/fourier), MIT License). The applets here are an independent rebuild in his spirit, with the mathematics in the notation of this course.