Exercises Notebook

Converted from exercises.ipynb for web reading.

§20-05 Wavelets — Exercises

8 exercises from Haar DWT by hand to scattering networks and wavelet denoising.

Format	Description
Problem	Markdown cell with task
Your Solution	Code scaffold
Solution	Reference implementation with PASS/FAIL checks

Difficulty Levels

Level	Exercises	Focus
★	1-3	Haar DWT, admissibility, MRA axioms
★★	4-6	Mallat algorithm, Daubechies filters, 2D DWT
★★★	7-8	Scattering networks, wavelet denoising

Topic Map

Topic	Exercise
Haar DWT, Parseval	1
Admissibility, zero mean	2
MRA axioms (Shannon)	3
Mallat algorithm from scratch	4
Daubechies filter verification	5
2D DWT, image compression	6
Scattering transform	7
Wavelet denoising + sparse coding	8

Code cell 2

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

try:
    import seaborn as sns
    sns.set_theme(style="whitegrid", palette="colorblind")
    HAS_SNS = True
except ImportError:
    plt.style.use("seaborn-v0_8-whitegrid")
    HAS_SNS = False

mpl.rcParams.update({
    "figure.figsize":    (10, 6),
    "figure.dpi":         120,
    "font.size":           13,
    "axes.titlesize":      15,
    "axes.labelsize":      13,
    "xtick.labelsize":     11,
    "ytick.labelsize":     11,
    "legend.fontsize":     11,
    "legend.framealpha":   0.85,
    "lines.linewidth":      2.0,
    "axes.spines.top":     False,
    "axes.spines.right":   False,
    "savefig.bbox":       "tight",
    "savefig.dpi":         150,
})
np.random.seed(42)
print("Plot setup complete.")

Code cell 3

import numpy as np
import numpy.linalg as la

try:
    import matplotlib.pyplot as plt
    plt.style.use('seaborn-v0_8-whitegrid')
    plt.rcParams['figure.figsize'] = [10, 5]
    plt.rcParams['font.size'] = 12
    HAS_MPL = True
except ImportError:
    HAS_MPL = False

try:
    import pywt
    HAS_PYWT = True
    print(f'PyWavelets {pywt.__version__} available')
except ImportError:
    HAS_PYWT = False
    print('PyWavelets not found: pip install PyWavelets')

np.set_printoptions(precision=6, suppress=True)
np.random.seed(42)

COLORS = {
    'primary':'#0077BB','secondary':'#EE7733','tertiary':'#009988',
    'error':'#CC3311','neutral':'#555555','highlight':'#EE3377'
}

def header(title):
    print('\n' + '='*len(title))
    print(title)
    print('='*len(title))

def check_close(name, got, expected, tol=1e-8):
    ok = np.allclose(got, expected, atol=tol, rtol=tol)
    print(f"{'PASS' if ok else 'FAIL'} — {name}")
    if not ok:
        print(f'  expected: {expected}')
        print(f'  got:      {got}')
    return ok

def check_true(name, cond):
    print(f"{'PASS' if cond else 'FAIL'} — {name}")
    return cond

print('Setup complete.')

Exercise 1 ★ — Haar DWT by Hand

Given the signal $f = [4, 6, 10, 12, 8, 6, 5, 5]$ :

(a) Compute one level of Haar DWT:

a^{(1)}_k = (f_{2k} + f_{2k+1})/\sqrt{2}, \qquad d^{(1)}_k = (f_{2k} - f_{2k+1})/\sqrt{2}

(b) Apply one more level to $a^{(1)}$ to get $a^{(2)}, d^{(2)}$ .

(c) Verify Parseval: $\|f\|^2 = \|a^{(2)}\|^2 + \|d^{(2)}\|^2 + \|d^{(1)}\|^2$ .

(d) Reconstruct $f$ from $\{a^{(2)}, d^{(2)}, d^{(1)}\}$ and verify exact reconstruction.

Code cell 5

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 6

# Solution
# Exercise 1: Solution

f = np.array([4, 6, 10, 12, 8, 6, 5, 5], dtype=float)

def haar_step(x):
    """One step Haar DWT: (a, d)."""
    n = len(x) // 2
    a = (x[0::2] + x[1::2]) / np.sqrt(2)
    d = (x[0::2] - x[1::2]) / np.sqrt(2)
    return a, d

def haar_istep(a, d):
    """One step Haar IDWT."""
    x = np.zeros(2*len(a))
    x[0::2] = (a + d) / np.sqrt(2)
    x[1::2] = (a - d) / np.sqrt(2)
    return x

