Fourier and harmonic analysis

In 1807 Joseph Fourier told the Paris Academy that the temperature in a heated bar could be computed by writing any initial temperature profile as a sum of sines and cosines, each of which evolves independently: $e^{inx}$ simply decays like $e^{-n^2 t}$ under the heat equation. Lagrange objected that a discontinuous function could not possibly be a sum of smooth waves. Both were partly right, and settling the dispute took a century: Dirichlet's convergence theorem (1829), Riemann's integral (1854), Cantor's set theory (born from the study of sets where trigonometric series may misbehave), and Lebesgue's integral and the Riesz–Fischer theorem (1902–1907) all came out of it.

Today Fourier analysis is everywhere: it diagonalises every translation-invariant linear operation (differential operators, filters, convolutions), it is how MRI scanners and phone networks work, it proves theorems about primes, and its uncertainty principle is a law of physics. Harmonic analysis is the general theory: decompose functions on a group into characters, and study how size, smoothness and localisation on one side are reflected on the other.

Why exponentials? Translation $\tau_a f(x) = f(x - a)$ acts on functions, and $e^{inx}$ is an eigenvector of every translation at once: $\tau_a e^{inx} = e^{-ina} e^{inx}$. So any operator that commutes with translations — differentiation, the heat flow, convolution with a fixed kernel — is diagonal in the basis of exponentials. Differentiation becomes multiplication by $in$; convolution becomes multiplication of coefficients.

Three dictionaries follow, and recur throughout the subject:

• Smoothness ↔ decay. If $f$ has $k$ integrable derivatives, $|\hat f(n)| = O(|n|^{-k})$; jumps produce $1/n$ decay and the slow, oscillating convergence of the Gibbs phenomenon.
• Convolution ↔ multiplication. $\widehat{f * g} = \hat f \, \hat g$. Averaging a function against a peaked kernel smooths it, and on the Fourier side multiplies by a factor near $1$ at low frequencies.
• Localisation ↔ spread. A function concentrated near a point has a Fourier transform spread over many frequencies, and conversely. The Heisenberg inequality quantifies this, and the Gaussian is the perfect compromise.

The deep difficulties are about limits: in what sense does $\sum \hat f(n) e^{inx}$ converge to $f$? In $L^2$, always (a Hilbert space fact). Uniformly, after averaging (Fejér). Pointwise, the answer is subtle and was only settled in 1966.

Write $\mathbb T = \R / 2\pi\Z$ for the circle, with functions on $\mathbb T$ identified with $2\pi$-periodic functions on $\R$.

With this normalisation of the transform no factors of $2\pi$ appear in the inversion formula $f(x) = \int \hat f(\xi) e^{2\pi i x \xi} \, d\xi$ or in Plancherel's theorem. Other books use $e^{-ix\xi}$ and put $\frac{1}{2\pi}$ or $\frac{1}{\sqrt{2\pi}}$ elsewhere; always check the convention.

Applied to the square wave, Parseval gives $\sum_{n \text{ odd}} \frac{4}{\pi^2 n^2} = 1$, that is, $1 + \frac19 + \frac1{25} + \cdots = \frac{\pi^2}{8}$, and hence $\sum_{n \ge 1} \frac1{n^2} = \frac{\pi^2}{6}$.

Applied to $f(x) = e^{-\pi t x^2}$, Poisson summation gives Jacobi's identity $\theta(t) = t^{-1/2} \theta(1/t)$ for $\theta(t) = \sum_n e^{-\pi n^2 t}$. Riemann (1859) wrote $\pi^{-s/2}\Gamma(s/2)\zeta(s) = \int_0^\infty \frac{\theta(t) - 1}{2} t^{s/2} \frac{dt}{t}$ and used Jacobi's identity to prove the functional equation $\xi(s) = \xi(1 - s)$ of the zeta function: see analytic number theory.

In quantum mechanics a particle's position density is $|f|^2$ and its momentum density is $|\hat f|^2$ (in units with $p = h\xi$); the theorem says the product of the standard deviations of position and momentum is at least $h/(4\pi) = \hbar/2$. This is a theorem of Fourier analysis, applied to a physical model.

Pointwise convergence of $S_N f(x)$ is delicate, because the Dirichlet kernel is not an approximate identity. What is known:

• Positive results for nice functions. If $f$ is piecewise $C^1$, then $S_N f(x) \to \frac12\big(f(x^+) + f(x^-)\big)$ everywhere (Dirichlet, 1829), uniformly on closed intervals avoiding jumps. If $f$ is Hölder continuous, or of bounded variation and continuous, $S_N f \to f$ uniformly (Dini–Lipschitz, Jordan).
• Continuity is not enough. du Bois-Reymond (1873) constructed a continuous function whose Fourier series diverges at a point. A quick existence proof: the functionals $f \mapsto S_N f(0)$ on the Banach space $C(\mathbb T)$ have norms $\frac{1}{2\pi}\|D_N\|_1 \to \infty$, so by the uniform boundedness principle they cannot be pointwise bounded. Kahane and Katznelson (1966) showed that for every null set $E \subset \mathbb T$ there is a continuous $f$ whose Fourier series diverges at every point of $E$.
• $L^1$ is hopeless. Kolmogorov (1923) found $f \in L^1(\mathbb T)$ whose Fourier series diverges almost everywhere, and (1926) one that diverges everywhere.
• $L^2$ is fine. Lusin conjectured (1915) and Carleson proved (1966) that the Fourier series of every $f \in L^2(\mathbb T)$ converges to $f$ almost everywhere; Hunt (1968) extended this to $L^p$, $p > 1$. The proof is among the hardest in classical analysis; it controls the maximal operator $\sup_N |S_N f|$ by decomposing time–frequency space into tiles.
• Averaging rescues everything. For $f \in L^1$, the Fejér means converge to $f$ in $L^1$ and almost everywhere (at every Lebesgue point), and uniformly for continuous $f$.

So the right summary is not "Fourier series converge" but: in $L^2$, always; after averaging, always; pointwise, almost everywhere exactly in the range $p > 1$, with the Gibbs overshoot near jumps never going away.

The three settings above are one theory. For a locally compact abelian group $G$, a character is a continuous homomorphism $\chi : G \to \mathbb S^1$; the characters form the dual group $\widehat G$, and the Fourier transform is $\hat f(\chi) = \int_G f(x) \overline{\chi(x)} \, dx$ with respect to Haar measure. Pontryagin duality (1934) says $\widehat{\widehat G} \cong G$, and Plancherel holds: with suitably normalised Haar measures, $f \mapsto \hat f$ is unitary from $L^2(G)$ to $L^2(\widehat G)$.

Orthogonality of characters on $\Z/N$ is the geometric-series identity $\sum_{k < N} z^k = 0$ for an $N$-th root of unity $z \ne 1$ (formalised below). For compact non-abelian groups, characters are replaced by irreducible representations: the Peter–Weyl theorem (1927) decomposes $L^2(G)$ into matrix coefficients, and for the rotation group acting on the sphere it yields spherical harmonics. Harmonic analysis on non-compact Lie groups such as $SL_2(\R)$ was developed by Gelfand, Bargmann and above all Harish-Chandra, and leads to the theory of automorphic forms.

The finite model first: orthogonality of the characters of $\Z/N$, proved from the geometric series.

/-- **Orthogonality of characters on `ℤ/N`.** If `z^N = 1` and `z ≠ 1`, then
`1 + z + ⋯ + z^{N-1} = 0`. Proof: `(∑ zᵏ)(z − 1) = z^N − 1 = 0` and `z − 1 ≠ 0`. With
`z = e^{2πi(m−n)/N}` this says the discrete exponentials are orthogonal — the finite model of
Fourier series and the basis of the fast Fourier transform. -/
theorem roots_of_unity_sum {z : ℂ} {N : ℕ} (hz : z ^ N = 1) (hz1 : z ≠ 1) :
    ∑ k ∈ Finset.range N, z ^ k = 0 := by
  have h : (∑ k ∈ Finset.range N, z ^ k) * (z - 1) = 0 := by
    rw [geom_sum_mul, hz, sub_self]
  exact (mul_eq_zero.mp h).resolve_right (sub_ne_zero.mpr hz1)

Mathlib works on AddCircle T $= \R / T\Z$ with Haar measure normalised to total mass $1$, so its coefficients are $\hat f(n) = \frac1T \int_0^T f(x) e^{-2\pi i n x/T} dx$. The exponentials are orthonormal (and, by Stone–Weierstrass rather than Fejér, a Hilbert basis), and Parseval's identity holds:

/-- **The exponentials are orthonormal** in `L²` of the circle `ℝ/Tℤ` with normalised Haar
measure: `⟪eₘ, eₙ⟫ = δₘₙ` for `eₙ(x) = e^{2πinx/T}`. This is Mathlib's `orthonormal_fourier`;
`span_fourierLp_closure_eq_top` adds that their span is dense, so they form a Hilbert basis. -/
theorem orthonormal_exponentials (T : ℝ) [Fact (0 < T)] :
    Orthonormal ℂ (fun n : ℤ => fourierLp (T := T) 2 n) :=
  orthonormal_fourier

/-- **Parseval's identity.** For `f ∈ L²(ℝ/Tℤ)`, `∑ₙ |f̂(n)|² = ∫ |f|²` (normalised measure).
This is Mathlib's `hasSum_sq_fourierCoeff`, a consequence of the exponentials forming a Hilbert
basis. -/
theorem parseval (T : ℝ) [Fact (0 < T)] (f : Lp ℂ 2 (@AddCircle.haarAddCircle T _)) :
    HasSum (fun n : ℤ => ‖fourierCoeff f n‖ ^ 2)
      (∫ t : AddCircle T, ‖f t‖ ^ 2 ∂AddCircle.haarAddCircle) :=
  hasSum_sq_fourierCoeff f

From Parseval we derive the Riemann–Lebesgue lemma for $L^2$ functions, exactly as in the proof above:

/-- **Riemann–Lebesgue lemma for `L²` functions on the circle**, from Parseval: the series
`∑ |f̂(n)|²` converges, so its terms tend to `0` (along the cofinite filter, i.e. as `|n| → ∞`),
and hence so do the coefficients `f̂(n) = √(|f̂(n)|²)`. -/
theorem fourier_coeff_tendsto_zero (T : ℝ) [Fact (0 < T)]
    (f : Lp ℂ 2 (@AddCircle.haarAddCircle T _)) :
    Tendsto (fun n : ℤ => fourierCoeff f n) cofinite (𝓝 0) := by
  have hsq : Tendsto (fun n : ℤ => ‖fourierCoeff f n‖ ^ 2) cofinite (𝓝 0) :=
    (parseval T f).summable.tendsto_cofinite_zero
  have hnorm : Tendsto (fun n : ℤ => ‖fourierCoeff f n‖) cofinite (𝓝 0) := by
    have h := (Real.continuous_sqrt.tendsto 0).comp hsq
    rw [Real.sqrt_zero] at h
    refine h.congr fun n => ?_
    simp only [Function.comp_apply, Real.sqrt_sq (norm_nonneg _)]
  exact tendsto_zero_iff_norm_tendsto_zero.mpr hnorm

On $\R$, Mathlib writes 𝐞 t for $e^{2\pi i t}$ and 𝓕 for the transform $\hat f(\xi) = \int f(x) e^{-2\pi i x\xi} dx$, the convention of this lesson. The trivial bound $|\hat f(\xi)| \le \|f\|_1$, proved directly; the Riemann–Lebesgue lemma; Plancherel on $L^2(\R)$; and the self-dual Gaussian:

/-- **The Fourier transform is bounded from `L¹` to `L^∞`:** `|f̂(ξ)| ≤ ∫ |f|` for every `ξ`,
because `|∫ g| ≤ ∫ |g|` and the character `e^{−2πixξ}` has modulus `1`. (For non-integrable `f`
both sides are `0` under Lean's convention for the Bochner integral.) -/
theorem fourier_bound (f : ℝ → ℂ) (ξ : ℝ) :
    ‖∫ x : ℝ, 𝐞 (-(ξ * x)) • f x‖ ≤ ∫ x : ℝ, ‖f x‖ := by
  calc ‖∫ x : ℝ, 𝐞 (-(ξ * x)) • f x‖ ≤ ∫ x : ℝ, ‖𝐞 (-(ξ * x)) • f x‖ :=
        norm_integral_le_integral_norm _
    _ = ∫ x : ℝ, ‖f x‖ := by
        congr 1
        ext x
        simp

/-- **Riemann–Lebesgue lemma** on `ℝ`: for integrable `f`, `f̂(ξ) = ∫ e^{−2πixξ} f(x) dx → 0`
as `|ξ| → ∞`. This applies Mathlib's `tendsto_integral_exp_inner_smul_cocompact` (stated on any
finite-dimensional real inner product space; for `ℝ` the cocompact filter is `|ξ| → ∞`). -/
theorem riemann_lebesgue (f : ℝ → ℂ) :
    Tendsto (fun ξ : ℝ => ∫ x : ℝ, 𝐞 (-(ξ * x)) • f x) (cocompact ℝ) (𝓝 0) := by
  simpa using tendsto_integral_exp_inner_smul_cocompact f

