LLGA · Course notes

Signal Processing

APM_51055_EP · Year 1 · Period 1

Contents

Introduction and overview

Module Signal Processing — SynapseS 2026-2027 sheet: “Programme d’approfondissement — Traitement du Signal : de Fourier à l’Apprentissage Machine” (occurrence: “Signal processing: from Fourier to Machine Learning”, X3A VE1, applied mathematics department)
Code APM_51055_EP
Team and workload Responsible (sheet and 2026-2027 occurrence): Rémi Flamary · teaching team: Mansour Benbakoura, Jules Désir · 36 h in person (no lecture / tutorial / lab split displayed; max enrollment 60) · 4.5 ECTS (value kept by these notes for the LLGA, whose section of the sheet displays none; the sheet shows 5 ECTS for the École Polytechnique degree, the DAIIG and MJH masters, exchange programs and non-degree students, 4.5 for ViCAI and TRAI, 4 for the MScT Artificial Intelligence and Advanced Visual Computing) · assessment: lab reports (40%) and a final written theoretical and practical exam (60%), as stated in the description · teaching languages: French and English (“in English or in French depending on the audience, with course materials in English”) · signals, Fourier, sampling, filtering, time-frequency analysis, spectral estimation — up to the spectrograms consumed by ML and the spectral distances of optimal transport (the same instructor’s M1P2 module).
Status in the program P1 electiveAI-systems route M1 · S1 / period 1 — X-3A-P1 (fall semester, from 17/09/26 to 18/12/26) · one of the “Choose 3” electives of Year 1 · Period 1, preselected by the AI-systems route of the home page

Audience: an engineer who consumes audio / series / images in ML and wants the chain before the network. Empirical first: an unreadable spectrogram from a bad window; aliasing that invents a frequency — then we name LTI, Fourier, Shannon, STFT.

Bridge to speech, spectral OT, sensors: the raw material of modern AI.

What you will be able to do

Official objectives and syllabus (2026-2027 catalog)

Official objectives (paraphrased from the description published on SynapseS): historically driven by telecommunications, signal processing is today present in every field that analyzes, extracts and transforms digital information. The course introduces the field — Fourier analysis and analog filtering (with modulation and Fourier optics as examples), sampling and digital filtering which have become the standard of applications, the fast Fourier transform algorithm, image filtering, then the random aspects of signals and the optimal linear filtering of the signal/noise pair modeled as a stochastic process (autoregressive models, Wiener filtering).

It concludes with the common representations of signals — the discrete cosine transform and wavelets (JPEG encoding, image reconstruction), the short-time Fourier transform for non-stationary signals — and with recent approaches coming from machine learning: learned dictionaries and signal reconstruction by deep learning. All of it is complemented by Python/NumPy labs (generating an audio signal, filtering). Official assessment: lab reports (40%) and a final theoretical and practical written exam (60%).

Block of the official syllabus Course section that covers it
Fourier analysis and analog filtering (Fourier transform, convolution and filtering, applications: modulation, Fourier optics) ch. 3 (transform: definition 3.1, propositions 3.2–3.3), ch. 1 (convolution and LTI systems), ch. 2 (series); applications: example 3.6 (amplitude modulation and demodulation) and example 3.7 (Fourier optics: diffraction, telescope resolution, interferometry)
Digital signal processing (sampling and discrete signals, Z-transform and transfer function, fast Fourier transform, applications to signals and images) ch. 4 (Shannon, aliasing, acquisition chain), ch. 5: DFT and FFT (theorem 5.3), the Z-transform and transfer function with a first-order recursive filter worked out, FIR/IIR filters and stability by the poles (theorem 5.5), filter design (example 5.6); applications: example 5.7 (denoising a signal) and example 5.8 (image filtering: 2D DFT, blur, edge detection)
Random signals (correlations and representation, filtering and prediction of stationary signals, autoregressive process and Wiener filtering) ch. 7: WSS processes and autocorrelation (definition 7.1, proposition 7.2), Wiener-Khinchin and filtered noise (theorem 7.3), AR(1), periodogram and Bartlett/Welch (theorem 7.5), linear prediction, AR(p) and the Yule–Walker equations with a worked AR(2) fit, and the Wiener filter (note “Complement to the official syllabus”), and the Kalman filter, its recursive state-space form
Signal representation and learning dictionary (non-stationary signals and short-time FT, common representations: Fourier, wavelets; source separation and learning dictionary, machine learning for signal processing) ch. 6: STFT (definition 6.1), wavelets and the Haar transform (example 6.4), JPEG 2000, the DCT and JPEG (definition 6.5, example 6.6); ch. 8: Mel spectrograms and MFCC, networks on spectrograms, spectral distances, Fourier in architectures, source separation and dictionary learning (note “Complement to the official syllabus” and example 8.2: sparse coding, NMF, soft masks, unrolled networks)
Python/NumPy labs (generating an audio signal, filtering) examples and NumPy code excerpts of ch. 4 to 8, and the exercises

Our materials cover every topic of the four blocks of the official syllabus. Several were added after checking them against the sheet: amplitude modulation and Fourier optics (examples 3.6–3.7), the Z-transform, image filtering (example 5.8), the DCT and JPEG (definition 6.5, example 6.6), linear prediction and Yule–Walker (ch. 7), Wiener filtering (ch. 7) and dictionary learning and source separation (ch. 8). In return, these materials go beyond the sheet: ch. 0 (L²/Fubini review) and ch. 9 (the instructor’s research: optimal transport) are deepenings specific to the LLGA track, and ch. 8 explicitly connects spectral representation and ML architectures.

Official field Published value (SynapseS 2026-2027)
Fields Applied Mathematics
Workload 36 h in person (no lecture / tutorial / lab split displayed) · max enrollment 60
Period / language X-3A-P1, fall semester, period 1 (17/09/26 – 18/12/26) · French, English (course materials in English)
Responsible Rémi Flamary (also responsible of the occurrence “X3A VE1”) · teaching team: Mansour Benbakoura, Jules Désir
Prerequisites The prerequisites tab is empty; the description states that the course requires basic knowledge of analysis (Fourier transform), probability (random variables, random processes) and linear algebra, and strongly recommends Python/NumPy for the labs
Assessment Lab reports (40%) and a final theoretical and practical written exam (60%), stated in the description and in the sections of the X, AI&AVC, DAIIG, exchange and non-degree students · grade format: numeric out of 20 (and US letter grade) · the section for the LLGA track is empty: no assessment form, resit rule, pass threshold or ECTS is displayed for LLGA · for the other degrees: resit allowed (max of the two grades, or resit grade kept for the DAIIG and MJH masters), unit acquired if final grade ≥ 10, grade counted in the GPA; ECTS 5 (X, DAIIG, MJH, exchange programs, non-degree students), 4.5 (ViCAI, TRAI), 4 (AI&AVC) · these notes keep 4.5 ECTS for LLGA

Where this course fits in the program

One of the “Choose 3” electives of Year 1 · Period 1 (36 h, 4.5 ECTS as kept by these notes), preselected by the AI-systems route of the home page.

The official prerequisites of the description — basic knowledge of analysis (Fourier transform), probability (random variables and processes) and linear algebra, plus Python/NumPy for the labs — correspond exactly to the foundations of the two refresher modules: Refresher in Statistics for the probability theory useful for random signals, and Refresher in Computer Science for NumPy and the implementation of the labs; the review chapter (ch. 0) moreover takes up Fourier theory at the required level.

The course unit is shared with the engineering cycle (X), the DAIIG and MJH master’s programs and the other MScT tracks (AI&AVC, ViCAI, TRAI).

Students of the AI-systems route take this course before any deep-learning course, since they take Deep Learning in its P2 variant. The few uses of convolutional networks, transformers, autoencoders, graph convolutional networks and optimal transport in these notes are therefore labelled as previews, and each is defined in two lines in the box that opens ch. 8.

Downstream, three modules directly consume this course. Speech Technology for Conversational AI (M2) starts from the speech signal representations introduced here; Multimodal Generative AI leverages spectrograms as the interface between audio and models; Computational Optimal Transport for ML and Fairness, a module by the same instructor, uses the spectral distances (sinkhorn on spectrograms) whose foundations are laid in the chapter “The signal for modern ML”. The course also sheds light on the filtering and convolution that Image Analysis and Computer Vision mobilizes from its first chapters.

In Period 2, Graph ML, ch. 4 carries the Fourier transform and the filters of ch. 3–5 over to signals on the nodes of a graph, and the Kalman filter of ch. 7 is the motion model of the trackers of Computer Vision Applications, ch. 5.

Course outline

Session plan

36 h in class = 12 sessions of 3 h. The volume is the one of the SynapseS 2026–2027 sheet (36 h in person, no lecture/tutorial/lab split displayed).

The sheet gives no session-by-session timetable, but its detailed programme lists four blocks in order, and the split below follows it — analog Fourier analysis (sessions 1–4), digital processing (sessions 5–7), then representations and random signals, and machine learning last (sessions 11–12) — with one exception: time-frequency representations (block 4 on the sheet, ch. 6) are taught before random signals (block 3, ch. 7), because the STFT reuses the windows of ch. 5 directly and the two chapters are otherwise independent.

Level 1 = foundations (reviews, definitions, first worked examples), 2 = core methods (transforms, sampling, filtering, spectral estimation), 3 = advanced material (the bridge to current machine learning); the level never goes down from one session to the next.

Session Hours Level Content Practice
1 3 h 1 ch. 0 Review: L^2, orthogonality of the exponentials, Fejér lab: partial sums and Fejér means of a square wave
2 3 h 1 ch. 1 Signals and LTI systems: convolution, BIBO stability, frequency response lab: generating and filtering an audio signal
3 3 h 1 ch. 2 Fourier series: projection, Parseval, decay of the coefficients, Gibbs Ex. 1, 2, 8
4 3 h 2 ch. 3 Fourier transform: properties, Gaussian, sinc, Heisenberg-Gabor, Dirac comb; applications: amplitude modulation, Fourier optics Ex. 9, 10
5 3 h 2 ch. 4 Sampling: Shannon theorem, aliasing, acquisition chain Ex. 3
6 3 h 2 ch. 5 DFT, circular convolution, FFT, windows Ex. 4, 5
7 3 h 2 ch. 5 (continued) Z-transform and transfer function, FIR and IIR filters, stability by the poles, filter design, image filtering (2D DFT, blur, edges) Ex. 6
8 3 h 2 ch. 6 STFT, Gabor frames, inversion by overlap-add (COLA), multi-resolution STFT losses, wavelets and the Haar transform, DCT and JPEG lab: spectrogram of a chirp, effect of the window length
9 3 h 2 ch. 7 Random signals: stationarity, autocorrelation, Wiener-Khinchin, AR(1), the periodogram and its bias lab: simulating an AR(1) noise and its spectrum
10 3 h 2 ch. 7 (continued) Bartlett and Welch, linear prediction and Yule–Walker, Wiener filtering, the Kalman filter as recursive Gaussian conditioning Ex. 11
11 3 h 3 ch. 8 Preview of the learned models, Mel spectrogram, MFCC, networks on spectrograms Ex. 7
12 3 h 3 ch. 8 (continued) Spectral distances (Fréchet, optimal transport), Fourier in architectures, dictionary learning and NMF source separation Ex. 12 (project)

Self-study, outside the 36 h (about 2 h in total): ch. 9 The teacher’s research — POT, transport between spectrograms, sliced transport (≈ 2 h).

How to study this module

Independent note The conversational and generative models at the heart of the LLGA master’s program consume signals: speech, music, sensors. This module provides the spectral culture shared by the engineers who will deploy these systems and the researchers who will invent their representations — from the spectrogram to transport distances.

0. Review: Fourier series in the L² sense, Fubini

Every chapter of this course decomposes a signal on the complex exponentials e^{ik\omega_0 t}, and almost every proof exchanges a sum with an integral. This review fixes the two tools that make these manipulations legitimate: the geometry of the Hilbert space L^2, in which the exponentials form an orthonormal basis, and Fubini’s theorem.

L3 background used

  • The Hilbert space L^2: hermitian inner product \langle f, g\rangle = \int f(t)\,\overline{g(t)}\,dt, Cauchy-Schwarz inequality, orthogonal projection onto a finite-dimensional subspace — Fourier = a change of basis in an infinite orthonormal basis;
  • Fubini’s theorem: if \iint |f(x,y)|\,dx\,dy < \infty (or for positive functions), the order of integration is free — it justifies all the convolution computations of chapter 3;
  • Integration by parts and changes of variable;
  • Complex geometric sequences, which give the orthogonality of the exponentials.

Proposition 0.1 — The exponentials form an orthogonal system On [0, T_0], set e_k(t) = e^{ik\omega_0 t} with \omega_0 = 2\pi/T_0. Then \langle e_k, e_\ell\rangle = \int_0^{T_0} e^{i(k-\ell)\omega_0 t}\,dt = \begin{cases} T_0 & \text{if } k = \ell,\\[2pt] \dfrac{e^{i(k-\ell)\omega_0 T_0} - 1}{i(k-\ell)\omega_0} = 0 & \text{if } k \neq \ell,\end{cases} since e^{i(k-\ell)\omega_0 T_0} = e^{2i\pi(k-\ell)} = 1. The e_k / \sqrt{T_0} therefore form an orthonormal system of L^2([0,T_0]).

Definition 0.2 — The Dirichlet and Fejér kernels The partial sum of order N of a Fourier series can be written as a convolution with the Dirichlet kernel D_N(\theta) = \sum_{k=-N}^{N} e^{ik\theta} = \frac{\sin\big((N+\tfrac12)\theta\big)}{\sin(\theta/2)}, while the Cesàro mean of the partial sums (arithmetic mean S_0, \dots, S_{N-1}) is written with the Fejér kernel F_N = \frac{1}{N}\sum_{n<N} D_n, which has the decisive property of being positive (F_N \geq 0) and of mass 1.