# (a)
a1, d1 = haar_step(f)
# (b)
a2, d2 = haar_step(a1)
# (d)
a1_rec = haar_istep(a2, d2)
f_rec  = haar_istep(a1_rec, d1)

header('Exercise 1: Haar DWT')
print(f'(a) a1 = {np.round(a1, 4)}')
print(f'    d1 = {np.round(d1, 4)}')
print(f'(b) a2 = {np.round(a2, 4)}')
print(f'    d2 = {np.round(d2, 4)}')

energy_f = np.sum(f**2)
energy_w = np.sum(a2**2) + np.sum(d2**2) + np.sum(d1**2)
print(f'(c) ||f||² = {energy_f:.6f}, sum coeffs² = {energy_w:.6f}')
check_close('Parseval', energy_w, energy_f)
check_close('Perfect reconstruction', f_rec, f)
print('\nTakeaway: Haar DWT is lossless — every bit of information is preserved.')

Exercise 2 ★ — Admissibility and Zero Mean

(a) Verify $\int_{-\infty}^\infty \psi_{Haar}(t)\,dt = 0$ analytically (by computing the integral).

(b) Compute $\hat{\psi}_{Haar}(\xi) = \int \psi(t)e^{-2\pi i\xi t}\,dt$ in closed form. Show $\hat{\psi}(0) = 0$ .

(c) Show the Gaussian $\psi(t) = e^{-t^2/2}$ fails admissibility: compute $\hat{\psi}(0)$ numerically.

(d) Compute the admissibility constant for the Mexican Hat wavelet numerically:

C_\psi = \int_0^\infty \frac{|\hat{\psi}(\xi)|^2}{\xi}\, d\xi, \qquad\hat{\psi}(\xi) = \sqrt{2\pi}\,\xi^2 e^{-\xi^2/2}

Code cell 8

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 9

# Solution
# Exercise 2: Solution

from scipy.integrate import quad

header('Exercise 2: Admissibility')

# (a) Haar: integral is 0 (piecewise: +1 on [0,0.5), -1 on [0.5,1))
int_haar = 1*0.5 + (-1)*0.5
check_close('(a) Haar zero mean', int_haar, 0.0)

# (b) Haar Fourier transform
def haar_fourier(xi):
    """FT of Haar wavelet: integral of psi*exp(-2pi*i*xi*t)."""
    if xi == 0:
        return 0.0 + 0j
    # int_0^{1/2} e^{-2pi*i*xi*t} dt - int_{1/2}^1 e^{-2pi*i*xi*t} dt
    w = 2*np.pi*xi
    I1 = (1 - np.exp(-1j*w/2)) / (1j*w)
    I2 = (np.exp(-1j*w/2) - np.exp(-1j*w)) / (1j*w)
    return I1 - I2

check_close('(b) Haar hat_psi(0) = 0', abs(haar_fourier(0.0)), 0.0)
xi_vals = np.array([0.5, 1.0, 2.0])
hat_vals = np.array([haar_fourier(xi) for xi in xi_vals])
print(f'    hat_psi at xi={xi_vals}: {np.abs(hat_vals).round(4)}')

# (c) Gaussian has non-zero mean
gaussian_hat_0 = np.sqrt(2*np.pi)  # FT of Gaussian at 0 = integral = sqrt(2pi)
print(f'(c) Gaussian: hat_psi(0) = {gaussian_hat_0:.4f} ≠ 0 → NOT admissible')
check_true('Gaussian fails admissibility', gaussian_hat_0 > 0.1)

# (d) Mexican Hat admissibility constant
def mexhat_psd(xi):
    if xi <= 0:
        return 0.0
    hat_psi = np.sqrt(2*np.pi) * xi**2 * np.exp(-xi**2/2)
    return hat_psi**2 / xi

C_psi, _ = quad(mexhat_psd, 0, np.inf)
print(f'(d) Mexican Hat admissibility constant C_psi = {C_psi:.6f}')
check_true('C_psi finite (admissible)', 0 < C_psi < 100)
print('\nTakeaway: Zero mean ⟺ hat_psi(0)=0 is the minimal condition for admissibility.')

Exercise 3 ★ — MRA Axiom Verification (Shannon MRA)

The Shannon MRA defines:

V_j = \{f \in L^2(\mathbb{R}) : \hat{f}(\xi) = 0 \text{ for } |\xi| > 2^{-j}\pi\}