/-- **Plancherel's theorem.** The Fourier transform, defined on `L¹ ∩ L²` and extended by
continuity, is an isometry of `L²(ℝ)`: `‖𝓕 f‖₂ = ‖f‖₂` (indeed a unitary map, with inverse the
conjugate transform). This is Mathlib's `MeasureTheory.Lp.norm_fourier_eq`. -/
theorem plancherel (f : Lp (α := ℝ) ℂ 2) : ‖𝓕 f‖ = ‖f‖ :=
  Lp.norm_fourier_eq f

/-- **The Gaussian `e^{−πx²}` is its own Fourier transform**, with the normalisation
`f̂(ξ) = ∫ f(x) e^{−2πixξ} dx`. Obtained from Mathlib's `fourier_gaussian_pi`, which computes the
transform of `e^{−πbx²}` for `Re b > 0`, by setting `b = 1` and simplifying `1/1^{1/2} = 1`. -/
theorem gaussian_self_dual :
    (𝓕 fun x : ℝ => cexp (-π * x ^ 2)) = fun ξ : ℝ => cexp (-π * ξ ^ 2) := by
  have h := fourier_gaussian_pi (b := 1) (by simp)
  simpa using h

Poisson summation for Schwartz functions, the theta functional equation that Riemann used, and $\zeta(2) = \pi^2/6$, which Mathlib proves with Fourier series of Bernoulli polynomials:

/-- **Poisson summation formula.** For a Schwartz function `f` on `ℝ`,
`∑ₙ f(x + n) = ∑ₖ f̂(k) e^{2πikx}`; at `x = 0`, `∑ f(n) = ∑ f̂(k)`. The left side is a
`1`-periodic function whose Fourier coefficients are the `f̂(k)`. This is Mathlib's
`SchwartzMap.tsum_eq_tsum_fourier`. -/
theorem poisson_summation (f : SchwartzMap ℝ ℂ) (x : ℝ) :
    ∑' n : ℤ, f (x + n) = ∑' k : ℤ, 𝓕 f k * fourier k (x : UnitAddCircle) :=
  SchwartzMap.tsum_eq_tsum_fourier f x

/-- **Jacobi's theta functional equation** `θ(−1/τ) = (−iτ)^{1/2} θ(τ)` for
`θ(τ) = ∑ₙ e^{πin²τ}` on the upper half-plane — Poisson summation applied to a Gaussian. Riemann
derived the functional equation of `ζ(s)` from it. This is Mathlib's `jacobiTheta_S_smul`. -/
theorem theta_functional_equation (τ : UpperHalfPlane) :
    jacobiTheta ↑(ModularGroup.S • τ) = (-I * τ) ^ (1 / 2 : ℂ) * jacobiTheta τ :=
  jacobiTheta_S_smul τ

/-- **`∑ 1/n² = π²/6`**, which Mathlib proves with Fourier series (Parseval applied to Bernoulli
polynomials). The `n = 0` term is `1/0 = 0` in Lean's conventions. This is Mathlib's
`hasSum_zeta_two`. -/
theorem basel : HasSum (fun n : ℕ => (1 : ℝ) / (n : ℝ) ^ 2) (π ^ 2 / 6) :=
  hasSum_zeta_two

Not formalised (here or, as far as we know, in Mathlib): Fejér's theorem and approximate identities on the circle in this form (Mathlib has related results for convolutions with bump functions), the Gibbs phenomenon, the Heisenberg inequality, Carleson's theorem, and Pontryagin duality for general locally compact abelian groups (Mathlib has the Pontryagin dual group and duality for finite abelian groups).

The first experiment samples the partial sums of the square wave near the jump: the maximum stays at about $1.179$ for $N = 9, 49, 199, 999$, matching $\frac2\pi \operatorname{Si}(\pi)$, while the Fejér means never exceed $1$. It also computes the exact $L^2$ error of $S_N$ from Parseval and watches $\sum_{k \text{ odd}} k^{-2}$ approach $\pi^2/8$.

