Fourier Series

The series

For a function f periodic with period 2π, the Fourier series is:

f(x) = a₀/2 + ∑_{n=1}^∞ [aₙ cos(nx) + bₙ sin(nx)]

where the coefficients are:

a₀ = (1/π) ∫_{−π}^π f(x) dx                    (twice the average)
aₙ = (1/π) ∫_{−π}^π f(x) cos(nx) dx            (cosine coefficients)
bₙ = (1/π) ∫_{−π}^π f(x) sin(nx) dx            (sine coefficients)

Each term in the sum is a sinusoid — a cosine of frequency n with amplitude aₙ, plus a sine of frequency n with amplitude bₙ. As n increases, the frequencies get higher.

For a signal of period T (other than 2π), substitute (2π/T)x for x in the formulas. The frequencies become 2πn/T — the harmonics of the fundamental frequency 1/T.

The orthogonality miracle

The reason this decomposition is unique is that sines and cosines at different frequencies are orthogonal under the inner product (1/π) ∫_{−π}^π f(x)g(x) dx:

∫_{−π}^π cos(mx) cos(nx) dx = π · δ_{mn}     (zero unless m = n)
∫_{−π}^π sin(mx) sin(nx) dx = π · δ_{mn}     (same)
∫_{−π}^π cos(mx) sin(nx) dx = 0              (always — sin and cos are orthogonal)

So the cosines and sines form an orthogonal basis. Each Fourier coefficient is the inner product of f with the corresponding basis function — that's why the formula has the integral of f against cos(nx) or sin(nx).

Example — square wave

The "square wave" — f(x) = 1 for x ∈ (0, π), f(x) = −1 for x ∈ (−π, 0). Compute Fourier coefficients:

• aₙ = 0 for all n (f is odd, cos is even, integral is zero by symmetry).
• bₙ = 4/(nπ) for odd n, 0 for even n.

So:

square wave = (4/π) [sin(x) + sin(3x)/3 + sin(5x)/5 + sin(7x)/7 + ...]

An infinite sum of odd-frequency sines reconstructs the square wave. Visualizing — the first term is a smooth sine. Adding the next term sharpens the corners. Each additional term refines further. The infinite sum is the exact square wave.

Approximating a square wave with N sine terms — visible Gibbs overshoot near the discontinuities, ~9% above the actual value. The overshoot persists no matter how many terms you add; it just gets narrower.

Complex form

Using Euler's formula e^(iθ) = cos θ + i sin θ, the Fourier series compactifies to:

f(x) = ∑_{n=−∞}^∞ cₙ · e^(inx)

with complex coefficients:

cₙ = (1/(2π)) ∫_{−π}^π f(x) · e^(−inx) dx

This is the more elegant form for theoretical work. Real coefficients aₙ and bₙ are recovered as 2 Re(cₙ) and −2 Im(cₙ) for positive n. The complex form treats positive and negative frequencies symmetrically.

Discrete Fourier Transform and FFT

For a discrete signal of N samples x₀, x₁, ..., x_{N−1}, the Discrete Fourier Transform is:

X_k = ∑_{j=0}^{N−1} x_j · e^(−2πijk/N)

Each X_k is the amplitude (and phase) of frequency component k. The inverse DFT recovers the original signal from frequency coefficients.

Naively, DFT is O(N²) — for each of N output components, sum N input samples. The Fast Fourier Transform (Cooley-Tukey 1965) computes the same DFT in O(N log N) by recursive divide-and-conquer. For N = 1 million, that's 20 million operations vs a trillion — 50,000× speedup. FFT enabled real-time digital signal processing.

JavaScript — DFT and Fourier series

// Naive DFT (O(N²) — for educational purposes)
function dft(signal) {
  const N = signal.length;
  const X = new Array(N);
  for (let k = 0; k < N; k++) {
    let real = 0, imag = 0;
    for (let n = 0; n < N; n++) {
      const angle = -2 * Math.PI * k * n / N;
      real += signal[n] * Math.cos(angle);
      imag += signal[n] * Math.sin(angle);
    }
    X[k] = { re: real, im: imag, magnitude: Math.sqrt(real*real + imag*imag) };
  }
  return X;
}

// Fourier series approximation of square wave with first N terms
function squareWaveApprox(x, terms) {
  let sum = 0;
  for (let n = 1; n < 2 * terms; n += 2) {  // only odd harmonics
    sum += Math.sin(n * x) / n;
  }
  return (4 / Math.PI) * sum;
}

// Visualize: more terms → closer approximation, but Gibbs overshoot persists
const xs = Array.from({length: 1000}, (_, i) => -Math.PI + i * 2 * Math.PI / 999);
const approxN10 = xs.map(x => squareWaveApprox(x, 10));
// Plot approxN10 — you'll see the square wave shape with ~9% overshoot at the jumps

Where Fourier analysis shows up

• Audio signal processing. MP3, AAC, Opus all use Fourier-related transforms (or modified discrete cosine transforms) to identify frequency components and discard inaudible ones. Real-time pitch detection, equalizers, noise cancellation — all FFT-based.
• Image compression. JPEG uses the Discrete Cosine Transform. JPEG-2000 uses wavelets (Fourier's spiritual successor). Compression works by aggressively quantizing high-frequency components.
• Medical imaging. MRI raw data is in Fourier space (k-space). Inverse Fourier transform reconstructs the image. CT scans use related Radon transforms.
• Radio and telecommunications. Modulation, demodulation, channel coding — all rely on frequency-domain analysis. OFDM (used in 4G, 5G, Wi-Fi) sends data on many subcarriers via FFT.
• Spectroscopy. Emission and absorption spectra are essentially Fourier transforms of atomic/molecular response functions. NMR, IR, Raman — all interpret data in the frequency domain.
• Differential equations. Heat, wave, diffusion equations on periodic domains — solutions are Fourier series. The PDE in physical space becomes an ODE in frequency space (separation of variables).
• Quantum mechanics. Position-space wavefunction and momentum-space wavefunction are Fourier transforms of each other. Plane wave decomposition uses Fourier series. Time-evolution involves Fourier on time.

Common mistakes

• Forgetting the period adjustment. The standard Fourier series formulas assume period 2π. For period T, frequencies become 2πn/T and the integration is over (−T/2, T/2). Forgetting the substitution gives wrong frequencies.
• Not noticing Gibbs phenomenon. Truncated Fourier series of discontinuous functions overshoot near jumps by ~9%. This is real and can't be eliminated by adding more terms — it just gets narrower. Visible in image compression artifacts.
• Confusing series with transform. Series for periodic functions (discrete sum); transform for non-periodic functions (continuous integral). They're related but different mathematical objects.
• Using DFT when FFT is available. Naive DFT is O(N²); FFT is O(N log N). Using DFT for large N is a serious performance bug. NumPy and most libraries default to FFT.
• Wrong normalization conventions. Different texts and libraries use different normalizations (some divide by N, some by √N, some not at all). Inverse Fourier should give back the original signal; if it doesn't, check the convention.
• Not zero-padding for finer frequency resolution. An N-point signal has N frequency bins. To get more bins (interpolated frequencies), zero-pad before FFT — but the underlying frequency resolution is still 1/T (where T is signal duration), not 1/(NT_padded).