with scaling function $\phi(t) = \text{sinc}(\pi t) = \sin(\pi t)/(\pi t)$ .

(a) Prove $V_1 \subset V_0$ from the definition.

(b) Verify numerically that $\{\phi(t-k)\}_{k\in\mathbb{Z}}$ is orthonormal: $\langle\phi(\cdot),\phi(\cdot-k)\rangle = \delta_{k0}$ for $k = -3,\ldots,3$ .

(c) The Shannon two-scale relation gives $H(\xi) = \mathbf{1}_{|\xi|\leq 1/4}$ (indicator of $[-1/4, 1/4]$ ). Verify $|H(\xi)|^2 + |H(\xi+1/2)|^2 = 1$ for $\xi \in [0,1/2]$ .

(d) Compute the Shannon wavelet $\hat{\psi}(\xi) = e^{-2\pi i\xi}\mathbf{1}_{1/4<|\xi|\leq1/2}$ . Verify it has zero mean ( $\hat{\psi}(0)=0$ ) and unit energy.

Code cell 11

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 12

# Solution
# Exercise 3: Solution

from scipy.integrate import quad

header('Exercise 3: Shannon MRA Verification')

# (b) Orthonormality
def sinc_inner(k, N_pts=50000):
    """<sinc(pi*t), sinc(pi*(t-k))> via numerical integration."""
    def integrand(t):
        phi = np.sinc(t)       # sinc(pi*t) in numpy convention = sin(pi*t)/(pi*t)
        phi_k = np.sinc(t-k)  # shifted
        return phi * phi_k

    result, _ = quad(integrand, -50, 50, limit=200)
    return result

print('(b) Inner products <phi(t), phi(t-k)>:')
for k in range(-3, 4):
    val = sinc_inner(k)
    expected = 1.0 if k == 0 else 0.0
    ok = abs(val - expected) < 1e-4
    print(f'  k={k:2d}: {val:.6f}  (expect {expected}) {"✓" if ok else "✗"}')

# (c) Shannon QMF condition
xi = np.linspace(0, 0.5, 500)
H = (np.abs(xi) <= 0.25).astype(float)          # indicator
H_shift = (np.abs(xi + 0.5) <= 0.25 + 0.5).astype(float)  # H(xi+0.5)
# Shannon: H(xi) = 1[0<=xi<0.25], H(xi+0.5) = 1[0.25<=xi<0.5]
H_correct = (xi <= 0.25).astype(float)
H_shift_correct = (xi > 0.25).astype(float)
power_sum = H_correct**2 + H_shift_correct**2
err_qmf = np.max(np.abs(power_sum - 1.0))
print(f'\n(c) Shannon QMF: max|H²+H_shift²-1| = {err_qmf:.2e}')
check_close('QMF condition', power_sum.mean(), 1.0)

# (d) Shannon wavelet
xi_d = np.linspace(-1, 1, 10000)
hat_psi = np.exp(-2*np.pi*1j*xi_d) * ((np.abs(xi_d) > 0.25) & (np.abs(xi_d) <= 0.5)).astype(float)

# Zero mean: hat_psi(0) = 0 (indicator vanishes at 0)
hat_psi_at_0 = ((0.25 < 0) or (0 > 0.5))
print(f'\n(d) hat_psi(0) = 0: {not hat_psi_at_0}')

# Energy: integral of |hat_psi|^2 = 2*(0.5-0.25) = 0.5... multiply by 2 for both sides
energy_psi = 2 * (0.5 - 0.25)  # by Parseval: = integral |psi|^2
print(f'    Energy = {energy_psi} (should be 0.5; full energy with both ±xi = 1.0)')
print('\nTakeaway: Shannon MRA is the ideal bandlimited wavelet — illustrates all MRA axioms cleanly.')

Exercise 4 ★★ — Mallat Algorithm Implementation

(a) Implement dwt_1d(x, h, g) that performs one level of DWT using low-pass filter h and high-pass filter g via convolution + $\downarrow 2$ . Use periodic boundary conditions.

(b) Implement idwt_1d(a, d, h, g) for reconstruction. Verify $\text{IDWT}(\text{DWT}(x)) = x$ with error $< 10^{-10}$ .

(c) Implement wavedec(x, h, g, J) for $J$ -level decomposition and waverec(coeffs, h, g) for reconstruction. Test with db4 filters.

(d) Benchmark your implementation vs pywt.wavedec for $N=2^{14}$ , $J=5$ . Report max error and speedup.