Theorem 0.3 — Fejér convergence, density, completeness

  • Fejér: for f continuous and T_0-periodic, \sigma_N f = f * F_N \to f uniformly. The proof uses F_N \geq 0, \int F_N = 1 and the concentration of the mass of F_N near 0.
  • Consequence (density): the trigonometric polynomials \sum_{|k| \leq N} c_k e^{ik\omega_0 t} are dense in L^2([0,T_0]).
  • Consequence (Hilbert basis): the orthonormal system \{e_k/\sqrt{T_0}\} is complete: it is a Hilbert basis of L^2. The whole theory of Fourier series “in the L^2 sense” follows from it.

Dirichlet versus Fejér — the message to remember The partial sums themselves (Dirichlet kernel, which changes sign) do not converge uniformly for every continuous function: they “oscillate” near discontinuities (Gibbs phenomenon, quantified in ch. 2). The Cesàro average smooths these oscillations. Moral: in L^2, everything converges nicely; uniformly, one must regularize (averaging, windows) — exactly the role of the Hann windows of ch. 5.

1. Signals and linear systems

Concrete scene: a fine-tuned speech model whose accuracy drops because the Mel spectrogram was badly parameterized upstream. Before any spectrogram, one needs a precise language for signals and for the systems that transform them. This chapter sets it up with linear time-invariant (LTI) systems, which are entirely described by one convolution.

What breaks Three classic failures: sampling below the Nyquist rate, choosing an analysis window too short for the frequency resolution one needs, and treating filtering as a black box inside the ML pipeline.

Definition 1.1 — Signal, LTI system A signal is a function x : \mathbb{R} \to \mathbb{R} (or \mathbb{C}), or a sequence x : \mathbb{Z} \to \mathbb{R} (discrete signal).

A system transforms one signal into another; it is LTI (linear time-invariant) if x \mapsto y = T(x) satisfies: linearity, and x(\cdot - \tau) \mapsto y(\cdot - \tau) for every delay \tau.

We write \delta_\tau for the Dirac mass at \tau (discrete: \delta_{n} = the sequence equal to 1 at n, 0 elsewhere).

Proposition 1.2 — Convolution representation (proved) Every LTI system (discrete, and continuous subject to a continuity condition) is entirely described by its impulse response h = T(\delta_0): y(t) = (x * h)(t) = \int_{\mathbb{R}} x(\tau)\, h(t - \tau)\, d\tau \qquad \Big(\text{discrete: } y[n] = \sum_k x[k]\,h[n-k]\Big).

Proof. Discrete version (the continuous argument replaces the sums by integrals against Dirac masses, x = \int x(\tau)\,\delta_\tau\,d\tau, which requires a minimum of distribution theory — admitted here).

In discrete time, x = \sum_k x[k]\,\delta_k (a sum of “column” sequences). By linearity, T(x) = \sum_k x[k]\,T(\delta_k). By invariance, T(\delta_k) = \tau_k h, where (\tau_k h)[n] = h[n-k]. Hence y[n] = \sum_k x[k]\,h[n-k]. Commutativity follows from the change of variable \tau' = t - \tau (Fubini: the integral is absolutely convergent under the stability assumptions below).∎

Theorem 1.3 — BIBO stability and causality (proved) The LTI system with impulse response h is:

  • BIBO stable (bounded input ⇒ bounded output) if and only if \sum_n |h[n]| < \infty (continuous: \int |h| < \infty);
  • causal (the output at time t depends only on the past) if and only if h[n] = 0 for n < 0.

Proof.

  • Stability (⇐). If \|x\|_\infty \leq M then |y[n]| \leq \sum_k |x[k]|\,|h[n-k]| \leq M \sum_k |h[k]|: the convolution series converges absolutely.
  • Stability (⇒). Suppose \sum |h[n]| = \infty and take the bounded input x[k] = \mathrm{sign}(h[-k]) (\|x\|_\infty \leq 1). Then y[0] = \sum_k x[k]\,h[-k] = \sum_k |h[k]| = \infty: the output explodes.
  • Causality (⇒). If h[n_0] \neq 0 for some n_0 < 0, the response to \delta_0 (a causal input, zero for t < 0) would be nonzero at n_0 < 0, before the impulse: contradiction.
  • Causality (⇐). If h[n] = 0 for n < 0, the sum y[n] = \sum_{k \leq n} x[k]\,h[n-k] involves only x[0], \dots, x[n].∎

Proposition 1.4 — Frequency response: the exponentials are eigenvectors (proved) For a stable LTI system, the input x(t) = e^{i2\pi\nu t} (a pure frequency) produces y(t) = H(\nu)\,e^{i2\pi\nu t} with H(\nu) = \hat h(\nu) = \int h(\tau)\, e^{-i2\pi\nu\tau}\,d\tau. H is the frequency response: the system merely multiplies each frequency by a complex number (gain |H(\nu)|, phase \arg H(\nu)).

Proof. y(t) = \int h(\tau)\,x(t-\tau)\,d\tau = \int h(\tau)\,e^{i2\pi\nu(t-\tau)}\,d\tau = e^{i2\pi\nu t} \int h(\tau)\,e^{-i2\pi\nu\tau}\,d\tau. The exponential comes out unchanged up to a complex factor: it is an eigenvector, with eigenvalue H(\nu). This computational triviality is the heart of the whole course: filtering = multiplying the spectrum.∎

Example 1.5 — Elementary convolutions

  • Moving average (h = \frac1M \mathbf{1}_{[0,M]}): smooths, frequency response in \mathrm{sinc};
  • discrete differentiation (h = \delta_{1} - \delta_0): H(\nu) = e^{-i2\pi\nu} - 1, amplifies the high frequencies;
  • integrator (h = \mathbf{1}_{t \geq 0}): unstable in the BIBO sense (\int|h| = \infty).

A filter = a choice of h (or of H); the whole of chapter 5 consists in designing h from a desired H.

A computation by hand fixes the mechanics. Take the two-point average h = (\tfrac12, \tfrac12) (that is, h[0] = h[1] = \tfrac12) and the input x = (1, 2, 3). The formula y[n] = \sum_k x[k]\,h[n-k] gives y[0] = \tfrac12, \quad y[1] = \tfrac12(2 + 1) = 1.5, \quad y[2] = \tfrac12(3 + 2) = 2.5 \quad\text{and}\quad y[3] = \tfrac32: the output is one sample longer than the input, and each value is the mean of the current and previous inputs.

Its frequency response is H(\nu) = \tfrac12(1 + e^{-i2\pi\nu}), so |H(0)| = 1 (a constant passes unchanged) and |H(\tfrac12)| = 0 (the fastest alternation 1, -1, 1, \dots is erased exactly): a first low-pass filter.

2. The Fourier series

A periodic signal — a sustained note, the rotation of a motor, a periodic texture — is entirely described by the amplitudes and phases of its harmonics. This chapter computes these amplitudes (the Fourier coefficients), proves that the energy is the same in time and in frequency (Parseval), and links the smoothness of a signal to the speed at which its coefficients decay.

Theorem 2.1 — Fourier series decomposition (coefficients by projection) Let f be T_0-periodic, square integrable over one period. Then f = \sum_{k \in \mathbb{Z}} c_k\, e^{ik\omega_0 t} \quad \text{in } L^2, \qquad c_k = \frac{1}{T_0}\int_0^{T_0} f(t)\, e^{-ik\omega_0 t}\, dt, the convergence taking place in L^2 norm (and pointwise toward \frac{f(x^+)+f(x^-)}{2} at discontinuity points if f is piecewise C^1: Dirichlet’s theorem).

The c_k form the spectrum: |c_k| magnitude, \arg c_k phase.

Proof. The e_k = e^{ik\omega_0t} are orthogonal (prop. 0.1) and form a Hilbert basis (completeness, th. 0.3). The decomposition in an infinite orthonormal basis follows the logic of finite dimension: the partial sum S_N f = \sum_{|k|\leq N} c_k e_k is the orthogonal projection of f onto \mathrm{Vect}\{e_k\}_{|k|\leq N}, since \langle f - S_Nf, e_j\rangle = 0 for |j| \leq N.

Since the basis is complete, \|f - S_N f\|_{L^2} \to 0. Explicitly: for every trigonometric polynomial p = \sum_{|k|\leq N} d_k e_k, Pythagoras gives \|f - p\|^2 = \|f - S_Nf\|^2 + \|S_Nf - p\|^2 \geq \|f - S_Nf\|^2: the projection is the best approximation — and the d_k that cancel the cross terms are exactly d_k = \langle f, e_k\rangle / \|e_k\|^2 = c_k.∎

Theorem 2.2 — Parseval’s identity (proved by density) \frac{1}{T_0}\int_0^{T_0} |f(t)|^2\,dt = \sum_{k \in \mathbb{Z}} |c_k|^2.

Proof. Pythagoras on f = S_N f + (f - S_N f) with \langle S_Nf, f - S_Nf\rangle = 0: \|f\|^2 = \|S_Nf\|^2 + \|f - S_Nf\|^2. Now \|S_Nf\|^2 = \sum_{|k|\leq N}|c_k|^2\|e_k\|^2 = T_0 \sum_{|k|\leq N}|c_k|^2, and \|f - S_Nf\|^2 \to 0 by theorem 2.1 (whose proof rests on Fejér density, th. 0.3). Dividing by T_0 and passing to the limit, the “time” energy equals the “frequency” energy.∎

Lemma 2.3 — Regularity and decay of the coefficients (by integrations by parts) If f is T_0-periodic, of class C^{p-1}, and f^{(p-1)} is piecewise C^1 (so that f^{(p)} is defined except at finitely many points, and integrable), then c_k = \frac{1}{(ik\omega_0)^p}\, c_k\big(f^{(p)}\big), \qquad k \neq 0. In particular:

  • a piecewise C^1 signal with a jump [f] ⇒ decay in 1/k: c_k \sim \frac{[f]}{T_0\,ik\omega_0} = \frac{[f]}{2i\pi k};
  • f continuous and piecewise C^2 ⇒ c_k = O(1/k^2);
  • f of class C^p ⇒ c_k = o(1/k^p).

Proof. Integrate by parts the coefficient of f' over one period. The boundary term \frac{1}{T_0}\big[ f(t)\, e^{-ik\omega_0t}\big]_0^{T_0} vanishes because f is continuous and periodic, so c_k(f') = \frac{1}{T_0}\int_0^{T_0} f'(t)\,e^{-ik\omega_0t}dt = \frac{ik\omega_0}{T_0}\int_0^{T_0} f(t)\,e^{-ik\omega_0t}dt = ik\omega_0\, c_k(f). Iterating p times (every derivative up to f^{(p-1)} is continuous) gives the formula. Since f^{(p)} is integrable, its coefficients tend to 0 (Riemann–Lebesgue lemma), hence c_k = o(k^{-p}).