"""Gibbs overshoot, Fejer's cure and Parseval, for the square wave.

The square wave f = sign(sin x) has Fourier series (4/pi) sum over odd k of
sin(k x)/k. Its partial sums S_N overshoot the jump at 0 by about 9% of the
jump however large N is (the Gibbs phenomenon, limit (2/pi) Si(pi)). The
Fejer means (averages of S_0, ..., S_N) never exceed 1, because the Fejer
kernel is a positive approximate identity. Parseval's identity gives the
exact L^2 error of S_N and the sum 1 + 1/9 + 1/25 + ... = pi^2/8.
"""

import math

import numpy as np

def partial_sum(n, x):
    """S_n(x) for the square wave: (4/pi) * sum of sin(kx)/k, odd k <= n."""
    x = np.asarray(x, dtype=float)
    total = np.zeros_like(x)
    for k in range(1, n + 1, 2):
        total += np.sin(k * x) / k
    return 4 / math.pi * total

def fejer_mean(n, x):
    """sigma_n(x) = (S_0 + S_1 + ... + S_n)(x) / (n + 1)."""
    x = np.asarray(x, dtype=float)
    current = np.zeros_like(x)
    total = np.zeros_like(x)
    for m in range(n + 1):
        if m % 2 == 1:
            current = current + 4 / math.pi * np.sin(m * x) / m
        total += current
    return total / (n + 1)

def sine_integral(t, steps=2000):
    """Si(t) = integral of sin(u)/u from 0 to t, by Simpson's rule."""
    if steps % 2:
        steps += 1
    h = t / steps
    def g(u):
        return 1.0 if u == 0 else math.sin(u) / u
    s = g(0) + g(t)
    s += sum((4 if i % 2 else 2) * g(i * h) for i in range(1, steps))
    return s * h / 3

def l2_error_squared(n):
    """(1/pi) * integral of (f - S_n)^2 over (-pi, pi), by Parseval:
    2 minus the sum of (4/(pi k))^2 over odd k <= n."""
    return 2 - sum((4 / (math.pi * k)) ** 2 for k in range(1, n + 1, 2))

def odd_reciprocal_squares(n):
    """1 + 1/9 + 1/25 + ... over odd k <= n."""
    return sum(1 / k ** 2 for k in range(1, n + 1, 2))

def main():
    grid = np.linspace(1e-6, math.pi / 2, 20001)
    gibbs = 2 / math.pi * sine_integral(math.pi)
    print("Square wave sign(sin x): jump of size 2 at x = 0")
    print(f"Gibbs limit (2/pi) Si(pi) = {gibbs:.6f}")
    print(f"{'N':>5} {'max S_N':>10} {'overshoot/jump':>15} {'max Fejer':>10}")
    for n in [9, 49, 199, 999]:
        top = float(np.max(partial_sum(n, grid)))
        fej = float(np.max(fejer_mean(n, grid)))
        print(f"{n:>5} {top:>10.6f} {(top - 1) / 2:>15.4f} {fej:>10.6f}")
    print("The overshoot does not shrink; it moves towards the jump.")
    print("The Fejer means stay below sup |f| = 1 (positive kernel).")
    print()
    print("Parseval: L^2 error of S_N, and sum over odd k of 1/k^2")
    for n in [1, 9, 99, 999]:
        print(f"  N = {n:>3}: error^2 = {l2_error_squared(n):.6f},"
              f" partial sum = {odd_reciprocal_squares(n):.6f}")
    print(f"  pi^2/8 = {math.pi ** 2 / 8:.6f}")
    print()
    print("Epistemic status: numerical evidence. Maxima are sampled on a grid of")
    print("20001 points and Si(pi) is computed by Simpson's rule; finitely many N")
    print("are shown. That the overshoot tends to (2/pi) Si(pi) - 1 and that")
    print("Fejer means of continuous functions converge uniformly are theorems")
    print("proved in the lesson, not consequences of this table.")

if __name__ == "__main__":
    main()

The second checks Jacobi's identity $\theta(t) = t^{-1/2}\theta(1/t)$ to machine precision, and computes the uncertainty ratio $\big(\int x^2|f|^2\big)\big(\int \xi^2 |\hat f|^2\big) \big/ \big(\|f\|_2^4/16\pi^2\big)$ for three functions with known transforms: $1$ for the Gaussian, $2$ for $e^{-|x|}$ and $6/5$ for a tent.