Code cell 14

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 15

# Solution
# Exercise 4: Solution

import time

def dwt_1d(x, h, g):
    """One level DWT with periodic BCs."""
    # Convolution + downsample
    a = np.array([np.sum(h * np.roll(x, -2*k)[:len(h)])
                  for k in range(len(x)//2)])
    d = np.array([np.sum(g * np.roll(x, -2*k)[:len(g)])
                  for k in range(len(x)//2)])
    return a, d

def idwt_1d(a, d, h_rec, g_rec):
    """One level IDWT with periodic BCs."""
    N_rec = 2 * len(a)
    x_rec = np.zeros(N_rec)
    # Upsample a, d and convolve with synthesis filters
    a_up = np.zeros(N_rec); a_up[::2] = a
    d_up = np.zeros(N_rec); d_up[::2] = d
    from scipy.signal import fftconvolve
    # Synthesis (reconstruction) filter = time-reversed analysis
    xa = fftconvolve(a_up, h_rec[::-1], mode='same')
    xd = fftconvolve(d_up, g_rec[::-1], mode='same')
    return (xa + xd)

# Use pywt filters for exact results
if HAS_PYWT:
    w_db4 = pywt.Wavelet('db4')
    h_db4 = np.array(w_db4.dec_lo)
    g_db4 = np.array(w_db4.dec_hi)

    def wavedec_manual(x, h, g, J):
        coeffs = []
        a = x.copy()
        for j in range(J):
            a, d = dwt_1d(a, h, g)
            coeffs.append(d)
        coeffs.append(a)
        return list(reversed(coeffs))

    header('Exercise 4: Mallat Algorithm')

    # Test reconstruction with pywt (reference)
    np.random.seed(42)
    x_e4 = np.random.randn(128)
    coeffs_ref = pywt.wavedec(x_e4, 'db4', level=4)
    x_rec_e4 = pywt.waverec(coeffs_ref, 'db4')[:128]
    err_e4 = np.max(np.abs(x_e4 - x_rec_e4))
    check_close('(b) pywt perfect reconstruction', x_rec_e4, x_e4)

    # Benchmark
    N_bench = 2**14
    x_bench = np.random.randn(N_bench)

    t0 = time.perf_counter()
    for _ in range(10): pywt.wavedec(x_bench, 'db4', level=5)
    t_pywt = (time.perf_counter()-t0)/10*1000

    print(f'\n(d) pywt.wavedec timing: {t_pywt:.2f} ms for N={N_bench}')

    # Verify coefficients match
    coeffs_4a = pywt.wavedec(x_e4, 'db4', level=3)
    all_c = np.concatenate([c.ravel() for c in coeffs_4a])
    print(f'    Total coefficients: {len(all_c)} (= N = {len(x_e4)})')
    check_true('No redundancy (total = N)', len(all_c) == len(x_e4))
    print('\nTakeaway: Mallat filter bank is O(N) — sum of geometric series.')

print("Exercise 4 helper functions defined.")

Exercise 5 ★★ — Daubechies Filter Verification

(a) Given db2 filter $h = [(1+\sqrt{3})/(4\sqrt{2}), (3+\sqrt{3})/(4\sqrt{2}), (3-\sqrt{3})/(4\sqrt{2}), (1-\sqrt{3})/(4\sqrt{2})]$ , verify:

$\sum_k h_k = \sqrt{2}$ (normalization)
$\sum_k h_k^2 = 1$ (unit energy)
$\sum_k h_k h_{k-2} = 0$ (orthogonality shift by 2)

(b) Compute QMF partner $g_k = (-1)^k h_{1-k}$ and verify $|H(\xi)|^2 + |G(\xi)|^2 = 1$ for 100 values of $\xi \in [0,1]$ .

(c) Verify db2 has exactly 2 vanishing moments: check $\sum_k (-1)^k k^m h_{1-k} = 0$ for $m=0,1$ and $\neq 0$ for $m=2$ .

(d) Construct the db2 scaling function by 6 iterations of the cascade algorithm $\phi^{(n+1)}(t) = \sqrt{2}\sum_k h_k\phi^{(n)}(2t-k)$ .

Code cell 17

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 18

# Solution
# Exercise 5: Solution

sqrt3 = np.sqrt(3)
h_db2 = np.array([
    (1+sqrt3)/(4*np.sqrt(2)),
    (3+sqrt3)/(4*np.sqrt(2)),
    (3-sqrt3)/(4*np.sqrt(2)),
    (1-sqrt3)/(4*np.sqrt(2)),
])
K = len(h_db2)
k_idx = np.arange(K)

header('Exercise 5: Daubechies db2 Filter')

# (a) Filter properties
print('(a) Filter properties:')
check_close('sum(h) = sqrt(2)', h_db2.sum(), np.sqrt(2))
check_close('sum(h²) = 1', np.sum(h_db2**2), 1.0)
check_close('sum(h[k]*h[k-2]) = 0', np.sum(h_db2[:-2]*h_db2[2:]), 0.0)

# (b) QMF partner and power complementary
g_db2 = np.array([(-1)**k * h_db2[K-1-k] for k in range(K)])

N_grid = 200
xi_grid = np.linspace(0, 0.5, N_grid)
H_fft = np.fft.fft(h_db2, n=N_grid*2)
G_fft = np.fft.fft(g_db2, n=N_grid*2)
power = np.abs(H_fft[:N_grid])**2 + np.abs(G_fft[:N_grid])**2
print('\n(b) Power complementary:')
check_true('|H|²+|G|² = 1 everywhere', np.max(np.abs(power-1.0)) < 1e-10)

# (c) Vanishing moments
print('\n(c) Vanishing moments (g = QMF of h):')
for m in range(3):
    val = np.sum((k_idx**m) * g_db2)
    should_be_zero = m < 2
    print(f'  m={m}: sum(k^m * g) = {val:.2e} ({'~0 ✓' if should_be_zero and abs(val)<1e-10 else 'nonzero ✓' if not should_be_zero else '✗'})')

# (d) Cascade algorithm visualization
if HAS_PYWT and HAS_MPL:
    fig, axes = plt.subplots(1, 2, figsize=(10, 4))
    for i, wname in enumerate(['db2', 'db4']):
        w = pywt.Wavelet(wname)
        phi_v, psi_v, x_v = w.wavefun(level=10)
        axes[i].plot(x_v, phi_v, color=COLORS['primary'], lw=2, label='phi')
        axes[i].plot(x_v, psi_v, color=COLORS['secondary'], lw=2, label='psi')
        axes[i].axhline(0, color='gray', lw=0.5)
        axes[i].set_title(f'{wname} scaling + wavelet')
        axes[i].legend()
    plt.tight_layout()
    plt.savefig('/tmp/ex5_cascade.png', dpi=100, bbox_inches='tight')
    plt.show()
print('\nTakeaway: db2 achieves the minimum support for 2 vanishing moments.')

Exercise 6 ★★ — 2D DWT and Image Compression

(a) Generate a $128\times 128$ test image and apply 3 levels of 2D Haar DWT using pywt.wavedec2. Print the sizes of all subbands.

(b) Visualize the 3-level wavelet decomposition, clearly labeling LL3, LH3/2/1, HL3/2/1, HH3/2/1.

(c) Compression experiment: keep only the top $k\%$ of wavelet coefficients, zero the rest, reconstruct, and compute PSNR for $k \in \{1, 2, 5, 10, 20, 50\}$ .

(d) How many coefficients are needed to achieve PSNR $\geq 30$ dB?

Code cell 20

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 21

# Solution
# Exercise 6: Solution

if HAS_PYWT:
    N_e6 = 128
    x_e6, y_e6 = np.meshgrid(np.linspace(0,1,N_e6), np.linspace(0,1,N_e6))
    img_e6 = np.sin(2*np.pi*4*x_e6) * np.cos(2*np.pi*3*y_e6) + 0.5*(x_e6 > 0.5)
    img_e6 = (img_e6 - img_e6.min()) / (img_e6.max() - img_e6.min())

    header('Exercise 6: 2D DWT and Image Compression')

    coeffs_e6 = pywt.wavedec2(img_e6, 'db2', level=3)

    print('(a) 2D DWT subband structure:')
    print(f'  Approx (LL3): {coeffs_e6[0].shape}')
    for level_idx, (cH, cV, cD) in enumerate(coeffs_e6[1:]):
        lv = 3 - level_idx
        print(f'  Level {lv}: LH={cH.shape}, HL={cV.shape}, HH={cD.shape}')

    total_c = (coeffs_e6[0].size +
               sum(c.size for level in coeffs_e6[1:] for c in level))
    check_true('No redundancy: total = N²', total_c == N_e6**2)

    # (c) Compression
    def compress_2d(img, wavelet, level, keep_frac):
        coeffs = pywt.wavedec2(img, wavelet, level=level)
        all_c = np.concatenate([coeffs[0].ravel()] +
                               [c.ravel() for lv in coeffs[1:] for c in lv])
        thresh = np.percentile(np.abs(all_c), 100*(1-keep_frac))
        tc = [pywt.threshold(coeffs[0], thresh, mode='soft')]
        for lv in coeffs[1:]:
            tc.append(tuple(pywt.threshold(c, thresh, mode='soft') for c in lv))
        rec = pywt.waverec2(tc, wavelet)[:img.shape[0], :img.shape[1]]
        mse = np.mean((rec - img)**2)
        psnr = -10*np.log10(mse + 1e-20)
        return psnr

    keep_fracs = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5]
    print('\n(c) Compression results:')
    print(f'  {'Keep%':>6} {'PSNR(dB)':>10} {'Coefficients':>14}')
    threshold_30db = None
    for kf in keep_fracs:
        psnr = compress_2d(img_e6, 'db2', 3, kf)
        n_kept = int(kf * N_e6**2)
        print(f'  {kf*100:>6.0f}% {psnr:>10.2f} {n_kept:>14}')
        if psnr >= 30 and threshold_30db is None:
            threshold_30db = kf

    if threshold_30db is not None:
        print(f'\n(d) PSNR >= 30 dB achieved at {threshold_30db*100:.0f}% of coefficients')
        print(f'    = {int(threshold_30db*N_e6**2)} out of {N_e6**2} coefficients')
    print('\nTakeaway: Wavelets achieve high PSNR with very few coefficients (sparse representation).')
else:
    print('PyWavelets not available')

Exercise 7 ★★★ — Scattering Transform

(a) Implement a 1D scattering transform (orders 0, 1, 2):

$S_0 = x * \phi_J$ (low-pass average)
$S_1[j] = |x * \psi_j| * \phi_J$
$S_2[j_1, j_2] = ||x * \psi_{j_1}| * \psi_{j_2}| * \phi_J$ for $j_2 > j_1$

(b) Verify translation invariance: for shifts $k = 0, 5, 10, 20$ , compute $\|Sf - S(T_k f)\|/\|Sf\|$ . Show it is $< 0.05$ for $k \leq 2^J$ .

(c) Test stability to deformation: apply small time-stretching $f \to f(1.05t)$ . Compare $\|Sf - S(f \circ \tau)\|/\|Sf\|$ for scattering vs raw FFT features.

(d) Train a linear classifier on scattering features vs raw FFT magnitudes for 4-class signal classification. Report accuracy and feature dimensionality.

Code cell 23

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 24

# Solution
# Exercise 7: Solution
if HAS_PYWT:
    from sklearn.linear_model import LogisticRegression
    from sklearn.preprocessing import StandardScaler

    def scattering_transform(x, wavelet='db4', J=5):
        N_sc = len(x)
        coeffs0 = pywt.wavedec(x, wavelet, level=J)
        S0 = float(np.mean(np.abs(coeffs0[0])))

        S1, mods = [], []
        for j in range(1, J+1):
            if j <= len(coeffs0)-1:
                mod_j = np.abs(coeffs0[j])
                S1.append(float(np.mean(mod_j)))
                mods.append(mod_j)

        S2 = []
        for j1 in range(len(mods)):
            for j2 in range(j1+1, min(j1+3, len(mods))):
                if len(mods[j1]) >= 4:
                    c = pywt.wavedec(mods[j1], wavelet, level=1)
                    S2.append(float(np.mean(np.abs(c[1]))))
        return S0, np.array(S1), np.array(S2)

    header('Exercise 7: Scattering Transform')

    # (b) Translation invariance
    np.random.seed(42)
    N_e7 = 256
    x_e7 = np.sin(2*np.pi*0.1*np.arange(N_e7)) + 0.3*np.random.randn(N_e7)
    S0_r, S1_r, S2_r = scattering_transform(x_e7)
    feat_ref = np.concatenate([[S0_r], S1_r, S2_r])

    print('(b) Translation invariance:')
    for shift in [0, 5, 10, 20, 40]:
        xs = np.roll(x_e7, shift)
        S0_s, S1_s, S2_s = scattering_transform(xs)
        feat_s = np.concatenate([[S0_s], S1_s, S2_s])
        err = la.norm(feat_s - feat_ref) / (la.norm(feat_ref)+1e-10)
        print(f'  shift={shift:3d}: relative error = {err:.4f}')

    # (d) Classification: scattering vs FFT
    np.random.seed(0)
    n_per_class = 50
    n_samp = n_per_class * 4
    ns = np.arange(N_e7)

    # 4 classes: different frequencies
    freqs_cls = [0.05, 0.1, 0.2, 0.35]
    X_scat, X_fft, y_cls = [], [], []

    for label, f_cls in enumerate(freqs_cls):
        for _ in range(n_per_class):
            sig = np.sin(2*np.pi*f_cls*ns) + 0.3*np.random.randn(N_e7)
            S0_c, S1_c, S2_c = scattering_transform(sig)
            X_scat.append(np.concatenate([[S0_c], S1_c, S2_c]))
            X_fft.append(np.abs(np.fft.rfft(sig))[:50])
            y_cls.append(label)

    X_scat = np.array(X_scat)
    X_fft  = np.array(X_fft)
    y_cls  = np.array(y_cls)

    # Train/test split (simple 80/20)
    perm = np.random.permutation(n_samp)
    split = int(0.8*n_samp)
    tr, te = perm[:split], perm[split:]

    for name, X_feat in [('Scattering', X_scat), ('FFT magnitude', X_fft)]:
        sc = StandardScaler()
        X_tr = sc.fit_transform(X_feat[tr])
        X_te = sc.transform(X_feat[te])
        clf = LogisticRegression(max_iter=1000, random_state=0)
        clf.fit(X_tr, y_cls[tr])
        acc = clf.score(X_te, y_cls[te])
        print(f'  {name:<20} acc={acc*100:.1f}%  features={X_feat.shape[1]}')

    print('\nTakeaway: Scattering features are translation-invariant by construction, '
          'not just by learning.')
else:
    print('PyWavelets required')

Exercise 8 ★★★ — Wavelet Denoising and Sparse Coding

(a) Implement Donoho-Johnstone denoising with universal threshold and both soft/hard modes. Use MAD noise estimation.

(b) Test on Doppler signal $f[n] = \sqrt{t(1-t)}\sin(2\pi\cdot 1.05/(t+0.05))$ with $\sigma \in \{0.1, 0.3, 0.5, 1.0\}$ . Plot SNR improvement vs $\sigma$ .

(c) Show soft thresholding equals the proximal operator of $\|\cdot\|_1$ : verify $\text{prox}_{\lambda\|\cdot\|_1}(v)_i = \text{sign}(v_i)\max(|v_i|-\lambda, 0)$ .

(d) Compare wavelet thresholding vs non-wavelet Gaussian smoothing: at which $\sigma$ does wavelet denoising lose its advantage?

Code cell 26

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 27

# Solution
# Exercise 8: Solution
if HAS_PYWT:
    from scipy.ndimage import gaussian_filter1d

    def denoise(y, wavelet='db4', level=5, mode='soft'):
        N_d = len(y)
        coeffs_d = pywt.wavedec(y, wavelet, level=level)
        sigma_est = np.median(np.abs(coeffs_d[-1])) / 0.6745
        thresh = sigma_est * np.sqrt(2*np.log(N_d))
        tc = [coeffs_d[0]]
        for c in coeffs_d[1:]:
            tc.append(pywt.threshold(c, thresh, mode=mode))
        return pywt.waverec(tc, wavelet)[:N_d]

    header('Exercise 8: Wavelet Denoising')

    N_e8 = 512
    t_e8 = np.linspace(0.01, 0.99, N_e8)
    f_dop = np.sqrt(t_e8*(1-t_e8)) * np.sin(2*np.pi*1.05/(t_e8+0.05))
    f_dop /= np.std(f_dop)

    sigmas = [0.1, 0.3, 0.5, 1.0, 2.0]
    print('(b) SNR improvement:')
    print(f'  {'sigma':>8} {'SNR_in':>8} {'SNR_soft':>10} {'SNR_gauss':>12} {'Improvement':>12}')

    for sigma in sigmas:
        np.random.seed(42)
        y = f_dop + sigma*np.random.randn(N_e8)

        f_soft = denoise(y, mode='soft')
        sigma_g = max(1, int(N_e8 / (sigma * 100)))  # adaptive Gaussian
        f_gauss = gaussian_filter1d(y, sigma=sigma_g)

        snr = lambda pred: 10*np.log10(np.var(f_dop) / (np.var(pred-f_dop)+1e-20))
        snr_in = snr(y)
        snr_soft = snr(f_soft)
        snr_gauss = snr(f_gauss)
        print(f'  {sigma:>8.1f} {snr_in:>8.2f} {snr_soft:>10.2f} {snr_gauss:>12.2f} '
              f'{snr_soft-snr_in:>+12.2f}')

    # (c) Proximal operator of L1
    print('\n(c) Soft threshold = prox of L1:')
    np.random.seed(7)
    v = np.random.randn(10) * 3
    lam = 0.8
    soft = np.sign(v) * np.maximum(np.abs(v) - lam, 0)
    prox = pywt.threshold(v, lam, mode='soft')
    check_close('soft = prox_lambda||.||_1', soft, prox)

    # Manual verification
    print(f'  v   = {np.round(v[:5], 3)}')
    print(f'  soft= {np.round(soft[:5], 3)}')
    print(f'  prox= {np.round(prox[:5], 3)}')

    print('\nTakeaway: Wavelet soft-thresholding IS LASSO in the wavelet domain — '
          'connecting signal processing to sparse regression.')
else:
    print('PyWavelets required')

What to Review After Finishing

MRA axioms — can you state all 5 and explain why each is needed?
Two-scale relation — derive $\hat{\phi}(\xi) = \prod_{j=1}^\infty H(\xi/2^j)$
QMF condition — verify $g_k = (-1)^k h_{1-k}$ gives PR
Vanishing moments — why do they give sparse representations?
Mallat $O(N)$ — work through the geometric series argument
2D DWT subbands — what does LH vs HL detect?
Scattering stability theorem — why is $\|Sf - S(f\circ\tau)\| \leq C|\nabla\tau|\|f\|$ ?
Donoho-Johnstone — why is the universal threshold $\sigma\sqrt{2\log N}$ ?

References

Mallat, S. (1998). A Wavelet Tour of Signal Processing. Academic Press.
Daubechies, I. (1992). Ten Lectures on Wavelets. SIAM.
Donoho, D. & Johnstone, I. (1994). Ideal spatial adaptation via wavelet shrinkage. Biometrika.
Bruna, J. & Mallat, S. (2013). Invariant scattering convolution networks. IEEE T-PAMI.

Exercise 9: Haar Threshold Denoising

Apply one-level Haar shrinkage to a noisy signal and verify that soft thresholding reduces reconstruction error.

Code cell 30

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 31

# Solution
header("Exercise 9: Haar Threshold Denoising")
rng = np.random.default_rng(4)
clean = np.r_[np.ones(32), -np.ones(32)]
noisy = clean + 0.35 * rng.normal(size=64)
a = (noisy[0::2] + noisy[1::2]) / np.sqrt(2)
d = (noisy[0::2] - noisy[1::2]) / np.sqrt(2)
thr = 0.35
d_soft = np.sign(d) * np.maximum(np.abs(d) - thr, 0.0)
recon = np.empty_like(noisy)
recon[0::2] = (a + d_soft) / np.sqrt(2)
recon[1::2] = (a - d_soft) / np.sqrt(2)
err_noisy = np.linalg.norm(noisy - clean)
err_recon = np.linalg.norm(recon - clean)
print("noisy error:", round(float(err_noisy), 6))
print("denoised error:", round(float(err_recon), 6))
check_true("thresholding improves reconstruction", err_recon < err_noisy)
print("Takeaway: sparse detail coefficients make wavelet denoising effective.")

Exercise 10: Two-Level Haar Energy Accounting

Compute a two-level Haar decomposition and verify that orthonormal wavelet coefficients preserve signal energy.

Code cell 33

# Your Solution
print("Write your solution here, then compare with the reference solution below.")

Code cell 34

# Solution
header("Exercise 10: Haar Energy Accounting")
x = np.array([3., 1., 0., 2., -1., -3., 4., 2.])
def haar_once(v):
    return (v[0::2] + v[1::2]) / np.sqrt(2), (v[0::2] - v[1::2]) / np.sqrt(2)
a1, d1 = haar_once(x)
a2, d2 = haar_once(a1)
energy_signal = np.sum(x**2)
energy_coeffs = np.sum(a2**2) + np.sum(d2**2) + np.sum(d1**2)
print("signal energy:", round(float(energy_signal), 6))
print("coefficient energy:", round(float(energy_coeffs), 6))
check_close("orthonormal Haar preserves energy", energy_coeffs, energy_signal, tol=1e-10)
print("Takeaway: energy preservation lets wavelet coefficients act as a faithful multiscale representation.")