If instead f has a jump [f] = f(0^+) - f(0^-) at t = 0, the boundary term no longer vanishes and contributes exactly \frac{[f]}{T_0\,ik\omega_0} to c_k. The remaining term is \frac{1}{ik\omega_0}c_k(f') = o(1/k), so the jump dictates the 1/k decay.∎

Example 2.4 — The square wave: explicit coefficients and the Gibbs phenomenon For the odd square wave f = \mathrm{sign}(\sin) (value ±1, jump of height 2): the direct computation (exercise 1, with solution) gives c_k = 0 \ (k \text{ even}), \qquad c_k = \frac{-2i}{\pi k} \ (k \text{ odd}) \quad\Longleftrightarrow\quad f(t) = \frac{4}{\pi}\sum_{m \geq 0} \frac{\sin\big((2m+1)\omega_0 t\big)}{2m+1}. Decay in 1/k, consistent with lemma 2.3 (discontinuity).

The Gibbs phenomenon: near the jump, the partial sum S_N f overshoots the target by about \Big(\frac{1}{\pi}\int_0^{\pi} \frac{\sin u}{u}\,du\Big) - \frac12 \approx 8.95\,\% of the height of the jump, whatever N: the oscillation concentrates (width in O(1/N)) but its amplitude does not tend to 0.

Origin: the Dirichlet kernel changes sign (ch. 0); by averaging (Fejér) or by windowing (ch. 5), the overshoot disappears — at the price of a widened transition.

Why these three results structure the whole course Parseval: the energy can be read in the spectrum (noise, filtering, SNR).

Lemma 2.3: a smooth signal has a fast-decaying spectrum — the foundation of compression (one throws away the high frequencies) and of sampling (a fast-decaying spectrum is quickly negligible: ch. 4).

Gibbs: any abrupt cutoff in frequency costs oscillations in time — hence the soft windows of ch. 5.

3. The Fourier transform

Most signals are not periodic: a spoken word, a transient, an image patch. The Fourier transform extends the Fourier coefficients to such signals, replacing the discrete set of harmonics by a continuum of frequencies. This chapter lists its properties — the most important being that convolution becomes a product — and proves the uncertainty principle, which limits how concentrated a signal can be in time and in frequency at once.

Definition 3.1 — The Fourier transform For x \in L^1(\mathbb{R}) (then L^2 by density): \hat x(\nu) = \int_{\mathbb{R}} x(t)\, e^{-i 2\pi \nu t}\, dt, \qquad x(t) = \int_{\mathbb{R}} \hat x(\nu)\, e^{i 2\pi \nu t}\, d\nu. \nu in Hz (or \omega in rad/s with a factor 2\pi); the “continuous” transform generalizes the Fourier coefficients to non-periodic signals: a signal of infinite duration has a continuous spectrum.

Proposition 3.2 — Fundamental properties (proved)

Time Frequency Proof (sketch)
delay x(t - \tau) e^{-i2\pi\nu\tau}\hat x(\nu) change of variable u = t-\tau
dilation x(at) \frac1{|a|}\hat x(\nu/a) change of variable + axis reversal if a < 0
modulation x(t)e^{i2\pi\nu_0t} \hat x(\nu - \nu_0) the exponentials compose
derivative x' i2\pi\nu\,\hat x integration by parts, vanishing boundary terms
convolution x * h product \hat x\,\hat h Fubini (see proof)
product x\,y convolution \hat x * \hat y duality of the previous one
Gaussian e^{-\pi t^2} e^{-\pi\nu^2} see proposition 3.3

Proof. Convolution (via Fubini). \widehat{x*h}(\nu) = \iint x(\tau)h(t-\tau)e^{-i2\pi\nu t}dt\,d\tau. The modulus of the integrand is |x|\,|h| \in L^1(\mathbb{R}^2) (a product of integrable functions), so Fubini allows the exchange; set u = t - \tau: \int x(\tau)e^{-i2\pi\nu\tau}d\tau \cdot \int h(u)e^{-i2\pi\nu u}du = \hat x\,\hat h.

Differentiation. By parts, \widehat{x'}(\nu) = [x e^{-i2\pi\nu t}] + i2\pi\nu\int x e^{-i2\pi\nu t}; the boundary terms vanish for x the derivative of an L^1 function.

Duality. \mathcal{F}(\mathcal{F}x)(t) = \int \hat x(\nu) e^{-i2\pi\nu t}d\nu = x(-t) (inversion formula evaluated at -t): Fourier² reverses time — hence product ↔︎ convolution by the previous property applied twice.

Dilation. With u = at: \int x(at)\,e^{-i2\pi\nu t}dt = \frac{1}{|a|}\int x(u)\,e^{-i2\pi(\nu/a)u}du = \frac{1}{|a|}\hat x(\nu/a). ∎

Proposition 3.3 — The Gaussian is its own transform (computed) \mathcal{F}\big(e^{-\pi t^2}\big)(\nu) = e^{-\pi \nu^2}.

Proof. Set g(t) = e^{-\pi t^2} and consider \hat g. We exploit the differential equation g' = -2\pi t\,g. Integrating by parts, \int t\,g(t)\,e^{-i2\pi\nu t}dt = -\frac{1}{2\pi}\int g'(t)\,e^{-i2\pi\nu t}dt = -\frac{1}{2\pi}\,i2\pi\nu\,\hat g = -i\nu\,\hat g. Hence \hat g'(\nu) = -i2\pi \int t g e^{-i2\pi\nu t} = -i2\pi(-i\nu)\hat g = -2\pi\nu\,\hat g — the same differential equation. Solution: \hat g(\nu) = \hat g(0)\,e^{-\pi\nu^2}, and \hat g(0) = \int g = 1 (Gaussian integral).∎

The key duo: convolution ↔︎ multiplication Filtering = convolving in time = multiplying in frequency:

Key result — filtering is a product in frequency y = x*h \iff \hat y = \hat x\,H.

Here H(\nu) = \hat h(\nu) is the frequency response (prop. 1.4). Low-pass: H \approx 1 for |\nu| < \nu_c, \approx 0 beyond.

This is also why the FFT speeds up large convolutions in convolutional networks: complexity \mathcal{O}(n \log n) instead of \mathcal{O}(n^2). Such networks (CNNs) are only a preview here, defined in two lines in the box that opens ch. 8 and taught in Deep Learning, ch. 4 or its P2 variant. With proposition 3.3, a “narrow” Gaussian in time is “wide” in frequency — a foretaste of the next theorem.

Worked example — the rectangular pulse and its sinc Let x = \mathbf{1}_{[-T/2,\, T/2]} be a pulse of duration T. Direct integration gives \hat x(\nu) = \int_{-T/2}^{T/2} e^{-i2\pi\nu t}\,dt = \frac{e^{-i\pi\nu T} - e^{i\pi\nu T}}{-i2\pi\nu} = \frac{\sin(\pi\nu T)}{\pi\nu} = T\,\mathrm{sinc}(\nu T). Three readings.

  1. The spectrum first vanishes at \nu = \pm 1/T: a pulse of T = 1 ms occupies a main lobe of \pm 1 kHz, and shortening it to 0.1 ms widens the lobe to \pm 10 kHz — duration and bandwidth vary inversely, as theorem 3.4 will make precise.
  2. The value at zero, \hat x(0) = T, is the area of the pulse, as it must be since \hat x(0) = \int x.
  3. |\hat x(\nu)| decays only like 1/(\pi|\nu|), because x jumps: this is the continuous analogue of lemma 2.3, and the reason why the rectangular window of chapter 5 leaks so much.

By duality, the same computation gives the impulse response of the ideal low-pass filter, used in the proof of the Shannon theorem (ch. 4).

Theorem 3.4 — Heisenberg-Gabor uncertainty (proved) For x of energy 1, with standard deviations \sigma_t (time, around t_0) and \sigma_\nu (frequency): \sigma_t \, \sigma_\nu \ \geq\ \frac{1}{4\pi}, with equality if and only if x is a Gaussian (modulated and translated).

Proof. Time translations and frequency modulations do not change the product \sigma_t\sigma_\nu (delay = pure phase; modulation = a translation of the spectrum): one assumes t_0 = \nu_0 = 0, \|x\|_2 = 1.

Let us derive the key identity by parts (vanishing boundary terms): \int t\,\frac{d}{dt}|x|^2 dt = -\int |x|^2 = -1, \qquad\text{and}\qquad \frac{d}{dt}|x|^2 = 2\mathrm{Re}(x\,\overline{x'}), hence 2\,\mathrm{Re}\langle t\,x, x'\rangle = -1, i.e. |\langle tx, x'\rangle| \geq 1/2.

Cauchy-Schwarz: \tfrac12 \leq \|t x\|_2\,\|x'\|_2 = \sigma_t\,\|x'\|_2. Finally \|x'\|_2^2 = \|\widehat{x'}\|_2^2 = \|i2\pi\nu\,\hat x\|_2^2 = (2\pi\sigma_\nu)^2 (Parseval + differentiation). Hence \sigma_t\sigma_\nu \geq 1/(4\pi).

Equality: Cauchy-Schwarz must be saturated, x' = \lambda t x: the solutions are the Gaussians e^{-\pi t^2/\sigma^2} (up to a linear phase and a translation).∎

Definition 3.5 — Useful distributions: Dirac and comb Without formalism: the Dirac mass \delta is “the function” of mass 1 concentrated at 0, defined by \int f\,\delta = f(0) (a distribution, i.e. a continuous linear form on test functions).

We admit: \hat\delta = 1, \hat 1 = \delta (duality), and the Dirac comb \,\text{Ш}_{T_e}(t) = \sum_n \delta(t - nT_e), whose transform is a comb: \widehat{\text{Ш}_{T_e}}(\nu) = \frac{1}{T_e}\,\text{Ш}_{1/T_e}(\nu). These objects make rigorous: the spectrum of a sinusoid (two \delta lines), sampling (a product with a comb, ch. 4), white noise (ch. 7).

Applications of analog processing: modulation and Fourier optics. The official syllabus illustrates Fourier analysis and analog filtering with two applications, modulation and Fourier optics in astronomy. Both are direct consequences of the properties above.

Example 3.6 — Amplitude modulation and coherent demodulation A voice signal x of bandwidth B = 4 kHz cannot be radiated efficiently as is. Multiply it by a carrier at f_c = 1 MHz: y(t) = x(t)\cos(2\pi f_c t). Since \cos(2\pi f_c t) = \tfrac12(e^{i2\pi f_c t} + e^{-i2\pi f_c t}), the modulation property gives \hat y(\nu) = \tfrac12\big[\hat x(\nu - f_c) + \hat x(\nu + f_c)\big]: the spectrum is copied around \pm f_c and occupies [996, 1004] kHz. Stations with carriers at least 2B = 8 kHz apart do not overlap (frequency-division multiplexing).

To demodulate, multiply again by the carrier: y(t)\cos(2\pi f_c t) = \tfrac12 x(t) + \tfrac12 x(t)\cos(4\pi f_c t), whose second term lives around \pm 2f_c; a low-pass filter with cutoff between B and 2f_c - B keeps \tfrac12 x. The receiver is thus a product followed by a filter, two operations that the Fourier transform makes transparent.

Example 3.7 — Fourier optics: a telescope is a low-pass filter In the far field, or in the focal plane of a lens of focal length f, the complex amplitude of light diffracted by an aperture with transmission p(x, y) is the two-dimensional Fourier transform of p, evaluated at the spatial frequencies (u, v) = (x'/\lambda f,\ y'/\lambda f) (Fraunhofer diffraction); the observed intensity is its squared modulus.

Slit. A slit of width a is a rectangle, so by the worked example above the intensity is proportional to \mathrm{sinc}^2(a x'/\lambda f), with first zeros at x' = \pm\lambda f/a: the narrower the slit, the wider the pattern, as in the uncertainty principle.

Telescope. A circular pupil of diameter D gives the Airy pattern, whose first dark ring is at the angle \theta \approx 1.22\,\lambda/D. A star, a point source, is imaged as this pattern, the point spread function, and the image of an extended object is the convolution of the object with it: the telescope is a two-dimensional linear time-invariant system (chapter 1).

Its transfer function, the Fourier transform of the point spread function, is the autocorrelation of the pupil and vanishes beyond the spatial frequency D/\lambda: no detail finer than about \lambda/D survives. For D = 8 m and \lambda = 500 nm, \theta \approx 1.22 \times 500 \times 10^{-9}/8 \approx 7.6 \times 10^{-8} \ \text{rad}, about 0.016 arcsecond.

Radio interferometers push the same logic further: each pair of antennas separated by a baseline b measures one Fourier coefficient of the sky image at spatial frequency b/\lambda, and the image is reconstructed by an inverse Fourier transform from these samples.

4. Sampling: the Shannon theorem

A computer stores only the values of a signal at regularly spaced instants. This chapter answers the question every acquisition chain must face: when do these samples determine the continuous signal exactly, and what goes wrong — aliasing — when they do not?

Theorem 4.1 — Shannon-Nyquist (proved) Let x be band-limited: \hat x(\nu) = 0 for |\nu| \geq B, and f_s = 1/T_e > 2B. Then the samples x_n = x(nT_e) suffice to reconstruct exactly x:

Key result — Shannon interpolation x(t) = \sum_{n \in \mathbb{Z}} x_n\, \mathrm{sinc}\Big(\frac{t}{T_e} - n\Big), \qquad \mathrm{sinc}(u) = \frac{\sin(\pi u)}{\pi u}.

Proof.

  • Step 1 (sampling = periodizing the spectrum). Ideal sampling is written x_s(t) = x(t)\cdot\text{Ш}_{T_e}(t) = \sum_n x_n\,\delta(t - nT_e). Product in time ⇒ convolution in frequency (prop. 3.2): \hat x_s(\nu) = \hat x * \widehat{\text{Ш}_{T_e}}(\nu) = \frac{1}{T_e}\sum_{k \in \mathbb{Z}} \hat x\Big(\nu - \frac{k}{T_e}\Big). The spectrum of x_s is the periodization of the spectrum of x, of period f_s.
  • Step 2 (no overlap). If f_s > 2B, the copies (centered at 0, \pm f_s, \pm 2f_s, \dots, each of half-width B) are disjoint: in the central band |\nu| < f_s/2 we read \hat x = T_e\,\hat x_s.
  • Step 3 (ideal low-pass filtering). Let us multiply by R(\nu) = T_e\,\mathbf{1}_{|\nu| < f_s/2}: \hat x = R\,\hat x_s, that is, x = r * x_s where r = \check R is the inverse transform of the rectangle — computed by direct integration: r(t) = T_e\int_{-f_s/2}^{f_s/2} e^{i2\pi\nu t}d\nu = \frac{\sin(\pi f_s t)}{\pi f_s t}\cdot f_s T_e = \mathrm{sinc}(f_s t).
  • Step 4 (interpolation formula). Convolving sinc with the comb of samples gives, since \mathrm{sinc}\big(f_s(t - nT_e)\big) = \mathrm{sinc}(f_st - n) and since the sinc equals 1 at 0 and 0 at the other integers (the samples do not “talk” to each other): x(t) = \sum_n x_n\, \mathrm{sinc}(f_s t - n).

If f_s \leq 2B, the copies overlap: the information is lost without return (two different signals have the same samples).∎

Proposition 4.2 — Aliasing, in formulas After sampling at f_s, every frequency \nu is seen at the apparent frequency \nu_{\text{app}} = \min_{k \in \mathbb{Z}}\ \big|\nu - k f_s\big| \ \in\ \Big[0, \frac{f_s}{2}\Big], the image of \nu in the fundamental band [-f_s/2, f_s/2] modulo f_s.

Example: f_s = 80 Hz, \nu = 60 Hz ⇒ \nu_{\text{app}} = |60 - 80| = 20 Hz (the line “moves backward”); a rising chirp passes f_s/2, then goes back down — the ghost trace of the spectrogram of exercise 3. Wheels that seem to turn backward in the movies, moiré of images: the same phenomenon.

Bandpass sampling For a signal modulating only a narrow band [f_c - B, f_c + B] (radio, sensors), the useful condition bears on the width 2B and not on the highest frequency: it suffices to choose f_s such that the copies of the band do not overlap.

Example: band [19, 21] MHz (f_c = 20, B = 1), f_s = 4.2 MHz, so that 21 = 5 f_s.

Every frequency f of the band folds back to 21 - f \in [0, 2] MHz, inside the fundamental band [0, 2.1] MHz and without collision — the spectrum is transposed to baseband, here flipped (the spectral reversal depends on the parity of the multiple of f_s).

Beware that f_s = 4 MHz, although equal to 2 \times the width, would fail: 20 = 5 \times 4 sits in the middle of the band, and the band and its mirror image fold onto the same interval [0, 1] MHz. This is the basis of software receivers (SDR).

Example 4.3 — Acquisition chain and practical reconstruction Chain: analog signal → analog anti-aliasing filter (cuts at f_s/2, indispensable before sampling) → sampling f_s → quantization (b bits: quantization noise at \approx 6b dB of dynamic range) → digital. Audio CD: B = 20 kHz ⇒ f_s = 44.1 kHz, 16 bits.

The Shannon sinc is perfectly localized in frequency but infinitely long in time: one truncates and windows it (sinc × Hann over M points) — hence a slight aliasing and a soft transition, quantifiable with the tools of ch. 5; linear interpolation (triangle) is its cheap version.

import numpy as np
fs = 200                       # sampling frequency
t  = np.arange(0, 1, 1/fs)
f0 = 60                        # signal 60 Hz
x  = np.sin(2*np.pi*f0*t)

X = np.fft.rfft(x)             # spectrum (positive frequencies)
freqs = np.fft.rfftfreq(len(x), 1/fs)
peak = freqs[np.argmax(np.abs(X))]
print("dominant frequency:", peak, "Hz")

# Aliasing: sampling a 60 Hz signal at 80 Hz "sees" it at 20 Hz
fs2 = 80; t2 = np.arange(0, 1, 1/fs2)
x2 = np.sin(2*np.pi*f0*t2)
print("seen at:", np.fft.rfftfreq(len(x2),1/fs2)[np.argmax(np.abs(np.fft.rfft(x2)))], "Hz")

# Shannon reconstruction on a fine grid (sum over the 200 available samples;
# the error is largest near the ends, where the missing samples are felt)
tr = np.arange(0, 1, 1/2000); xr = np.zeros_like(tr)
for n in range(len(x)):
    xr += x[n]*np.sinc(fs*tr - n)
print("max error after reconstruction:", np.abs(np.sin(2*np.pi*f0*tr) - xr).max())

5. DFT, FFT and digital filtering

Once a signal is a finite list of N samples, its spectrum must be computed by a finite sum, and filters become finite recursions. This chapter defines the discrete Fourier transform, shows how the FFT computes it in about N \log_2 N operations instead of N^2, and explains how to design a digital filter with a prescribed frequency response.

Definition 5.1 — The DFT For x_0, \dots, x_{N-1}:

Key result — the DFT X_k = \sum_{n=0}^{N-1} x_n\, e^{-i 2\pi k n / N}, \qquad k = 0, \dots, N-1.

The DFT samples the transform of the implicitly periodized signal (k \leftrightarrow \nu_k = k f_s/N). Two links to know:

  • DFT ↔︎ Fourier series: if x is T_0-periodic sampled with N points per period and f_s > 2B: X_k = N\,c_k \ \ (\text{for } k < N/2), \qquad X_{N-k} = N c_{-k};
  • DFT ↔︎ Fourier transform: for a segment of duration T = NT_e of an aperiodic signal: \hat x(\nu_k) \approx T_e\,X_k (the so-called “rectangular” approximation, with the spectral leakage of prop. 5.4).

Proposition 5.2 — Circular convolution (proved) Let x \circledast h be the circular convolution (indices modulo N): y_n = \sum_{m=0}^{N-1} x_m\,h_{(n-m) \bmod N}. Then \mathrm{DFT}(x \circledast h)_k = X_k\,H_k.

Proof. Y_k = \sum_n \sum_m x_m h_{(n-m) \bmod N}\,\omega^{kn} \qquad\text{with } \omega = e^{-2\pi i/N}. For m fixed, n \mapsto p = (n-m) \bmod N is a bijection of \mathbb{Z}/N\mathbb{Z} and \omega^{kn} = \omega^{k(m+p)} (since \omega^{kN} = 1). Hence Y_k = \sum_m x_m\omega^{km}\sum_p h_p\omega^{kp} = X_k H_k. The linear convolution is obtained by zero-padding both signals to N' \geq n + m - 1 points: the “circular” wrap-around then no longer has any effect (the overflowing indices fall on zeros).∎

Theorem 5.3 — Radix-2 FFT (proved: the butterfly recurrence) If N = 2^q, the DFT is computed in T(N) = \mathcal{O}(N \log N) operations.

Proof. Separate the even and odd sums: X_k = \underbrace{\sum_{m=0}^{N/2-1} x_{2m}\,\omega^{2mk}}_{E_k} +\ \omega^{k}\,\underbrace{\sum_{m=0}^{N/2-1} x_{2m+1}\,\omega^{2mk}}_{O_k}. Now \omega^2 = e^{-4\pi i/N} is an N/2-th root of unity: E and O are two DFTs of size N/2, periodic with period N/2.

Moreover \omega^{k + N/2} = -\omega^k. Hence the butterfly: X_k = E_k + \omega^k O_k, \qquad X_{k+N/2} = E_k - \omega^k O_k, \qquad k = 0,\dots,\tfrac{N}{2}-1, which produces the N outputs from the N inputs in O(N) multiply-adds once E, O are known.

Recurrence: T(N) = 2T(N/2) + cN = cN\log_2 N. For N = 2^{20} (about 24 s of audio at 44.1 kHz): \approx 2\times10^7 units versus 10^{12} for the definition — a factor of 50,000.∎

Worked example — a 4-point DFT, directly and by butterflies Take x = (1, 2, 0, -1), so N = 4 and \omega = e^{-i\pi/2} = -i.

  • Directly: \begin{aligned} X_0 &= 1 + 2 + 0 - 1 = 2, \\ X_1 &= 1 + 2(-i) + 0 + (-1)(-i)^3 = 1 - 2i - i = 1 - 3i, \\ X_2 &= 1 - 2 + 0 + 1 = 0, \\ X_3 &= \overline{X_1} = 1 + 3i \end{aligned} (for a real signal, X_{N-k} = \overline{X_k}).
  • By butterflies (theorem 5.3): the even samples (1, 0) have the 2-point DFT E = (1, 1), the odd samples (2, -1) have O = (1, 3), and \begin{aligned} X_0 &= E_0 + O_0 = 2, & X_1 &= E_1 + \omega O_1 = 1 - 3i, \\ X_2 &= E_0 - O_0 = 0, & X_3 &= E_1 - \omega O_1 = 1 + 3i: \end{aligned} the same four numbers, assembled from two DFTs of size 2.
  • Check by Parseval: \sum_n |x_n|^2 = 6 and \frac1N\sum_k |X_k|^2 = \frac14(4 + 10 + 0 + 10) = 6.
  • Reading: X_2 = 0 says that x contains no component at the alternating frequency f_s/2, which can be checked directly since 1 - 2 + 0 - (-1) = 0.

Proposition 5.4 — Windows: quantified lobes Analyzing a segment of signal = multiplying it by a window w ⇒ convolving the spectrum with \hat w (prop. 3.2): a pure line becomes the profile \hat w. Characteristics (window of length M, bin = f_s/M; side-lobe levels as tabulated by Oppenheim & Schafer, decay rates from Harris 1978):

Window Main lobe 1st side lobe Decay
rectangular ±1 bin −13 dB −6 dB/octave
Hann ±2 bins −31 dB −18 dB/octave
Hamming ±2 bins −41 dB −6 dB/octave
Blackman ±3 bins −57 dB −18 dB/octave

Universal trade-off: low side lobes ⇒ a wide main lobe (resolution lost). The rectangular one “leaks” at −13 dB: a strong line masks a weak neighboring line at −60 dB; Hann is the default choice.

The Z-transform and the transfer function The Z-transform of a discrete signal (x_n) is the power series X(z) = \sum_{n \in \mathbb{Z}} x_n z^{-n}, defined on its region of convergence, an annulus of the complex plane. On the unit circle z = e^{i2\pi\nu} it reduces to the DTFT, when the circle lies in that region.

Three properties make it the algebra of digital filters: it is linear; a delay by one sample multiplies it by z^{-1}; and a convolution becomes a product, Y(z) = H(z)X(z), where H(z) = \sum_n h_n z^{-n} is the transfer function of the filter.

A recursion y_n + \sum_{j \geq 1} a_j y_{n-j} = \sum_k b_k x_{n-k} therefore becomes A(z)Y(z) = B(z)X(z), so H(z) = B(z)/A(z): a rational function whose poles are the roots of A.

A filter is causal when its region of convergence is the outside of a circle, |z| > r, and then stable exactly when r < 1, that is when all its poles lie inside the unit disk.

Worked example — a first-order recursive filter Take y_n = a\,y_{n-1} + x_n. The delay rule gives Y(z)(1 - a z^{-1}) = X(z), so H(z) = 1/(1 - az^{-1}), with a single pole at z = a.

Expanding the geometric series for |z| > |a| gives H(z) = \sum_{n \geq 0} a^n z^{-n}, so the impulse response is h_n = a^n for n \geq 0: it is summable, hence the filter is stable, if and only if |a| < 1, as theorem 5.5 states.

On the unit circle, |H(\nu)|^2 = 1/(1 - 2a\cos 2\pi\nu + a^2); for a = 0.9 the gain is 1/(1 - 0.9)^2 = 100 (20 dB) at \nu = 0 and 1/1.9^2 \approx 0.28 at \nu = 1/2: an exponential smoother, that is, a low-pass filter. Its output driven by white noise is exactly the AR(1) process of chapter 7.

Digital filters: FIR and IIR FIR (finite impulse response): y_n = \sum_{k=0}^{M} b_k x_{n-k} — always stable (\sum_k |b_k| < \infty: th. 1.3), linear phase if the b_k are symmetric, designed by windowing.

IIR (infinite): y_n = \sum b_k x_{n-k} - \sum_{j\geq1} a_j y_{n-j}, transfer function H(z) = B(z)/A(z) — more selective at equal coefficient counts, but stability must be checked.

Theorem 5.5 — Stability of an IIR: the poles (proved) For a rational IIR, BIBO stable ⇔ all the poles of H(z) lie in the open disk |z| < 1.

Proof. After cancelling the common factors of B and A, a partial-fraction expansion of H(z) writes the causal impulse response (beyond finitely many initial terms) as a combination of modes n^k p_j^n, with k smaller than the multiplicity of the pole p_j.

Each mode is summable if and only if |p_j| < 1, since then \sum_n n^k |p_j|^n < \infty. A pole on the unit circle gives a mode that does not decay (bounded if the pole is simple, growing like n^k if it is multiple); a pole outside gives a divergent exponential. The accumulator h = \mathbf{1}, with a simple pole at 1, illustrates the unit-circle case.∎

Example 5.6 — Designing a low-pass FIR by windowing, step by step Objective: H_d(\nu) = \mathbf{1}_{|\nu| \leq \nu_c} (ν normalized by f_s).

  1. Ideal response in time: h_d[n] = 2\nu_c\,\mathrm{sinc}(2\nu_c n) (inverse FT of the rectangle, cf. proof of th. 4.1) — infinite and non-causal.
  2. Delay of M/2 for causality (and symmetry ⇒ linear phase).
  3. Windowing: h[n] = h_d[n - M/2]\,w[n - M/2], n = 0,\dots,M.
  4. Analysis: H = H_d * \hat w; the transition band equals the width of the main lobe of \hat w — rule of thumb M \cdot \Delta\nu \approx 4 (Hann), \approx 6 (Blackman). The ripple, in the pass band and in the stop band alike, comes from the integrated side lobes and sits well below the peak side lobe of prop. 5.4: about −21 dB for the rectangular window (the 9% Gibbs overshoot of example 2.4), −44 dB for Hann, −53 dB for Hamming and −74 dB for Blackman (Oppenheim & Schafer).

The effective frequency response is the DTFT of h, H(\nu) = \sum_n h[n]e^{-2\pi i\nu n}, evaluated by zero-padding the FFT (trick: a DFT with N_{\text{fft}} \gg M points interpolates the DTFT).

A finer alternative: frequency-domain least squares — minimize \sum_{\omega \in \Omega} |H_b(\omega) - D(\omega)|^2 over a frequency grid, a quadratic problem in b, hence a linear system (Parks-McClellan equiripple method: the same spirit in infinite norm).

For IIRs, one projects an analog prototype via the bilinear transform s \mapsto \frac{2}{T_e}\frac{1 - z^{-1}}{1 + z^{-1}} (Butterworth, Chebyshev…) while keeping the poles inside the disk — the “smooth” filter of example 5.7.

Example 5.7 — Low-pass denoising

rng = np.random.default_rng(0)
sig = x + 0.5*rng.standard_normal(len(x))
SIG = np.fft.rfft(sig)
SIG[np.abs(freqs) > 80] = 0            # abrupt cutoff at 80 Hz (zero-phase)
clean = np.fft.irfft(SIG, n=len(sig))
# Abrupt cutoff in frequency = Gibbs in time (example 2.4): prefer a soft
# profile (4th-order Butterworth via scipy.signal.butter + filtfilt, or a Hann mask)

Images: two-dimensional convolution and DFT A grey-level image is a two-dimensional discrete signal x[m, n]. A linear shift-invariant operation on it is a 2D convolution y[m, n] = \sum_{k,l} h[k, l]\,x[m-k, n-l] with a small kernel h, and the 2D DFT X[k, l] = \sum_{m,n} x[m, n]\,e^{-i2\pi(km/M + ln/N)} turns it into a product, exactly as in one dimension (proposition 5.2).

The 2D DFT is separable: it is computed by 1D FFTs along the rows, then along the columns. Low spatial frequencies carry the smooth regions, high ones the edges and textures, and sampling an image too coarsely aliases fine periodic patterns into moiré (proposition 4.2).

Example 5.8 — Blur and edge detection, read in frequency The 3 \times 3 box blur averages each pixel with its 8 neighbours. Its kernel is separable, the product of (\tfrac13, \tfrac13, \tfrac13) along each axis, so its frequency response is H(\nu_1, \nu_2) = g(\nu_1)g(\nu_2) \quad\text{with}\quad g(\nu) = (1 + 2\cos 2\pi\nu)/3. It equals 1 at \nu = 0 (flat regions pass unchanged), vanishes at \nu = 1/3, but equals -1/3 at \nu = 1/2: the finest checkerboard is attenuated but survives with inverted sign, an artefact that a Gaussian kernel avoids because its transform stays positive (proposition 3.3).

The discrete Laplacian kernel, with -4 at the centre and 1 at the four neighbours, has response H(\nu_1, \nu_2) = 2\cos 2\pi\nu_1 + 2\cos 2\pi\nu_2 - 4 = -4\big(\sin^2\pi\nu_1 + \sin^2\pi\nu_2\big). It is 0 at the origin, so constant regions are erased, and its magnitude is largest, 8, at (\tfrac12, \tfrac12): a high-pass filter that responds to edges. The first layer of a convolutional network (a preview, defined in the box that opens ch. 8; taught in Deep Learning, ch. 4, or in its P2 variant, ch. 1) learns kernels of exactly this kind.

import numpy as np
img = np.random.default_rng(0).random((64, 64))
nu1, nu2 = np.meshgrid(np.fft.fftfreq(64), np.fft.fftfreq(64), indexing="ij")
lap = -4 * (np.sin(np.pi * nu1)**2 + np.sin(np.pi * nu2)**2)   # Laplacian frequency response
edges = np.real(np.fft.ifft2(np.fft.fft2(img) * lap))           # circular convolution with the kernel

6. Time-frequency analysis and wavelets

The Fourier transform of a whole recording says which frequencies are present, but not when: a melody and the same melody played backward have spectra of the same magnitude. This chapter localizes the analysis in time, first with the short-time Fourier transform (one fixed window), then with wavelets (windows whose length adapts to the frequency).

Definition 6.1 — STFT and spectrogram \mathrm{STFT}_x(t, \nu) = \int x(\tau)\, w(\tau - t)\, e^{-i2\pi\nu\tau}\, d\tau, with w a short window (Hann). The spectrogram |\mathrm{STFT}_x|^2 describes the energy across time and frequency: it is the image of sound that a convolutional network “sees” (ch. 8).

Proposition 6.2 — Exact resolution of the STFT At each instant, the STFT is the FT of the windowed segment x\cdot w(\cdot - t): its frequency resolution is the width of the main lobe of \hat w (prop. 5.4); its time resolution is the duration T of the window. For Hann: main lobe \pm 2/T in frequency, hence \Delta t = T, \qquad \Delta\nu \approx \frac{2}{T} \ \ (\text{Hann}), \qquad \Delta t\,\Delta\nu \approx 2, to be compared with the floor 1/(4\pi) \approx 0.08 (standard deviations) of th. 3.4: the STFT works at a resolution fixed in advance, far from the floor, but with time-frequency cells all identical.

Long window ⇒ good frequency resolution, poor time resolution — and conversely: the trade-off is irreversible for a given window.

Gabor frames: sampling the STFT A continuously stored STFT is redundant; one samples it on a grid (t,\nu) = (m\,a,\ n\,b). The family of translated-modulated windows w_{m,n} is a Gabor frame when it allows a stable reconstruction, by a dual window. Results (admitted):

  • the density condition ab \leq 1 is necessary, but it is not sufficient for a general window;
  • for the Gaussian window, ab < 1 suffices (Lyubarskii, 1992; Seip, 1992);
  • for a window that is positive on an interval of length T and zero outside, a < T and b \leq 1/T suffice: this “painless” case is that of every practical spectrogram, with a hop shorter than the window and an FFT at least as long as the window;
  • at the critical density ab = 1, no such frame that is also a basis has a window well localized in both time and frequency (the Balian-Low theorem): one must accept redundancy, or change the geometry (wavelets below).

This is exactly the structure of spectrograms: frames every 10 ms, N frequencies — a frame, not a basis.

STFT inversion by overlap-add (proved) In discrete time, frame m of the STFT multiplies the signal by a window w of length L, shifted by m hops of R samples, and takes a DFT on N \geq L points: X_m[k] = \sum_{j=0}^{L-1} x[mR + j]\,w[j]\,e^{-i2\pi kj/N}. If the shifted windows add up to a constant, \sum_{m} w[n - mR] = C \quad\text{for every } n \qquad (\text{COLA: constant overlap-add}), then inverting the DFT of every frame and adding the frames at their positions gives back C\,x exactly.

Proof. Since N \geq L, the inverse DFT of frame m returns the windowed segment y_m[j] = x[mR + j]\,w[j] with no circular wrap-around (prop. 5.2). Placing each segment at time mR and summing, \sum_m y_m[n - mR] = x[n]\sum_m w[n - mR] = C\,x[n]. ∎

The periodic Hann window satisfies COLA with C = 1 at the hop R = L/2, because w[n] + w[n + L/2] = 1, and with C = 2 at R = L/4.

Other hops, such as the 25 ms windows every 10 ms of ch. 8 (R = 0.4L), are not COLA. Libraries then use the weighted overlap-add of Griffin and Lim (1984): each inverted frame is multiplied by w once more and the sum is divided by \sum_m w^2[n - mR]. The result is exact as long as that sum never vanishes (the NOLA condition checked by scipy.signal.istft). When the STFT has been modified, by a mask for instance, it returns the signal whose STFT is closest in least squares.

Every mask-based denoiser and source separator (example 8.2, exercise 12) relies on this inversion.

Multi-resolution STFT losses A neural vocoder must output a waveform \hat y whose spectrogram matches that of the target y. One window cannot be sharp in both time and frequency (prop. 6.2), so the training loss sums spectral distances over J STFT configurations S_1, \dots, S_{J} (for instance windows of 240, 600 and 1200 samples): \mathcal{L}_{\mathrm{MR}}(y, \hat y) = \frac{1}{J}\sum_{r=1}^{J}\left(\frac{\big\|\,|S_r y| - |S_r \hat y|\,\big\|_F}{\big\|\,|S_r y|\,\big\|_F} + \frac{1}{K_r}\big\|\log|S_r y| - \log|S_r \hat y|\big\|_1\right), with K_r the number of time-frequency cells of S_r. The first term matches the strong components, the second the low-energy ones that the logarithm brings up.

Only magnitudes enter: the network remains free to choose a consistent phase. This loss trains Parallel WaveGAN (Yamamoto, Song and Kim, 2020) and many neural audio codecs; HiFi-GAN (Speech Technology, ch. 6) uses a single-resolution Mel version of it next to its discriminators.

Definition 6.3 — Wavelets: analysis at variable scale Mother wavelet \psi (zero mean, localized); \psi_{a,b}(t) = \frac1{\sqrt a}\psi\big(\frac{t-b}{a}\big); the CWT samples (a,b), often in octaves a = 2^{-j} (scalogram, Q = \Delta\nu/\nu constant): good time localization at high frequency, frequency localization at low frequency — suited to sounds (closely spaced harmonics at the bottom of the spectrum, brief transients at the top), where the STFT is uniform.

Example 6.4 — Haar, a complete example down to reconstruction The Haar wavelet: \varphi = \mathbf{1}_{[0,1)} (scale), \psi = \mathbf{1}_{[0,1/2)} - \mathbf{1}_{[1/2,1)} (detail). Signal x = (9, 7, 3, 5):

  • level 1: averages a_1 = (8, 4), details d_1 = (1, -1);
  • level 2: a_2 = (6), d_2 = (2);
  • coding: x \leftrightarrow (a_2, d_2, d_1) = (6, 2, 1, -1);
  • reconstruction: 9 = 6+2+1, 7 = 6+2-1, 3 = 6-2+1, 5 = 6-2-1.

The multiresolution equations: \varphi(t) = \sqrt2\sum_k h_k\varphi(2t-k), \qquad \psi(t) = \sqrt2\sum_k g_k\varphi(2t-k); for Haar, h = (\tfrac{1}{\sqrt2}, \tfrac{1}{\sqrt2}), g = (\tfrac{1}{\sqrt2}, -\tfrac{1}{\sqrt2}) — h = averaging filter (low-pass), g = differentiator (high-pass).

The Daubechies D_N compactify this scheme with 2N coefficients and N vanishing moments: on a piecewise regular signal, the detail coefficients decay like 2^{-Nj} — this is what makes the coding sparse (a polynomial of degree < N has exactly zero details: Haar = D_1 annihilates constants).

DWT = a filter bank, and the JPEG 2000 link The DWT is implemented without any integral: a cascade of low-pass (h) / high-pass (g) filterings followed by decimation by a factor 2, iterated on the low branch; the “quadrature mirror” conditions on (h,g) cancel the aliasing introduced by the decimations. JPEG 2000 codes the subbands of a biorthogonal DWT (5/3 lossless, 9/7 lossy): an image = a 2D signal, and compression exploits lemma 2.3 (regularity ⇒ small coefficients). Same bank, from denoising to compression: the theory does not change.

Definition 6.5 — The discrete cosine transform (DCT-II) For x_0, \dots, x_{N-1}, X_k = s_k\sum_{n=0}^{N-1} x_n \cos\Big(\frac{\pi}{N}\big(n + \tfrac12\big)k\Big), \qquad k = 0, \dots, N-1, with s_0 = \sqrt{1/N} and s_k = \sqrt{2/N} otherwise, which makes the transform orthonormal. Up to a phase factor, it is the DFT of the even, mirror-symmetric extension of x to length 2N.

The DFT implicitly repeats the block periodically and so creates a jump between x_{N-1} and x_0; the mirror extension has no jump. By lemma 2.3 (regularity implies fast decay of the coefficients), the DCT coefficients of a smooth block decay much faster: the energy is compacted into few coefficients.

For strongly correlated signals, such as an AR(1) with correlation close to 1, the DCT is close to the Karhunen–Loève transform, the PCA of the block.

Example 6.6 — Energy compaction, then JPEG in four steps For the ramp x = (1, 2, \dots, 8), the orthonormal DCT is approximately (12.73,\ -6.44,\ 0,\ -0.67,\ 0,\ -0.20,\ 0,\ -0.05): the first two coefficients hold 99.8% of the energy.

The unitary DFT of the same block has moduli 12.73, 3.70, 2.00, 1.53, 1.41, \dots, and its three lowest frequencies hold only 92.8% of the energy, because of the jump from 8 back to 1.

The JPEG standard exploits this.

  1. The image is cut into 8 \times 8 blocks.
  2. Each block is transformed by the 2D DCT, which is separable: a 1D DCT on the rows, then on the columns.
  3. Each coefficient X[k, l] is quantized to \mathrm{round}(X[k, l]/Q[k, l]), with steps Q that grow with frequency, because the eye is less sensitive to fine detail; most high-frequency coefficients become 0. This is the only lossy step.
  4. The coefficients are read in zig-zag order from low to high frequency, so that the zeros form long runs, which are coded by run lengths and Huffman codes.

At high compression the blocks become visible at their borders; the wavelet transform of JPEG 2000 (note above) avoids this blocking.

7. Random signals and spectral estimation

Noise, speech and sensor readings are never the same twice: each recording is one realization of a random process, and only averages are reproducible. This chapter defines stationary processes, their autocorrelation and their power spectral density, and shows how to estimate that density from a finite recording without being fooled by the variance of the periodogram.

Definition 7.1 — Process, wide-sense stationarity (WSS) A random signal is a family X(t, \omega) of random variables. It is wide-sense stationary if its mean is constant m = \mathbb{E}[X(t)] and if its autocorrelation R_X(t+\tau, t) = \mathbb{E}\big[(X(t+\tau)-m)\,\overline{(X(t)-m)}\big] depends only on the lag \tau (written R_X(\tau)). The conjugate sits on the earlier time, so that the pure frequency e^{i2\pi\nu_0 t} (with a random phase) has R_X(\tau) = e^{i2\pi\nu_0\tau} and its power appears at +\nu_0 in theorem 7.3; for real signals the order does not matter. This is the right framework for noise, speech, sensors: each realization is a signal of this course, but the interesting quantities are expectations.

Proposition 7.2 — Properties of the autocorrelation R_X(0) = \mathbb{E}|X - m|^2 is the average power; |R_X(\tau)| \leq R_X(0) \qquad \big(\text{Cauchy-Schwarz: } |\mathbb{E}[X_\tau\overline{X_0}]| \leq \sqrt{R(0)R(0)}\big); R_X is hermitian (R_X(-\tau) = \overline{R_X(\tau)}) and positive definite. “Noise” is decorrelated at a distance: R(\tau) \to 0.

Theorem 7.3 — Wiener-Khinchin (statement) and the spectral density For a WSS process, there exists a positive measure S_X, the power spectral density, such that R_X(\tau) = \int e^{i2\pi\nu\tau}\, dS_X(\nu), \qquad\text{and if } R_X \in L^1:\quad S_X(\nu) = \widehat{R_X}(\nu). The power is distributed across the spectrum: \mathbb{E}|X-m|^2 = R_X(0) = \int dS_X: Parseval (th. 2.2) becomes a definition.

(The proof — Bochner’s theorem — goes beyond the scope of the course; the intuition: a WSS process is an incoherent superposition of sinusoids, each carrying its share of the power.)

  • White noise: decorrelated to all orders, R = \sigma^2\delta ⇒ S = \sigma^2 constant — “all frequencies on equal footing”. A band-limited white noise (physically realizable): S = \sigma^2\mathbf{1}_{|\nu| \leq B} ⇒ R(\tau) = 2B\sigma^2\,\mathrm{sinc}(2B\tau): decorrelation at the scale 1/B (direct computation of the FT of the rectangle, cf. th. 4.1).
  • Filtered noise: Y = h * X ⇒ S_Y = |H|^2\,S_X (colored noise): the filter shapes the spectrum squared.
  • AR(1) (auto-regressive model X_n = \rho X_{n-1} + W_n, |\rho| < 1): R[\tau] = \frac{\sigma_W^2}{1-\rho^2}\,\rho^{|\tau|} and S(\nu) = \frac{\sigma_W^2}{|1 - \rho e^{-i2\pi\nu}|^2} — the simplest model of colored noise, ubiquitous in finance and sensors.

Worked example — the spectrum of an AR(1) noise, in numbers Take \rho = 0.9 and \sigma_W^2 = 1. The power is R[0] = 1/(1 - 0.81) \approx 5.26, and the correlation between neighbours is R[1]/R[0] = \rho = 0.9: consecutive samples are strongly alike.

The density at zero frequency is S(0) = 1/|1 - 0.9|^2 = 100, and at the highest frequency \nu = \tfrac12 (where e^{-i\pi} = -1) it is S(\tfrac12) = 1/|1 + 0.9|^2 \approx 0.28: a ratio of about 360, that is 25.6 dB. The process is a low-pass (“red”) noise; with \rho = -0.9 the two values swap and the noise becomes high-pass.

A useful check: integrating S over one period of frequencies, \int_{-1/2}^{1/2} S(\nu)\,d\nu, gives back R[0] \approx 5.26, which is Wiener-Khinchin at \tau = 0.

Definition 7.4 — The periodogram “Natural” estimator of S_X from N samples: \hat S_N(\nu) = \frac{1}{N}\Big|\sum_{n<N} X_n\, e^{-i2\pi\nu n}\Big|^2 \qquad (\text{on the DFT grid: } |X_k|^2/N).

Theorem 7.5 — Bias and variance of the periodogram

  • Bias (leakage): for a centred process, \mathbb{E}[\hat S_N] = S * F_N, \qquad F_N(\nu) = \frac{1}{N}\left|\frac{\sin(\pi N\nu)}{\sin(\pi\nu)}\right|^2, the Fejér kernel of Definition 0.2 (at \theta = 2\pi\nu): non-negative, of mass 1 over one period, with a main lobe of ±1 bin. The expectation “smears” the spectrum over this lobe and its side lobes — a strong line masks a weak line (hence the windows, prop. 5.4). The bias tends to 0 as N \to \infty (F_N \to \delta).
  • Variance: for a Gaussian process, \mathrm{Var}[\hat S_N(\nu)] \approx S_X(\nu)^2 for every N — the relative error is 100% and does not decrease: the periodogram alone is a poor estimator, whatever the observation duration.

Proof. Sketch. Bias: expanding the squared modulus, with N - |\tau| pairs (n, m) at lag n - m = \tau, \mathbb{E}[\hat S_N(\nu)] = \frac{1}{N}\sum_{n,m<N} R[n-m]\,e^{-i2\pi\nu(n-m)} = \sum_{|\tau|<N}\Big(1 - \frac{|\tau|}{N}\Big) R[\tau]\,e^{-i2\pi\nu\tau}. This is the transform of R multiplied by a triangular (Bartlett) window. The transform of that window is F_N = |\hat w_N|^2/N, with \hat w_N the transform of the rectangular window of length N, and a product in time is a convolution in frequency (prop. 3.2). Since F_N \geq 0, the expected power is never negative, as it must be.

Variance: on the DFT grid, the \hat S_N(\nu_k) are asymptotically decorrelated and, for a Gaussian, \hat S_N(\nu) approximately follows S(\nu)\cdot \chi^2_2/2 (squared modulus of two independent Gaussians). Now \mathrm{Var}[\chi^2_2/2] = 1, hence a variance of order S^2, independent of N.∎

Bartlett and Welch methods — reducing the variance by averaging

  • Bartlett: split the N points into L disjoint segments of M = N/L points and average the periodograms: \mathrm{Var} \approx S^2/L…

    …at the price of an increased bias: the windows are L times shorter than the record, so the main lobe is L times wider and the resolution drops from f_s/N to f_s/M — the central bias/variance trade-off of spectral estimation.

  • Welch: overlapping segments (50%) and Hann windows: the same variance gain (Hann segments overlap “for free”), leakage controlled at −31 dB; this is the default method of instrumentation.

from scipy.signal import welch
n  = 20000; fs = 500.0
t  = np.arange(n)/fs
x  = np.sin(2*np.pi*50*t) + 0.02*np.sin(2*np.pi*80*t) + 0.5*rng.standard_normal(n)

f, S = welch(x, fs=fs, nperseg=1024, window='hann', noverlap=512)
# L = 2*20000/1024 - 1 segments: the line at 80 Hz (-34 dB below the line at 50 Hz)
# emerges from the noise, whereas the raw periodogram drowns it (th. 7.5)

Linear prediction and the Yule–Walker equations Let X be a centred WSS process with autocorrelation R. The best linear predictor of X_n from its p previous values, \hat X_n = \sum_{k=1}^{p} a_k X_{n-k}, minimizes \mathbb{E}|X_n - \hat X_n|^2.

Setting the derivatives with respect to the a_j to zero shows that the error must be orthogonal to each past value, \mathbb{E}[(X_n - \hat X_n)X_{n-j}] = 0, which gives the Yule–Walker equations \sum_{k=1}^{p} a_k\,R[j - k] = R[j], \qquad j = 1, \dots, p, a Toeplitz linear system, with prediction-error variance \sigma^2 = R[0] - \sum_k a_k R[k].

If X is an AR(p) process, X_n = \sum_{k=1}^p a_k X_{n-k} + W_n with white noise W, its coefficients solve these equations: the AR model is the optimal predictor, and the innovation W_n is the part of X_n that no linear function of the past can predict.

Replacing R by the empirical autocorrelation fits an AR model to data, and the fitted model gives a parametric spectral estimate \hat S(\nu) = \sigma^2/|1 - \sum_k a_k e^{-i2\pi\nu k}|^2, smooth by construction, an alternative to the averaged periodograms above.

Worked example — fitting an AR(2) by Yule–Walker Suppose the empirical autocorrelations of a recording are \hat R[0] = 2.0, \hat R[1] = 1.2 and \hat R[2] = 0.5. The equations read 2.0\,a_1 + 1.2\,a_2 = 1.2 \qquad\text{and}\qquad 1.2\,a_1 + 2.0\,a_2 = 0.5. The determinant is 4 - 1.44 = 2.56, so a_1 = (2.4 - 0.6)/2.56 \approx 0.703 \qquad\text{and}\qquad a_2 = (1.0 - 1.44)/2.56 \approx -0.172. The predictor is \hat X_n = 0.703\,X_{n-1} - 0.172\,X_{n-2}, and the error variance is 2.0 - 0.703 \times 1.2 + 0.172 \times 0.5 \approx 1.24: prediction removes 38% of the power.

The poles of 1/(1 - 0.703z^{-1} + 0.172z^{-2}) are complex with modulus \sqrt{0.172} \approx 0.41 < 1, so the fitted model is stable. For an AR(1) with \rho = 0.9, the single equation gives back a_1 = R[1]/R[0] = 0.9. Linear prediction of this kind is the basis of classical speech coding (LPC).

Complement to the official syllabus: Wiener filtering The detailed syllabus of the course unit explicitly mentions “filtering and prediction of stationary random signals” and “Wiener filtering”: here is the object. We observe Y = X + W (signal X WSS, noise W WSS, with spectra S_X, S_W, assumed uncorrelated) and we seek the LTI filter h that minimizes the mean squared error \mathbb{E}|X - h * Y|^2. The (non-causal) solution is the Wiener filter:

Key result — the Wiener filter H(\nu) = \frac{S_X(\nu)}{S_X(\nu) + S_W(\nu)}.

Reading: at frequencies where the signal dominates (S_X \gg S_W) the filter transmits (H \approx 1); where the noise dominates, it cancels out — the spectrum decides, optimally in the L^2 sense.

This is exactly the “filtered noise” of Theorem 7.3 (S_Y = |H|^2 S_X) inverted: instead of computing the effect of a filter on a spectrum, we design the filter that best separates signal and noise, from their estimated spectral densities (Definition 7.4 and Theorem 7.5 — estimated in practice by Bartlett/Welch).

The causal filter is built by spectral factorization (via the transfer function H(z) of ch. 5); the linear prediction of an AR process (the X_n = \rho X_{n-1} + W_n of Theorem 7.3) is its “zero noise” case: predicting is filtering toward the future. The Kalman filter, its recursive state-space variant, follows.

The Kalman filter: recursive Gaussian conditioning

The Wiener filter needs a stationary process, and its frequency-domain form above uses the whole recording. A tracker or a navigation system needs an estimate at every instant, updated in constant time when a measurement arrives, under a model that may change over time. The Kalman filter provides it, and each of its updates is the Gaussian conditioning of Refresher in Statistics, theorem 3.2.

Definition 7.6 — Linear-Gaussian state-space model A hidden state x_t \in \mathbb{R}^d evolves linearly and is observed through noisy linear measurements z_t: x_t = F x_{t-1} + w_t, \qquad z_t = H x_t + v_t, with w_t \sim \mathcal{N}(0, Q), v_t \sim \mathcal{N}(0, R) and x_0 \sim \mathcal{N}(\hat x_0, P_0), all independent.

In this section only, H is the observation matrix and R the covariance of the measurement noise — the letters of the tracking literature and of Computer Vision Applications, ch. 5 — not a frequency response or an autocorrelation.

Theorem 7.7 — The Kalman filter (proved) Write \hat x_{t|s} and P_{t|s} for the mean and covariance of x_t given z_1, \dots, z_s. These conditional laws are Gaussian, and each new measurement updates them in two steps.

  • Predict: \hat x_{t|t-1} = F\hat x_{t-1|t-1}, \qquad P_{t|t-1} = F P_{t-1|t-1} F^\top + Q.
  • Update, with the innovation e_t = z_t - H\hat x_{t|t-1} and the gain K_t: K_t = P_{t|t-1}H^\top\big(HP_{t|t-1}H^\top + R\big)^{-1}, \hat x_{t|t} = \hat x_{t|t-1} + K_t e_t, \qquad P_{t|t} = (I - K_t H)P_{t|t-1}.

Proof. By induction on t, everything being conditioned on z_1, \dots, z_{t-1}.

  1. Predict. x_t = Fx_{t-1} + w_t is an affine image of the Gaussian x_{t-1} plus an independent Gaussian, hence Gaussian, with mean F\hat x_{t-1|t-1} and covariance FP_{t-1|t-1}F^\top + Q.
  2. Joint law. Since v_t is independent of x_t, the pair (x_t, z_t) is jointly Gaussian, with \mathrm{Cov}(x_t, z_t) = P_{t|t-1}H^\top and \mathrm{Cov}(z_t) = HP_{t|t-1}H^\top + R.
  3. Condition. Gaussian conditioning, with Y = x_t and X = z_t, gives the mean \hat x_{t|t-1} + K_t e_t and the covariance P_{t|t-1} - K_t H P_{t|t-1}.∎

The gain weighs the prediction against the measurement according to their uncertainties: in one dimension with H = 1, K = P_{t|t-1}/(P_{t|t-1} + R) is close to 1 for a precise sensor and close to 0 for a noisy one.

Worked example — tracking a drifting level Take a scalar random walk observed in noise: F = H = 1, Q = 1 and R = 2. Whatever P_0, the variances converge. At the fixed point, the predicted variance p = P_{t|t-1} solves p = \frac{pR}{p + R} + Q \iff p^2 - p - 2 = 0 \iff p = 2, so the gain is K = p/(p + R) = 0.5 and the filtered variance is P_{t|t} = 1.

The steady-state filter is therefore \hat x_{t|t} = 0.5\,\hat x_{t-1|t-1} + 0.5\,z_t: the first-order recursive filter of chapter 5 with a = 0.5 (and an input gain of 0.5, so that constants pass unchanged), an exponential smoother whose memory is set by the ratio of the two noise variances. A noisier sensor (larger R) gives a smaller gain and a longer memory.

Wiener, Kalman and the trackers downstream When the model is stationary, the gain converges and the steady-state Kalman filter is the causal Wiener filter of the note above, computed recursively.

With a constant-velocity state (position and velocity), the same equations are the motion model of the SORT and DeepSORT trackers of Computer Vision Applications, ch. 5. Nonlinear models are handled by linearizing F and H (the extended Kalman filter) or by particle filters.

Go further — the Y1.2 forecasting baseline In the Year 1 end-of-year project on fleet dispatch (Y1.2), the hourly trip counts of the NYC TLC Trip Record Data are first read with this chapter: Welch reveals the daily and weekly lines (periods 24 h and 168 h), and an AR model fitted by Yule–Walker once these lines are removed is the baseline that the spatio-temporal graph network must beat (solo project 3).

8. The signal for modern ML

Speech recognizers, audio generators and many sensor models do not consume raw samples but a time-frequency representation built with the tools of chapters 5 and 6. This chapter assembles the standard Mel pipeline, explains how generated signals are compared through distances between distributions of spectra, and points out where Fourier analysis hides inside modern architectures.

Preview — the learned models met in this chapter

Students of the AI-systems route reach this chapter before any deep-learning course. Two lines per model suffice here; each is taught later.

  • Convolutional network (CNN): a stack of learned convolutions (ch. 1) and pointwise nonlinearities, so each layer is a bank of FIR filters fitted to data. Taught in Deep Learning, ch. 4 or DL P2, ch. 1.
  • Transformer: a network whose layers mix the frames of a sequence by self-attention, with weights \mathrm{softmax}(QK^\top/\sqrt{d_k}) computed from the frames themselves. Taught in Deep Learning, ch. 6 or DL P2, ch. 0.
  • Autoencoder: an encoder maps a signal to a short code and a decoder maps the code back, both trained to reconstruct the input. Taught in Deep Learning, ch. 9 or, through its variational form, DL P2, section 7.3.
  • Graph convolutional network (GCN): a filter on signals carried by the nodes of a graph, a polynomial of the graph Laplacian taking the place of the convolution. Taught in Graph ML, ch. 4–5.
  • Optimal transport (OT): the cheapest way to move one distribution of mass onto another; its cost is a distance between distributions (Wasserstein). Taught in Computational Optimal Transport, ch. 1–3.

Mel spectrogram and MFCC: the standard pipeline

  1. Frames: 25 ms cuts every 10 ms, Hann window (a Gabor frame, ch. 6); STFT.
  2. Mel scale (perceptive, linear < 1 kHz then logarithmic): m = 2595\,\log_{10}\Big(1 + \frac{f}{700}\Big), \qquad f = 700\,\big(10^{m/2595} - 1\big).
  3. Mel filter bank: 20-80 triangular filters with edges uniformly spaced in Mel — band-pass filters of increasing width in Hz: E_m = \sum_\nu W_m(\nu)\,|\mathrm{STFT}(\nu)|^2 (a sparse weighting of the spectrogram).
  4. Log (amplitude perception); DCT along the Mel axis (a DCT = a “real FT”): the c_j are the MFCC; one keeps 13-20 coefficients + derivatives \Delta, \Delta\Delta, and normalizes in mean/variance per speaker.

Historically features for lightweight models (GMM-HMM), the MFCC remain the compact entry point; deep networks rather consume the raw Mel spectrogram.

Networks on spectrograms The Mel spectrogram is the image that 2D CNNs and audio transformers consume: the first layers of audio CNNs learn oriented frequency detectors — Gabor filters (ch. 6), see the Deep Learning module. The emblematic data augmentation (masking of frequency bands and time frames) operates directly in the time-frequency plane: the spectrogram is treated as a texture to be repaired. Speech technology, transcription, audio-multimodal transformers: this whole pipeline is the signal preprocessing of modern ML.

Distances between spectral distributions: FID and OT spectral loss Comparing two sets of signals = comparing two distributions of spectrograms (or of features). Standard: fit a Gaussian to the features of each set and measure the Fréchet distance d^2\big(\mathcal{N}(\mu_1,\Sigma_1), \mathcal{N}(\mu_2,\Sigma_2)\big) = \|\mu_1 - \mu_2\|^2 + \mathrm{tr}\Big(\Sigma_1 + \Sigma_2 - 2\big(\Sigma_1^{1/2}\Sigma_2\Sigma_1^{1/2}\big)^{1/2}\Big). This is nothing but W_2^2 between Gaussians (the Bures metric): the FID of image generators and the “audio Fréchet distance” of audio generators are Wasserstein distances in disguise.

The OT spectral loss of audio generative models generalizes the idea: comparing real/generated batches by Gaussian W_2 (closed form above, differentiable) or by Sinkhorn on the spectra. The complete proof of the Bures formula and the Sinkhorn algorithm are on the syllabus of the Computational Optimal Transport module (M1P2, ch. 2 and 3, same instructor) — this chapter is its preparation.

In 1D (spectra or one-dimensional signals), everything reduces to the quantile coupling (proved there): W_2^2 = \int_0^1\big(F_1^{-1}(u) - F_2^{-1}(u)\big)^2du.

Example 8.1 — 1D domain adaptation on signals Two sensors measure the same phenomenon with different drifts (offset, gain, frequency response): the distributions of the frames differ, not the semantics. Alignment by 1D OT: transport the quantiles of sensor A onto those of sensor B (monotone transport, computable by sorting); the same principle at d-dimensional scale (feature clouds) is the optimal transport domain adaptation detailed in the OT module (ch. 4), whose instructor is a co-author. The signal provides the spectrum; OT provides the geometry.

Fourier in the architectures: convolutions and positions

  • Convolution of CNNs = the discrete convolution of ch. 1: the learned kernel is a filter whose frequency response (DTFT, ch. 5) explains its behavior — CNNs learn banks of Gabor filters; convolution via FFT in \mathcal{O}(n\log n) when the kernels are large (th. 5.3).
  • Convolution on graphs (preview): the GCN layer H^{(\ell+1)} = \sigma\big(\tilde D^{-1/2}\tilde A \tilde D^{-1/2} H^{(\ell)} W^{(\ell)}\big) is a convolution in the spectral sense. Here \tilde A = A + I is the adjacency matrix with self-loops, \tilde D its degree matrix, and the rows of H^{(\ell)} are the node features (row-major, as in Graph ML), not a frequency response. The Kipf & Welling filter is a degree-1 polynomial of the normalized Laplacian of the graph with self-loops; Graph ML, ch. 4 builds the graph Fourier transform and its polynomial filters, and ch. 5 derives the GCN from them — the later deepening of this bullet.
  • Positional encodings: the sine/cosine of transformers (\mathrm{PE}(p, 2i) = \sin(p/10000^{2i/d})…) = Fourier bases; positional encodings on graphs (LapEig) reuse the eigenvectors of the Laplacian — Fourier again, see the Graph ML module.

Complement to the official syllabus: source separation and the learning dictionary The “Signal representation and learning dictionary” block of the detailed syllabus of the course unit names source separation and the learning dictionary: here is the principle. The idea — officially that of the course — is to move from fixed bases (Fourier, wavelets, ch. 2 and 6) to a learned dictionary: a signal is represented as x \approx D\alpha with D = [d_1, \dots, d_K] a dictionary of waveforms and \alpha sparse (few active coefficients).

Learning D (for instance by non-negative matrix factorization — NMF — or K-SVD-type sparse least squares) on spectrograms of mixtures separates the sources: each component d_k\alpha_k concentrates on one source (a speaker, an instrument), since a mixture is a sum of sources and sparsity forces the atoms to specialize.

This is the first historical bridge between signal processing and unsupervised learning: a learned dictionary is a latent representation, \alpha a code — the same grammar as the autoencoders previewed at the start of this chapter, with the added spectral interpretability of the atoms. The “deep learning” signal reconstruction announced by the sheet is the generalization of this scheme in which the decoder becomes a network trained on degraded/clean pairs (denoising, super-resolution).

Example 8.2 — Dictionary learning and NMF source separation, in numbers Sparse dictionary learning solves \min_{D,\,A}\ \tfrac12\|X - DA\|_F^2 + \lambda\|A\|_1 \quad\text{subject to } \|d_k\|_2 \leq 1, where the columns of X are signal frames, those of D the atoms and those of A the codes.

It alternates two convex steps: sparse coding of A with D fixed (a lasso, solved by iterative soft thresholding), then a least-squares update of D with A fixed.

For audio, non-negative matrix factorization approximates the magnitude spectrogram V \approx WH with W, H \geq 0: the columns of W are spectral templates and the rows of H their activations in time.

Take two notes with templates w_1 = (2, 1, 0) (energy in the low bins) and w_2 = (0, 1, 2) (high bins). A frame v = (4, 3, 2) is explained exactly by the activations h = (2, 1), since 2w_1 + w_2 = (4, 3, 2).

To separate, each time-frequency bin is shared in proportion to the contribution of each source, a Wiener-like soft mask. In the middle bin, source 1 contributes 2 \times 1 = 2 and source 2 contributes 1 \times 1 = 1, so it receives the fraction 2/3. Multiplying the complex STFT of the mixture by each mask and inverting the STFT by overlap-add (chapter 6) yields the two separated signals.

Unrolling a fixed number of soft-thresholding iterations into the layers of a network trained end to end (LISTA, Gregor & LeCun, 2010) is one of the routes from dictionaries to signal reconstruction by deep learning.

9. The teacher’s research: from signals to optimal transport

This self-study chapter shows where the tools of the course lead in current research: from comparing two spectra to comparing whole distributions of signals with optimal transport.

Who teaches this module, and why it is no accident Rémi Flamary (CMAP, École Polytechnique) is co-director of the LLGA master’s program and co-creator of the POT — Python Optimal Transport library (JMLR 2021): one of the most used tools in the field. His research connects precisely the two ends of this course: from signals and images (telecommunications, brain imaging) to distances between distributions of data — that is, from the spectrograms of this module to the OT losses of the M1P2 module.

POT: a library designed for lab practice Architecture (introductory session):

  • Exact solvers: ot.emd (network simplex, exact LP); regularized: ot.sinkhorn and its log-stabilized version; structured: ot.gromov (GW, FGW); barycenters; domain adaptation: ot.da.OTDA — the transported classifier of the TPAMI 2016 paper, packaged.
  • Native backpropagation: the NumPy/PyTorch/JAX backends — an OT loss (spectral loss included) fits into any training loop; the gradient of the regularized cost is the dual potential (proved in the OT module, ch. 3).
  • Impact: ~80 solvers from the literature, a gallery of examples, used from medical to generative audio: the “de facto standard” of numerical OT.

Bridging lab: optimal transport between two segments of a spectrogram

import numpy as np, ot
import scipy.signal as sig

fs = 8000; t = np.arange(2*fs)/fs
rng = np.random.default_rng(0)
x  = np.sin(2*np.pi*440*t) + 0.5*rng.random(2*fs)      # a note with noise
f, T, S = sig.stft(x, fs=fs, nperseg=512, window='hann')

A = np.abs(S[:, T < 0.5]).T                            # 16 frames t < 0.5 s, each a point of R^257
B = np.abs(S[:, T >= 0.5]).T                           # 48 frames t >= 0.5 s
a, b = ot.unif(len(A)), ot.unif(len(B))                # uniform mass per frame
M  = ot.dist(A, B)                                     # 16 x 48 squared distances between spectra
G  = ot.sinkhorn(a, b, M / M.max(), reg=0.05)          # transport plan between the two sets of frames
print("OT cost:", ot.emd2(a, b, M))                    # about 0.013 for this stationary note
# Switch the note to 660 Hz after t = 1 s: the cost rises to about 0.47, because
# two thirds of the frames of B (the 32 after t = 1 s) no longer resemble any frame of A.

Each frame is a magnitude spectrum, a point of \mathbb{R}^F with F = 257 frequency bins, and one “transports” the distribution of the frames of one half toward the other (flattening the array instead would only compare histograms of scalar magnitudes): this is exactly the basic building block of the spectral losses of generative audio ML (ch. 8), and the first lab of the M1P2 module.

A recent bridge: projections, sliced OT and reconstruction Part of the instructor’s recent work deals with sliced transport (random 1D projections): Tanguy, Flamary & Delon, “Reconstructing discrete measures from projections” (C. R. Mathématique, 2024 — arXiv:2304.12029) show how the projections of a discrete measure determine it, with consequences for the empirical sliced Wasserstein distance. The vocabulary is that of the course: projecting a cloud (a measure of spectrograms) onto directions means computing marginal spectra; 1D OT by quantiles (ch. 8) does the rest. Flamary’s other works (FGW, domain adaptation, COOT) are detailed in the Computational Optimal Transport module (M1P2, ch. 7).

Exercises

  1. ★ Square wave series (with solution). Compute the c_k of the odd square wave f = \mathrm{sign}(\sin), check the decay in 1/k and illustrate the Gibbs phenomenon (partial sums at 5, 21, 101 harmonics).

    Solution. f odd and real ⇒ c_k purely imaginary. For k \neq 0 even, the integral of \sin(k\omega_0 t) over the two half-periods cancels. For k odd: c_k = \frac{1}{T_0}\int_0^{T_0}f e^{-ik\omega_0t}dt = \frac{-2i}{T_0}\int_0^{T_0/2}\sin(k\omega_0 t)dt = \frac{-2i}{T_0}\cdot\frac{2\,(1-(-1)^k)}{2k\omega_0} = \frac{-2i}{\pi k}, that is |c_k| = 2/(\pi|k|).

    This is consistent with lemma 2.3: the period carries two jumps (of +2 at 0, of -2 at T_0/2), whose contributions \frac{2}{2\pi i k} and \frac{-2\,(-1)^k}{2\pi i k} add up to \frac{2(1-(-1)^k)}{2\pi i k} = \frac{-2i}{\pi k} \ \text{ for } k \text{ odd}, \qquad 0 \ \text{ for } k \text{ even}. Recombining: f(t) = \frac{4}{\pi}\sum_{m\geq0}\frac{\sin((2m+1)\omega_0t)}{2m+1}.

  2. ★ Numerical Parseval. Check Parseval’s identity on a sampled periodic signal (compare \sum|x_n|^2/N and \sum|X_k|^2/N^2), then explain why the DFT preserves the energy up to a factor N (a Vandermonde matrix, unitary up to a factor).

  3. ★ Aliasing. Generate a chirp f(t) = 100 + 400t Hz sampled at f_s = 500 Hz; compute the spectrogram and interpret the “ghost” trace that goes back down; what apparent frequency does f = 350 Hz correspond to? (prop. 4.2: |350 - 500| = 150 Hz, and the trace rises up to 250 Hz before going back down.)

  4. ★ Homemade FFT. Implement the butterfly recursion of th. 5.3 (even/odd recursive function), validate it against np.fft.fft, and plot the execution time as a function of q = \log_2 N for N = 2^q, q = 4,\dots,18; check the \log·N slope.

  5. ★ Convolution via FFT. Implement the linear convolution via zero-padding + FFT (prop. 5.2), measure the times against np.convolve for n = m = 2^{14}, and explain the crossing of the complexities.

  6. ★ Hann filter. Low-pass FIR design by windowing (example 5.6): ideal response H(\nu) = \mathbf{1}_{|\nu| \leq \nu_c}, multiplication by a Hann window, inverse transform. Plot the frequency response for M = 31 and comment on the transition/ripple trade-off; compare with Blackman at the same M.

  7. ★ Mel bank (with solution). Implement a bank of 26 triangular Mel filters between 0 and f_s/2 (formulas of ch. 8), apply it to a spectrogram, then the DCT; check that a vowel and a fricative have very different MFCCs.

    Solution (key). 26 triangles need 28 edge points equally spaced in Mel: m_i = i\,m_{\max}/27 for i = 0, \dots, 27, with m_{\max} = 2595\log_{10}(1+f_{\max}/700).

    Filter j = 1, \dots, 26 rises from m_{j-1} to m_j and falls back to zero at m_{j+1}. The edges are converted to Hz by the inverse formula, and the triangles are given a normalized area (the highest ones are wider in Hz).

    Energies E_m = \sum_\nu W_m(\nu)|S(\nu)|^2, log, DCT-II: the vowels (formants, low-Mel energy) concentrate on the first MFCCs, the fricatives (high frequencies) on the following ones.

  8. ★★ Quantitative Gibbs (with solution). Show that for the square wave of exercise 1, the maximum of the partial sum S_N converges to \frac{2}{\pi}\int_0^\pi \frac{\sin u}{u}du \approx 1.179, i.e. an overshoot of 0.179 above the level 1, about 9% of the jump (which has height 2).

    Solution (sketch). Near the jump at 0, S_N f(t) \approx \frac{2}{\pi}\,\mathrm{Si}(N\omega_0 t) \quad\text{with}\quad \mathrm{Si}(x) = \int_0^{x}\frac{\sin u}{u}\,du (substitution in the Dirichlet kernel D_N, Definition 0.2) — oscillations concentrated on a width O(1/N).

    The maximum of the sinc integral is reached at the first zero of the numerator, u = \pi, and equals \mathrm{Si}(\pi) \approx 1.852 versus \mathrm{Si}(\infty) = \pi/2. The relative overshoot is \frac{\mathrm{Si}(\pi) - \pi/2}{\pi/2} \approx 17.9\,\% of the half-jump, i.e. \approx 8.95\,\% of the full jump — independent of N; only the width \sim 1/N shrinks.

  9. ★★ Gaussian (with solution). Re-prove prop. 3.3 by completing the square: show that \int e^{-\pi(t + i\nu)^2}dt = \int e^{-\pi t^2}dt = 1 (move the contour, or differentiate \int e^{-\pi(t+i\nu)^2}dt with respect to \nu and conclude), and deduce the two half-widths (standard deviations) \sigma_t, \sigma_\nu and their product.

    Solution (key). -\pi t^2 - 2\pi i\nu t = -\pi(t+i\nu)^2 - \pi\nu^2 \quad\Rightarrow\quad \hat g = e^{-\pi\nu^2}\cdot\int e^{-\pi u^2}du with u = t + i\nu on the horizontal segment — the integral does not depend on \nu (holomorphy + decay at the ends).

    For g = e^{-\pi t^2}, the normalized energy density |g|^2/\|g\|_2^2 \propto e^{-2\pi t^2} is a Gaussian of variance 1/(4\pi): \sigma_t^2 = \sigma_\nu^2 = 1/(4\pi), \qquad\text{hence}\qquad \sigma_t\sigma_\nu = 1/(4\pi): the Gaussian saturates Heisenberg-Gabor (th. 3.4).

  10. ★★ Heisenberg (with solution). Check that equality \sigma_t\sigma_\nu = 1/(4\pi) requires x' = \lambda t x in the proof of th. 3.4, solve this differential equation and conclude that only the Gaussians (translated, modulated, dilated) saturate the inequality.

    Solution (key). Solve x'(t) = \lambda t\,x(t): x(t) = C e^{\lambda t^2/2}, which is in L^2 if and only if \lambda < 0. Writing \lambda = -2\pi/\sigma^2 gives x(t) = Ce^{-\pi t^2/\sigma^2}.

    One then adds a translation in time (center t_0) and a modulation (center \nu_0), which do not change the standard deviations.

  11. ★★ Line detection (with solution). A signal contains x(t) = \sin(2\pi\cdot50t) + A\sin(2\pi\cdot80t) + white noise, with A = 0.01 (a line at −40 dB). Show that the raw periodogram does not detect the line at 80 Hz (th. 7.5), then make it emerge with Welch (Hann, segments of 1024).

    Solution (key). Periodogram: leakage of the rectangular window at −13 dB (prop. 5.4) + 100% variance ⇒ the line at −40 dB is drowned.

    Welch: Hann window (first lobe −31 dB, decay −18 dB/octave) ⇒ the leakage of the 50 Hz line falls below −60 dB; L \approx 2N/1024 segments ⇒ the standard deviation of the estimator divided by \sqrt{L}: the line at 80 Hz exceeds the detection threshold. Code: scipy.signal.welch as in ch. 7.

  12. ★★★ Project — denoising and spectral loss. Denoise a recording (white noise added): spectrogram → time-frequency masking (attenuate the cells below an energy threshold) → reconstruction by inverse STFT (overlap-add, ch. 6); SNR before/after, musical artifacts. Then compare the distributions of original/denoised spectrogram frames by Gaussian W_2 (formula of ch. 8, POT: ot.gaussian.empirical_bures_wasserstein_distance) — direct preparation of the M1P2 module.

Projects and assessment

The official assessment of the course unit is “lab reports (40%) and a final written theoretical and practical exam (60%)”, as stated in the description of the 2026–2027 SynapseS sheet; the LLGA section of the sheet displays no assessment of its own. The subjects below are proposals of these notes, not official rules: a solo or group project can stand for the lab reports, and the exam proposal for the final written exam.

The projects exercise mainly O6 (signals of several modalities: speech, audio, sensors), O1 (the mathematics under the learned models) and O7 (hands-on work on real data), with O4 for the graph-signal subject.

Solo projects

  1. [Solo · ≈ 20 h] Spectral estimation on real ECG recordings. How many averaged segments does a weak line need, and when does a parametric AR spectrum beat Welch?
    • Data: the MIT-BIH Arrhythmia Database (Moody and Mark, 2001; PhysioNet), two-lead ECG at 360 Hz; targets are the heart-rate line, its harmonics, and a weak line added at a known level.
    • Work: periodogram, Bartlett and Welch from the definitions, validated against scipy.signal.welch; empirical bias and variance against the number of segments L, first on simulated AR(2) noise whose spectrum is known, then on ECG excerpts; AR(p) by Yule–Walker, with p chosen on held-out excerpts.
    • Success criterion: the measured variance follows S^2/L within a factor 1.5 (theorem 7.5); the weak line detected at a stated false-alarm rate.
    • Deliverables: code, a 4–6 page report, a 10-minute talk; chapters 5 and 7, exercise 11.
  2. [Solo · ≈ 20 h] Paper study: from the convolution theorem to the GCN (Shuman, Narang, Frossard, Ortega and Vandergheynst, 2013; Kipf and Welling, 2017). Which results of chapters 1–5 survive on a graph?
    • Work: the graph Fourier transform from the Laplacian eigenvectors; exact low-pass and band-pass filters, then the Chebyshev polynomial filters of Defferrard, Bresson and Vandergheynst (2016), which avoid the eigendecomposition; denoising of smooth node signals; the GCN layer of chapter 8 shown to be a degree-1 filter, with its frequency response plotted.
    • Data: the Minnesota road graph of the PyGSP library with synthetic smooth signals, then the Cora citation graph (Planetoid split).
    • Success criterion: a polynomial filter of degree at most 20 within 1 dB of the SNR gain of the exact filter, at a fraction of its cost.
    • Deliverables: code, report, talk; chapters 1, 3, 5 and 8. It leads to Graph ML, ch. 4 and to the Laplacian spectra of Geometric Data, ch. 0.
  3. [Solo · ≈ 20 h] Periodicities and AR baselines for taxi demand. Does an AR model on deseasonalized counts beat the seasonal naive forecast? This is the warm-up of the Y1.2 fleet-dispatch project.
    • Data: NYC TLC Trip Record Data (yellow taxis, three consecutive months), aggregated into hourly pickup counts per Manhattan taxi zone.
    • Work: Welch spectrum of the counts (lines at periods 24 h and 168 h and their harmonics); removal of these lines by a seasonal mean or notch filters (ch. 5); AR(p) by Yule–Walker; one-hour-ahead forecasts.
    • Success criterion: mean absolute error per zone on the held-out month, with a confidence interval over days, against the seasonal naive forecast (same hour of the previous week).
    • Deliverables: code, report, talk; chapters 5 and 7. The zone graph and the forecasts then feed Graph ML and the Y1.2 project.

Group projects

  1. [Group of 3 · ≈ 30 h each] From waveform to speech front end: Mel, MFCC and keyword spotting. Which input representation does a small keyword spotter need, and how fragile is it to the STFT parameters?
    • Data: Speech Commands v2 (Warden, 2018), 35 words in one-second clips at 16 kHz, official split.
    • Work: STFT, Mel bank and MFCC in NumPy (exercise 7), validated against torchaudio or librosa with matching options (HTK Mel formula, same normalization); a small CNN (preview, ch. 8) on the raw waveform, on log-Mel and on MFCC; ablations of the window length and hop (prop. 6.2) and of a naive resampling to 8 kHz without anti-aliasing filter (ch. 4).
    • Split: front end and its tests; models and training; ablations and evaluation.
    • Success criterion: test accuracy per representation over 3 seeds; the aliasing ablation measured and explained.
    • Deliverables: code, an 8–10 page joint report, a 20-minute talk; chapters 4–6 and 8. It continues in Speech Technology and in the voice-first assistant of the Y2.3 project.
  2. [Group of 3 · ≈ 30 h each] Speech denoising: Wiener filter, NMF and a learned mask. How far does a learned mask go beyond the two classical spectral methods, and do the metrics agree?
    • Data: VoiceBank+DEMAND (Valentini-Botinhao, Wang, Takaki and Yamagishi, 2016), paired noisy and clean speech, resampled from 48 to 16 kHz with an anti-aliasing filter.
    • Work: an STFT Wiener filter whose noise PSD is estimated by Welch on speech-free frames (ch. 7); NMF soft masks (example 8.2); a small learned mask network (preview, ch. 8); inversion by overlap-add (ch. 6).
    • Split: one member per method, with shared STFT and evaluation code.
    • Success criterion: SNR and SI-SDR (Le Roux, Wisdom, Erdogan and Hershey, 2019) gains on the test set, plus the Gaussian W_2 between clean and denoised spectrogram frames (ch. 8, POT); one case where the metrics disagree, explained.
    • Deliverables: code, joint report, talk; chapters 6–8 and exercise 12. It leads to Computational Optimal Transport.

Written exam (proposal)

Two hours, one handwritten A4 sheet allowed, a simple calculator allowed; four exercises. It stands for the official final written exam (60%).

  1. Question (5 points). Explain why the raw periodogram is not a consistent estimator of the power spectral density, and how the Bartlett and Welch methods fix it. A record has N = 8192 samples at f_s = 1 kHz; with L = 8 Bartlett segments, give the frequency resolution and the variance reduction.
    • Expected: bias S * F_N with the Fejér kernel, vanishing as N grows; variance about S^2 for every N (theorem 7.5).
    • Averaging L periodograms divides the variance by L, at the price of a resolution f_s/M instead of f_s/N.
    • Here M = 1024: a resolution of about 1 Hz instead of 0.12 Hz, and a variance divided by 8. Welch with Hann windows and 50% overlap uses 15 segments of the same length.
  2. Question (5 points). Consider the filter y_n = 0.5\,y_{n-1} + x_n - x_{n-1}. Give its transfer function, its poles and zeros, decide whether it is stable, compute its gain at \nu = 0 and \nu = 1/2, and give its impulse response.
    • Expected: H(z) = (1 - z^{-1})/(1 - 0.5z^{-1}), with a zero at 1 and a pole at 0.5, inside the unit circle: stable (theorem 5.5).
    • H = 0 at \nu = 0, so constants are removed, and H = 2/1.5 = 4/3 at \nu = 1/2: a high-pass, DC-blocking filter.
    • h_0 = 1 and h_n = -0.5^n for n \geq 1, so \sum_n |h_n| = 2 < \infty.
  3. Question (5 points). Prove the circular convolution theorem (proposition 5.2). Then explain how to compute the linear convolution of two signals of lengths 300 and 200 with radix-2 FFTs, and give the FFT size.
    • Expected: for fixed m, p = (n - m) \bmod N is a bijection of \mathbb{Z}/N\mathbb{Z}, and \omega^{kN} = 1 splits the double sum into X_k H_k.
    • Zero-pad both signals to N' \geq 300 + 200 - 1 = 499, so N' = 512: the wrap-around then falls on zeros.
  4. Question (5 points). A speech front end uses 25 ms Hann windows every 10 ms at 16 kHz, with a 512-point FFT. Give the window and the hop in samples, the bin spacing and the half-width of the main lobe, and decide whether plain overlap-add inverts this STFT. Then compute the Mel value of 1 kHz and say why MFCCs apply a DCT.
    • Expected: L = 400 and a hop of 160 samples; bins every 16000/512 = 31.25 Hz; a Hann main lobe of \pm 2/T = \pm 80 Hz (proposition 6.2).
    • The hop 0.4L is not COLA for Hann (the overlapped windows oscillate between about 1.19 and 1.31), but \sum_m w^2 never vanishes, so weighted overlap-add inverts it exactly (ch. 6).
    • 2595\log_{10}(1 + 1000/700) \approx 1000 Mel. The DCT decorrelates the log-Mel energies of neighbouring bands and compacts them into a few coefficients (definition 6.5).

Grading rubric (proposal)

Criterion Weight What earns the marks
Signal model and question 15 % the signal, the noise and the parameters (rate, window, hop, number of segments) stated and justified from the chapters
Correct implementation 25 % NumPy code written from the definitions and validated against a reference library with a stated tolerance
Experimental protocol 20 % held-out data, several excerpts or seeds, uncertainty reported, no leakage between splits
Analysis tied to the theory 25 % each result, failures included, explained by a statement of the course (resolution, aliasing, leakage, variance)
Communication 15 % figures with units that answer the question, a talk on time, honest limits

References