# 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. ## 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]- `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.