"""Poisson summation for the Gaussian, and the uncertainty principle.

Poisson summation applied to exp(-pi t x^2) gives Jacobi's identity
theta(t) = t^(-1/2) theta(1/t) for theta(t) = sum over n of exp(-pi n^2 t).
We check it numerically for several t. Then we compute the Heisenberg
uncertainty ratio (int x^2 |f|^2)(int xi^2 |f^|^2) / (||f||^4 / 16 pi^2) for
three functions with known Fourier transforms: it equals 1 for the Gaussian
and exceeds 1 otherwise, as the uncertainty inequality predicts.
"""

import math

import numpy as np

def theta(t, terms=60):
    """theta(t) = sum over |n| <= terms of exp(-pi n^2 t), for t > 0."""
    if t <= 0:
        raise ValueError("t must be positive")
    n = np.arange(-terms, terms + 1)
    return float(np.sum(np.exp(-math.pi * n.astype(float) ** 2 * t)))

def trapezoid(values, grid):
    """Trapezoidal rule on a uniform grid."""
    h = grid[1] - grid[0]
    return float(h * (np.sum(values) - 0.5 * (values[0] + values[-1])))

#: Functions f with their Fourier transforms f^(xi) = int f(x) e^(-2 pi i x xi) dx.
EXAMPLES = {
    "Gaussian exp(-pi x^2)": (
        lambda x: np.exp(-math.pi * x ** 2),
        lambda s: np.exp(-math.pi * s ** 2),
    ),
    "exp(-|x|)": (
        lambda x: np.exp(-np.abs(x)),
        lambda s: 2 / (1 + 4 * math.pi ** 2 * s ** 2),
    ),
    "tent max(1 - |x|, 0)": (
        lambda x: np.maximum(1 - np.abs(x), 0),
        lambda s: np.sinc(s) ** 2,
    ),
}

def uncertainty_ratio(f, fhat, half_width=400.0, points=400001):
    """(int x^2 f^2)(int xi^2 fhat^2) divided by ||f||_2^4 / (16 pi^2)."""
    grid = np.linspace(-half_width, half_width, points)
    fx, fs = f(grid), fhat(grid)
    spread_x = trapezoid(grid ** 2 * fx ** 2, grid)
    spread_s = trapezoid(grid ** 2 * fs ** 2, grid)
    norm_sq = trapezoid(fx ** 2, grid)
    return spread_x * spread_s / (norm_sq ** 2 / (16 * math.pi ** 2))

def main():
    print("Jacobi: theta(t) = t^(-1/2) theta(1/t), theta(t) = sum exp(-pi n^2 t)")
    for t in [0.1, 0.5, 1.0, 2.0, 7.0]:
        lhs, rhs = theta(t), theta(1 / t) / math.sqrt(t)
        print(f"  t = {t:>4}: {lhs:.15f}  {rhs:.15f}  diff {abs(lhs - rhs):.1e}")
    print("  (at t = 0.1 the left side needs many terms; the right side, few)")
    print()
    print("Uncertainty: (int x^2|f|^2)(int xi^2|f^|^2) >= ||f||^4 / (16 pi^2)")
    for name, (f, fhat) in EXAMPLES.items():
        print(f"  {name:<22} ratio = {uncertainty_ratio(f, fhat):.4f}")
    print("  Exact values: 1, 2 and 6/5. Equality holds only for Gaussians.")
    print()
    print("Epistemic status: numerical. Jacobi's identity is a theorem (Poisson")
    print("summation); here it is confirmed to about 1e-15 for five values of t.")
    print("The ratios are trapezoidal approximations on [-400, 400] (slowly")
    print("decaying transforms lose a little mass, about 1e-3); the inequality")
    print("itself, ratio >= 1, is proved in the lesson.")

if __name__ == "__main__":
    main()

These are numerical confirmations of theorems proved above, on finitely many cases; they cannot replace the proofs, but they make the constants concrete.

Build the square wave from sines and watch the overshoot near the jump. Increase the number of terms: the bump gets narrower, but its height does not shrink.

The kernels themselves explain the difference between partial sums and Fejér means. Move $N$: the Dirichlet kernel $D_N$ (whose convolution gives $S_N f$) has large negative lobes, while the Fejér kernel $F_N$ is positive and concentrates at $0$ — an approximate identity.