Série 04b-Wavelet-Scattering — bloc A.1 de l’issue #16055. Ce notebook réimplémente la transformée en ondelettes 1D à la main (Haar par paires, puis moteur général par filtrage + sous-échantillonnage), l’utilise pour débruitager par seuillage dur, et confronte le résultat à un passe-bas Fourier sur trois signaux de référence du domaine (Doppler, HeaviSine, Stationnaire — banc canonique Donoho–Johnstone).
Plan. 1. Le problème — pourquoi la base Fourier n’est pas parcimonieuse pour les signaux non stationnaires. 2. Haar à la main — la moyenne et la différence, la plus petite ondelette du monde. 3. Le moteur général — filtrage + sous-échantillonnage, Daubechies-4 par formes closes et constantes publiées, validation croisée contre PyWavelets. 4. Multi-résolution — l’énergie par échelle comme empreinte du signal. 5. Débruitage — seuillage dur/doux, seuil universel sans oracle, contre Fourier avec oracle, sur les trois bancs. 6. Trois exercices — entropie par échelle, seuillage par échelle, localisation de transitoire.
Cap de la série : from scratch — on réimplémente la transformée, on ne consomme pas une bibliothèque haut-niveau sans l’avoir démontée. Le seul usage de pywt (PyWavelets) est la validation croisée : deux implémentations indépendantes doivent s’accorder (écart maximal documenté).
Pédagogie sans image : ce notebook ne génère aucune figure ; toutes les lectures sont ancrées sur des nombres mesurés (erreurs de reconstruction, SNR, énergie par échelle). C’est un choix documenté de la tranche — la dimension visuelle des ondelettes (écriture des coefficients, damiers 2D) arrive avec le bloc A.2 (ondelettes 2D et compression d’image).
See #16055 (bloc A.1).
Débruitager, c’est projeter le signal observé \(y = x + \varepsilon\) sur une base où le signal \(x\) est parcimonieux (peu de coefficients grands) et le bruit uniformément étal (tous les coefficients petits), puis couper les petits. Toute la question est : quelle base ?
La transformée en ondelettes propose un troc : renoncer à la localisation fréquentielle parfaite (une ondelette couvre une bande, pas une fréquence) pour gagner la localisation temporelle (une ondelette vit autour d’une position). Le théorème clé pour ce notebook est plus modeste et plus mesurable : la DWT orthonormale conserve l’énergie exactement (Parseval), donc projeter puis couper des coefficients est une opération bien posée — le seuillage n’est licite que parce que la base est orthonormée.
Ce que nous allons mesurer, honnêtement, sur trois terrains : - Doppler (chirp décroissant) : terrain non stationnaire — l’ondelette doit dominer ; - Stationnaire (deux sinusoïdes + un burst) : terrain de Fourier — le passe-bas doit dominer ; - HeaviSine (sinusoïdes + deux sauts) : terrain mixte — le cas réaliste.
Et une asymétrie de protocole à garder en tête tout du long : le seuillage ondelette utilise un seuil universel \(\lambda = \hat\sigma\sqrt{2\ln N}\) calculé sans aucune connaissance du signal (l’estimateur \(\hat\sigma\) au MAD ne voit que le bruit), tandis que chaque variante Fourier reçoit une fréquence de coupure choisie en connaissant le signal (oracle). L’ondelette joue à armes inégées — contre elle.
import numpy as np
import pywt # UNIQUEMENT pour la validation croisee (allclose) -- jamais pour calculer
N = 2048 # puissance de 2 imposee par le mode periodization multi-niveaux
SIGMA = 0.35 # ecart-type du bruit additif gaussien
print("numpy", np.__version__, "| pywt", pywt.__version__, "| N =", N)numpy 2.2.6 | pywt 1.8.0 | N = 2048Imports justifiés (règle F). numpy est le seul moteur de calcul : convolution, FFT, médiane. pywt (PyWavelets 1.8.0) n’est jamais utilisé pour transformer : il sert de juge de paix — une seconde implémentation, indépendante, contre laquelle notre moteur doit s’accorder coefficient par coefficient. Si les deux s’accordent à \(2\cdot10^{-14}\) près, la nôtre est juste ; si elles divergent, l’une des deux a un défaut de convention (nous en verrons deux, mesurés, dans ce notebook).
Le bruit est gaussien blanc d’écart-type \(\sigma = 0{,}35\) : avec des signaux normalisés à puissance 1, cela place le SNR bruité autour de 9 dB — assez de bruit pour que le débruitage soit un problème, pas assez pour que la tâche soit désespérée.
def stationnaire(n):
"""Deux tons stationnaires (8 et 21 cycles sur la fenetre)."""
tt = np.arange(n) / n
return np.sin(2 * np.pi * 8 * tt) + 0.6 * np.sin(2 * np.pi * 21 * tt)
def doppler(n):
"""Chirp decroissant canonique (Donoho-Johnstone) : sinus de frequence (1+eps)/(eps+t)."""
tt = np.arange(1, n + 1) / n
eps = 0.05
return np.sqrt(n / 8) * np.sin(2 * np.pi * (1 + eps) / (eps + tt))
def heavisine(n):
"""Deux sinusoides + deux echelons (Donoho-Johnstone) : le cas mixte."""
tt = np.arange(1, n + 1) / n
return 4 * np.sin(4 * np.pi * tt) - np.sign(tt - 0.3) - np.sign(0.72 - tt)
# normalisation a puissance 1 (comparabilite des SNR entre bancs)
signals = {}
for name, fn in [("Stationnaire", stationnaire), ("Doppler", doppler), ("HeaviSine", heavisine)]:
s = fn(N).astype(float)
signals[name] = s / np.sqrt(np.mean(s ** 2))
# provenance : HeaviSine doit etre bit-exact contre la reference PyWavelets
h_ref = np.asarray(pywt.data.demo_signal("HeaviSine", N))
h_err = np.max(np.abs(heavisine(N) - h_ref))
d_ref = np.asarray(pywt.data.demo_signal("Doppler", N))
d_cos = float(
(signals["Doppler"] / np.linalg.norm(signals["Doppler"]))
@ (d_ref / np.linalg.norm(d_ref))
)
print(f"provenance HeaviSine vs pywt.data : maxdiff = {h_err:.2e} (bit-exact)")
print(f"provenance Doppler vs pywt.data : similarite cosinus = {d_cos:.4f}")provenance HeaviSine vs pywt.data : maxdiff = 0.00e+00 (bit-exact)
provenance Doppler vs pywt.data : similarite cosinus = 0.9747Lecture — provenance des signaux de test. HeaviSine est bit-exact (maxdiff = 0.0) contre la référence pywt.data.demo_signal : notre formule from scratch reproduit exactement le signal canonique de Donoho–Johnstone. Le Doppler, lui, n’est que la même famille (similarité cosinus \(\approx 0{,}97\)) : la convention de grille de PyWavelets diffère légèrement de la formule Wavelab historique que nous implémentons. C’est une leçon de rigueur plutôt qu’un échec : chaque fois qu’on compare deux implémentations, il faut distinguer erreur (les coefficients diffèrent alors que la convention est la même) et convention (chacun est cohérent, mais les choix de phase/grille diffèrent). Nous allons recontrer exactement cette distinction deux fois dans la suite — une fois sur l’ordre des coefficients de filtre, une fois sur le signe du filtre passe-haut.
Avant le moteur général, la plus petite ondelette du monde. Haar regarde le signal par paires adjacentes \((x_{2k}, x_{2k+1})\) et calcule :
\[a_k = \frac{x_{2k} + x_{2k+1}}{\sqrt{2}} \qquad d_k = \frac{x_{2k} - x_{2k+1}}{\sqrt{2}}\]
Inverser est immédiat : \(x_{2k} = (a_k + d_k)/\sqrt{2}\), \(x_{2k+1} = (a_k - d_k)/\sqrt{2}\). Notez la symétrie du geste : la même paire de formules, les signes échangés.
def dwt_haar(x):
"""DWT de Haar par paires : le geste fondateur, sans boucle sur les taps."""
a = (x[0::2] + x[1::2]) / np.sqrt(2.0)
d = (x[0::2] - x[1::2]) / np.sqrt(2.0)
return a, d
def idwt_haar(a, d):
x = np.empty(2 * len(a))
x[0::2] = (a + d) / np.sqrt(2.0)
x[1::2] = (a - d) / np.sqrt(2.0)
return x
x = signals["Stationnaire"]
a, d = dwt_haar(x)
rec = idwt_haar(a, d)
print(f"longueur : {len(x)} -> approximation {len(a)} + detail {len(d)}")
print(f"reconstruction : erreur max = {np.max(np.abs(rec - x)):.2e}")
print(f"Parseval : |x|^2 = {np.sum(x**2):.6f} |a|^2+|d|^2 = {np.sum(a**2)+np.sum(d**2):.6f}"
f" (ecart relatif {(np.sum(a**2)+np.sum(d**2)-np.sum(x**2))/np.sum(x**2):.2e})")
# demonstration sur 8 echantillons : ce que voient a et d
demo = np.array([3.0, 1.0, 3.0, 3.0, 1.0, 3.0, 3.0, 1.0])
ad, dd = dwt_haar(demo)
print("\npetite demo x =", demo)
print(" a (moyennes x racine2) =", np.round(ad, 4))
print(" d (differences x racine2) =", np.round(dd, 4))longueur : 2048 -> approximation 1024 + detail 1024
reconstruction : erreur max = 6.66e-16
Parseval : |x|^2 = 2048.000000 |a|^2+|d|^2 = 2048.000000 (ecart relatif -1.11e-16)
petite demo x = [3. 1. 3. 3. 1. 3. 3. 1.]
a (moyennes x racine2) = [2.8284 4.2426 2.8284 2.8284]
d (differences x racine2) = [ 1.4142 0. -1.4142 1.4142]Lecture — Haar en trois mesures. 1. Reconstruction exacte à \(7\cdot10^{-16}\) près : la transformation est bijective, aucune information n’est perdue — le seuillage sera donc le seul endroit où nous détruisons quelque chose, et ce sera un choix, pas un accident. 2. Parseval au niveau machine : \(\|a\|^2 + \|d\|^2 = \|x\|^2\) à \(1{,}1\cdot10^{-16}\) relatif près. L’énergie se répartit entre approximation et détail, elle ne se crée ni ne disparaît. Sur ce signal, l’essentiel part dans \(a\) (les deux tons sont lisses à l’échelle de la paire) ; un échelon, lui, produirait deux coefficients \(d\) énormes et localisés. 3. La petite démo montre le geste : là où la paire est constante (3,3), le détail est nul ; là où elle saute (3,1), le détail vaut \(2/\sqrt2 \approx 1{,}41\). L’ondelette détecte le changement — c’est toute son économie.
Ce que Haar ne sait pas faire : être régulière. Une constante par morceaux projette des résidus sur plusieurs niveaux ; pour représenter un signal lisse avec très peu de coefficients, il faut des ondelettes à moments d’annulation — c’est le rôle des Daubechies du §2.
Haar est un cas particulier d’une construction générale. Une ondelette orthonormale est définie par un seul filtre passe-bas \(h\) (le scaling filter, \(n\) coefficients, \(\sum_j h_j = \sqrt2\), \(\sum_j h_j^2 = 1\)) ; le filtre passe-haut d’analyse s’en déduit par miroir alterné :
\[g_j = (-1)^{j+1}\, h_{n-1-j}\]
(nous verrons mesuré pourquoi l’exposant est \(j+1\) et pas \(j\) : les deux choix reconstruisent parfaitement, mais un seul correspond à la convention PyWavelets). Un niveau de transformée, en mode périodisation (le signal est traité comme circulaire — indispensable pour une longueur exactement divisée par 2 à chaque niveau) :
\[a_k = \sum_{j=0}^{n-1} h_j \; x_{(2k + \lfloor n/2 \rfloor - j) \bmod N}, \qquad d_k = \sum_{j=0}^{n-1} g_j \; x_{(2k + \lfloor n/2 \rfloor - j) \bmod N}\]
et la synthèse est l’adjoint exact — la même expression d’indice, mais en accumulant au lieu de lire :
\[x_i = \sum_k a_k\, h_j + d_k\, g_j \quad \text{sur chaque indice } i = (2k + \lfloor n/2 \rfloor - j) \bmod N\]
Ce choix d’implémentation (gather/scatter symétriques) garantit la reconstruction par construction : la matrice de synthèse est littéralement la transposée de la matrice d’analyse. Nous écrivons les deux en boucles transparentes — c’est le prix pédagogique du from scratch : chaque indice est visible, aucune astuce de convolution cachée.
def qmf(lo):
"""Filtre passe-haut d'analyse, convention PyWavelets : g_j = (-1)^(j+1) * h_{n-1-j}."""
return ((-1.0) ** (np.arange(len(lo)) + 1)) * lo[::-1]
def dwt_1d(x, lo):
"""Un niveau de DWT, periodization. a[k] = somme_j lo[j] * x[(2k + n//2 - j) mod N]."""
N, n = len(x), len(lo)
hi = qmf(lo)
a = np.empty(N // 2)
d = np.empty(N // 2)
for k in range(N // 2):
sa = sd = 0.0
for j in range(n):
v = x[(2 * k + n // 2 - j) % N]
sa += lo[j] * v
sd += hi[j] * v
a[k] = sa
d[k] = sd
return a, d
def idwt_1d(a, d, lo):
"""Synthese = adjoint exact de dwt_1d : MEME expression d'indice, on accumule."""
N, n = 2 * len(a), len(lo)
hi = qmf(lo)
x = np.zeros(N)
for k in range(len(a)):
for j in range(n):
x[(2 * k + n // 2 - j) % N] += a[k] * lo[j] + d[k] * hi[j]
return x
# les trois ondelettes de ce notebook
HAAR = np.array([1.0, 1.0]) / np.sqrt(2.0)
print("HAAR :", HAAR, " somme =", HAAR.sum(), " somme^2 =", np.sum(HAAR**2))HAAR : [0.70710678 0.70710678] somme = 1.414213562373095 somme^2 = 0.9999999999999998Pourquoi ce découpage gather/scatter ? L’erreur d’implémentation classique de la DWT maison n’est pas dans les formules mais dans leur appariement : une analyse en phase « corrélation » \((x_{2k+j})\) couplée à une synthèse en phase « convolution » \((x_{2k-j})\) reconstruit… autre chose (notre première version mesurait une erreur de reconstruction de \(3{,}5\) sur un signal d’amplitude \(1{,}7\) — soit une transformation qui ne reconstruit pas du tout). Écrire la synthèse comme l’adjoint littéral de l’analyse — même expression d’indice, lecture remplacée par accumulation — ferme cette classe d’erreur par construction : la matrice de synthèse est la transposée de celle d’analyse, et pour un banc orthonormal, \(A^\top A = I\).
L’offset \(\lfloor n/2 \rfloor\) dans l’indice, discret et peu visible dans la littérature, est la convention de phase de PyWavelets en mode périodisation — nous l’avons déterminée empiriquement (sonde impulsion par impulsion) et elle est validée ci-dessous par l’écart maximal. Pour Haar (\(n=2\)), cet offset est invisible : les deux termes de la somme commutent. C’est précisément pourquoi Haar « marchait toujours » dans nos prototypes et pourquoi les filtres longs ont exposé le problème.
SQ2, SQ3 = np.sqrt(2.0), np.sqrt(3.0)
# Daubechies a 2 moments d'annulation : formes closes exactes (le "D4" des manuels)
D4 = np.array([1 - SQ3, 3 - SQ3, 3 + SQ3, 1 + SQ3]) / (4 * SQ2)
# Daubechies a 4 moments : constantes publiques, ordre PyWavelets
DB4 = np.array([
-0.010597401785069032, 0.0328830116668852, 0.030841381835560764,
-0.18703481171909309, -0.027983769416859854, 0.6308807679298589,
0.7148465705529157, 0.2303778133088965,
])
def wavedec_scratch(x, lo, level):
"""Decomposition multi-niveaux : [a_L, d_L, ..., d_1]."""
coeffs = []
a = x.astype(float).copy()
for _ in range(level):
a, d = dwt_1d(a, lo)
coeffs.append(d)
coeffs.append(a)
return coeffs[::-1]
def waverec_scratch(coeffs, lo):
a = coeffs[0].copy()
for d in coeffs[1:]:
a = idwt_1d(a, d, lo)
return a
# validation croisee : reconstruction + accord coefficient par coefficient avec pywt
x = signals["Stationnaire"]
print(f"{'ondelette':12s} {'taps':>4s} {'niv.max':>7s} {'recon':>10s} {'ecart pywt':>11s}")
for name, lo, wname in [("haar", HAAR, "db1"), ("D4 (db2)", D4, "db2"), ("db4", DB4, "db4")]:
L = pywt.dwt_max_level(N, len(lo))
mine = wavedec_scratch(x, lo, L)
rec = waverec_scratch(mine, lo)
ref = pywt.wavedec(x, wname, level=L, mode="periodization")
cv = max(np.max(np.abs(m - r)) for m, r in zip(mine, ref))
print(f"{name:12s} {len(lo):4d} {L:7d} {np.max(np.abs(rec - x)):10.2e} {cv:11.2e}")ondelette taps niv.max recon ecart pywt
haar 2 11 3.11e-15 1.95e-14
D4 (db2) 4 9 6.00e-15 1.78e-14
db4 8 8 2.00e-15 1.78e-15Lecture — deux implémentations indépendantes s’accordent. Pour les trois ondelettes : reconstruction à l’erreur machine (\(\le 6\cdot10^{-15}\)) et accord coefficient par coefficient avec PyWavelets à \(2\cdot10^{-14}\) près — y compris l’approximation et tous les niveaux de détail. Le moteur from scratch et la bibliothèque industrielle calculent la même transformée.
Deux précisions honnêtes sur ce que cette validation prouve et ne prouve pas. Elle prouve que notre implémentation est conforme à la convention de référence (périodisation, phase, signe du passe-haut). Elle ne prouve pas la théorie (l’existence des ondelettes de Daubechies, les moments d’annulation) — cela relève des mathématiques, pas du test. Et elle n’aurait rien détecté si nous avions choisi l’autre convention de signe : \((−1)^j\) au lieu de \((−1)^{j+1}\) reconstruit aussi parfaitement, mais détail et détail de référence sont alors opposés — c’est en comparant dec_hi directement qu’on a épinglé le signe. Le triplet (test de reconstruction / test d’accord / inspection des filtres) est le garde-fou complet ; chaque test seul laisse passer une classe d’erreur différente.
Sur les deux familles de coefficients : D4 (4 taps) s’écrit en formes closes avec \(\sqrt3\) — aucun chiffre magique, la calculatrice de l’étudiant la retrouve ; db4 (8 taps) n’a pas de forme close utilisable : ce sont des constantes numériques publiées, transcrites ici en 16 chiffres.
La décomposition complète récursive divise l’approximation encore et encore : après \(L\) niveaux, le signal est la superposition de \(L\) détails \(d_1, \dots, d_L\) (des bandes de fréquences de plus en plus basses, chacune localisée en temps) et d’une approximation résiduelle \(a_L\). Grâce à l’orthonormalité, l’énergie totale est la somme exacte des énergies par échelle — et le profil de cette répartition est une empreinte du signal : où vit son contenu, à quelle échelle ?
Mesurons cette empreinte sur nos trois signaux de référence, avec db4 à profondeur maximale.
L = pywt.dwt_max_level(N, len(DB4))
print(f"profondeur maximale db4 pour N={N} : L = {L} niveaux\n")
print(f"{'echelle':8s} {'Stationnaire':>14s} {'Doppler':>10s} {'HeaviSine':>11s}")
profils = {}
for name in ["Stationnaire", "Doppler", "HeaviSine"]:
coeffs = wavedec_scratch(signals[name], DB4, L)
tot = sum(np.sum(c ** 2) for c in coeffs)
profils[name] = [100 * np.sum(c ** 2) / tot for c in coeffs]
labels = [f"a{L}"] + [f"d{L - i}" for i in range(L)]
for i, lab in enumerate(labels):
row = f"{lab:8s}"
for name in ["Stationnaire", "Doppler", "HeaviSine"]:
row += f" {profils[name][i]:13.1f}%"
print(row)
print()
for name in ["Stationnaire", "Doppler", "HeaviSine"]:
p = profils[name]
top3 = sorted(p[1:], reverse=True)[:3]
print(f"{name:13s}: top-3 details = {sum(top3):5.1f}% (a{L} seul = {p[0]:.1f}%)")profondeur maximale db4 pour N=2048 : L = 8 niveaux
echelle Stationnaire Doppler HeaviSine
a8 0.0% 52.9% 97.1%
d8 38.2% 16.1% 2.4%
d7 36.5% 10.9% 0.2%
d6 23.5% 5.8% 0.1%
d5 1.7% 5.2% 0.1%
d4 0.0% 4.2% 0.0%
d3 0.0% 3.0% 0.0%
d2 0.0% 1.7% 0.0%
d1 0.0% 0.1% 0.0%
Stationnaire : top-3 details = 98.3% (a8 seul = 0.0%)
Doppler : top-3 details = 32.9% (a8 seul = 52.9%)
HeaviSine : top-3 details = 2.8% (a8 seul = 97.1%)Lecture — trois empreintes, trois structures temporelles. (les pourcentages exacts sont ceux de la sortie ci-dessus)
Cette lecture prépare le débruitage : le seuillage coupe les petits coefficients à toutes les échelles ; il préservera ce qui est concentré (les tons, les sauts) et combattra le bruit qui, lui, s’étale uniformément. Et le cas Doppler explique pourquoi Fourier va souffrir : son information est partagée entre les échelles, donc entre les bandes — couper une bande, c’est l’amputer.
Estimateur ondelette (sans oracle). On décompose \(y\), on estime le bruit par le MAD du détail le plus fin (\(\hat\sigma = \mathrm{MAD}(d_1)/0{,}6745\) — robuste car un signal parcimonieux y dépose peu d’énergie), et on coupe tous les détails sous le seuil universel \(\lambda = \hat\sigma\sqrt{2 \ln N}\). L’intuition : \(N\) coefficients de bruit gaussien ont un maximum espéré \(\approx \sigma\sqrt{2\ln N}\), donc ce seuil devrait tuer presque tout le bruit en épargnant les coefficients porteurs de signal. Le seuil dur garde les coefficients intacts ; le seuil doux les rétracte de \(\lambda\) (biais, mais pas de seuil de décision brutal).
Estimateur Fourier (avec oracle). On applique un passe-bas idéal à une coupure \(f_c\). Nous en mesurons trois variantes, toutes avec connaissance du signal propre : étroit (\(f_c\) sous la plus haute fréquence des tons — sacrifie le large bande), garde-tout (la plus petite \(f_c\) retenant 99,9 % de l’énergie propre), et libre (la \(f_c\) qui maximise le SNR global — la plus généreuse possible envers Fourier). Chaque variante incarne un choix de ce qu’on accepte de détruire ; l’ondelette, elle, ne choisit pas : un seul seuil, aucune coupure, aucune connaissance du signal.
Le SNR local (fenêtre autour du transitoire), le pic du burst reconstruit et l’amputation d’énergie propre complètent le SNR global — car nous allons voir qu’un SNR global peut récompenser une distorsion.
def shrink(coeffs, threshold, soft=False):
"""Seuillage des coefficients de detail (l'approximation est toujours gardee)."""
out = [coeffs[0]]
for c in coeffs[1:]:
if soft:
out.append(np.sign(c) * np.maximum(np.abs(c) - threshold, 0.0))
else:
out.append(np.where(np.abs(c) > threshold, c, 0.0))
return out
def sigma_hat_mad(coeffs):
d1 = coeffs[-1]
return np.median(np.abs(d1 - np.median(d1))) / 0.6745
def fourier_lowpass(y, fc):
Yf = np.fft.rfft(y)
freqs = np.fft.rfftfreq(len(y))
return np.fft.irfft(Yf * (freqs <= fc), n=len(y))
def snr(ref, est):
return 10 * np.log10(np.sum(ref ** 2) / np.sum((ref - est) ** 2))
def bench(name, clean, seed, burst_pos=None, burst_len=16, burst_amp=1.2):
"""Un banc de debruitage complet : ondelette dur/doux vs trois Fouriers."""
rng = np.random.default_rng(seed)
if burst_pos is not None:
b = np.zeros(N)
b[burst_pos:burst_pos + burst_len] = burst_amp
clean = clean + b
y = clean + SIGMA * rng.standard_normal(N)
w = slice(burst_pos - 50, burst_pos + burst_len + 64) if burst_pos else slice(0, N)
coeffs = wavedec_scratch(y, DB4, L)
lam = sigma_hat_mad(coeffs) * np.sqrt(2 * np.log(N))
xh_hard = waverec_scratch(shrink(coeffs, lam, False), DB4)
xh_soft = waverec_scratch(shrink(coeffs, lam, True), DB4)
kept = sum(np.count_nonzero(c) for c in shrink(coeffs, lam, False)[1:])
Cf = np.fft.rfft(clean)
freqs = np.fft.rfftfreq(N)
tot_f = np.sum(np.abs(Cf) ** 2)
cum = np.cumsum(np.abs(Cf) ** 2) / tot_f
fc_keep = freqs[int(np.searchsorted(cum, 0.999)) + 1]
xk = fourier_lowpass(y, fc_keep)
xn = fourier_lowpass(y, 0.05) # etroit : garde les tons, sacrifie le large-bande
best = (-np.inf, None)
for fc in np.linspace(0.005, 0.5, 200):
s = snr(clean, fourier_lowpass(y, fc))
if s > best[0]:
best = (s, fc)
snr_free, fc_free = best
amp = np.sum(np.abs(Cf[freqs > fc_free]) ** 2) / tot_f
res = {
"bruite": snr(clean, y), "bruite_local": snr(clean[w], y[w]),
"dur": snr(clean, xh_hard), "dur_local": snr(clean[w], xh_hard[w]),
"doux": snr(clean, xh_soft), "doux_local": snr(clean[w], xh_soft[w]),
"garde_tout": snr(clean, xk), "garde_tout_local": snr(clean[w], xk[w]),
"etroit": snr(clean, xn), "etroit_local": snr(clean[w], xn[w]),
"fc_keep": fc_keep, "libre": snr_free, "fc_free": fc_free,
"amputation": 100 * amp, "lam": lam, "gardes": kept,
"total_det": N - len(coeffs[0]),
"xh": xh_hard, "clean": clean, "y": y, "w": w,
}
if burst_pos is not None:
res["pic_dur"] = xh_hard[burst_pos:burst_pos + burst_len].max()
res["pic_clean"] = clean[burst_pos:burst_pos + burst_len].max()
res["pic_ft"] = xk[burst_pos:burst_pos + burst_len].max()
res["pic_narrow"] = xn[burst_pos:burst_pos + burst_len].max()
res["pic_doux"] = xh_soft[burst_pos:burst_pos + burst_len].max()
return res
R = {}
R["Doppler"] = bench("Doppler", signals["Doppler"], seed=101)
R["HeaviSine"] = bench("HeaviSine", signals["HeaviSine"], seed=202)
R["Stationnaire"] = bench("Stationnaire", signals["Stationnaire"], seed=303, burst_pos=1000)
for k, r in R.items():
print(f"{k:13s} bruite {r['bruite']:6.2f} dB | dur {r['dur']:6.2f} (local {r['dur_local']:6.2f})"
f" | doux {r['doux']:6.2f} | F garde-tout {r['garde_tout']:6.2f} @fc={r['fc_keep']:.4f}"
f" | F etroit {r['etroit']:6.2f} | F libre {r['libre']:6.2f} @fc={r['fc_free']:.4f}"
f" (amput {r['amputation']:.2f}%) | lam={r['lam']:.3f}, gardes {r['gardes']}/{r['total_det']}")Doppler bruite 8.97 dB | dur 17.83 (local 17.83) | doux 12.28 | F garde-tout 12.68 @fc=0.1997 | F etroit 11.95 | F libre 13.36 @fc=0.1120 (amput 1.70%) | lam=1.369, gardes 51/2040
HeaviSine bruite 9.03 dB | dur 22.43 (local 22.43) | doux 19.07 | F garde-tout 21.64 @fc=0.0195 | F etroit 18.94 | F libre 23.28 @fc=0.0100 (amput 0.19%) | lam=1.455, gardes 8/2040
Stationnaire bruite 9.38 dB | dur 16.37 (local 11.44) | doux 11.73 | F garde-tout 17.84 @fc=0.0776 | F etroit 19.92 | F libre 21.19 @fc=0.0249 (amput 0.37%) | lam=1.378, gardes 57/2040Lecture — Doppler : le terrain non stationnaire, victoire nette de l’ondelette (sans oracle). Le seuillage dur atteint 17,83 dB (+8,86 dB sur le bruité), contre 13,36 dB pour le meilleur Fourier — le cutoff libre, choisi en connaissant le signal — soit +4,47 dB d’avance sans oracle. Toutes les variantes Fourier s’effondrent ensemble (11,95 à 13,36 dB) parce que le problème est structurel : le §3 a montré que l’énergie du chirp est étalée sur toutes les échelles (top-3 détails = 32,9 % seulement). Un passe-bas doit donc choisir ce qu’il ampute : le cutoff libre sacrifie 1,70 % de l’énergie propre — le début du chirp, précisément la partie haute fréquence — et le cutoff étroit (11,95 dB) ampute plus fort encore tout en gardant du bruit dans sa bande.
Côté ondelette, la parcimonie est extrême : 51 coefficients de détail sur 2040 (2,5 %) survivent au seuil \(\lambda = 1{,}369\). Un chirp dont l’information occupe tout le spectre se représente par 51 nombres, parce que chaque portion temporelle du signal trouve son échelle — c’est exactement le troc temps-fréquence du §0.
Lecture — HeaviSine : le terrain mixte, quasi-égalité, et une lecture honnête. Le Fourier libre (23,28 dB @ fc=0,0100) devance ici le seuillage dur (22,43 dB) de 0,85 dB. Trois observations pour ne pas se mentir :
Verdict honnête : quasi-égalité sur un terrain qui favorise Fourier, avec un estimateur qui n’a rien su du signal.
Lecture — Stationnaire + burst : le terrain de Fourier, assumé — et le contrepoint du transitoire. Le passe-bas gagne nettement en SNR global : étroit 19,92 dB, libre 21,19 dB (amputation 0,37 %) contre 16,37 dB pour l’ondelette. C’est le résultat attendu — deux tons éternels sont l’objet parfait de Fourier, et couper à fc = 0,025–0,05 élimine le bruit partout ailleurs. L’asymétrie de protocole joue aussi : ces cutoffs savent où vivent les tons.
Le contrepoint est le pic du burst (valeur propre : 1,362). Toutes les méthodes le surestiment — le maximum d’un signal bruité est biaisé vers le haut — mais pas autant : dur 1,485 (écart +0,12), étroit 1,496 (+0,13), garde-tout 1,564 (+0,20). L’ondelette dure est la plus proche du vrai pic, le garde-tout la plus déformée. Surtout, le doux écrase le transitoire : pic 0,996, soit −27 % — les coefficients du burst, tout juste au-dessus du seuil, sont rétractés de \(\lambda\). Et le SNR local de l’ondelette dure (autour du burst) tombe à 11,44 dB : le transitoire reste le point dur des deux familles.
# tableau recapitulatif des trois bancs
print(f"{'banc':14s} {'bruite':>7s} {'db4 dur':>8s} {'db4 doux':>9s} {'F/tout':>7s} {'F/etroit':>8s} {'F/libre':>8s} {'amput%':>7s}")
for k, r in R.items():
print(f"{k:14s} {r['bruite']:7.2f} {r['dur']:8.2f} {r['doux']:9.2f} {r['garde_tout']:7.2f}"
f" {r['etroit']:8.2f} {r['libre']:8.2f} {r['amputation']:7.2f}")
print()
r = R["Stationnaire"]
print(f"burst (Stationnaire) : pic propre {r['pic_clean']:.3f} | dur {r['pic_dur']:.3f}"
f" | doux {r['pic_doux']:.3f} | Fourier etroit {r['pic_narrow']:.3f} | garde-tout {r['pic_ft']:.3f}")
print(f"gains vs bruite (dB) : " + " | ".join(
f"{k}: dur {r0['dur']-r0['bruite']:+.2f}" for k, r0 in R.items()))banc bruite db4 dur db4 doux F/tout F/etroit F/libre amput%
Doppler 8.97 17.83 12.28 12.68 11.95 13.36 1.70
HeaviSine 9.03 22.43 19.07 21.64 18.94 23.28 0.19
Stationnaire 9.38 16.37 11.73 17.84 19.92 21.19 0.37
burst (Stationnaire) : pic propre 1.362 | dur 1.485 | doux 0.996 | Fourier etroit 1.496 | garde-tout 1.564
gains vs bruite (dB) : Doppler: dur +8.86 | HeaviSine: dur +13.40 | Stationnaire: dur +6.99Lecture — dur contre doux : le biais, mesuré. L’écart est systématique et fort : −5,55 dB (Doppler), −3,36 dB (HeaviSine), −4,64 dB (Stationnaire). Deux mécanismes, l’un global, l’autre local :
C’est l’arbitrage biais/variance en action : doux = biais, dur = variance résiduelle. Les raffinements de la littérature (SureShrink, BayesShrink) attaquent précisément ce point en adaptant \(\lambda\) par échelle — l’exercice 2 en construit la version pédagogique.
| banc | db4 dur (sans oracle) | meilleur Fourier (avec oracle) | amputation propre | verdict |
|---|---|---|---|---|
| Doppler | 17,83 dB | 13,36 dB (libre) | 1,70 % | ondelette +4,47 dB |
| HeaviSine | 22,43 dB | 23,28 dB (libre) | 0,19 % | quasi-égalité, terrain Fourier |
| Stationnaire | 16,37 dB | 21,19 dB (libre) | 0,37 % | Fourier +4,82 dB |
Aucune base n’est universellement parcimonieuse — c’est la leçon centrale, mesurée trois fois. La base ondelette gagne là où le signal est étalé en fréquence mais localisé en temps (chirp) ; la base Fourier gagne là où le signal est stationnaire et à bande étroite ; le cas mixte rend le match serré. Choisir une base, c’est écrire une hypothèse sur la structure du signal dans la méthode elle-même.
Et l’asymétrie de protocole rend le score encore plus net qu’il n’y paraît : chaque Fourier cité a reçu sa coupure en connaissant le signal propre (l’amputation mesure ce que cette clairvoyance détruit : 0,19 à 1,70 %), tandis que l’ondelette n’a vu que le signal bruité — \(\lambda\) vient du MAD de \(d_1\), aucune autre information. Le seuil universel est d’ailleurs conservateur : 0,4 à 2,8 % des coefficients de détail gardés. On peut faire mieux (exercice 2) ; on ne peut pas faire plus honnête.
Contexte. Le §3 a mesuré l’énergie par échelle. L’énergie est dominée par les grands coefficients ; l’information se mesure mieux en entropie de Shannon : \(H_\ell = -\sum_k p_{k}^{(\ell)} \log_2 p_k^{(\ell)}\) où \(p_k^{(\ell)} = c_k^2 / \sum_{\ell'} \|c^{(\ell')}\|^2\) est la part d’énergie du coefficient \(k\) de l’échelle \(\ell\). Une échelle qui porte peu d’énergie mais beaucoup de coefficients non nuls a une entropie élevée : c’est du contenu distribué — typiquement le bruit. C’est aussi le fondement de la compression (codage entropique des coefficients) que le bloc A.2 exploitera en 2D.
Objectif. Compléter entropie_par_echelle ci-dessous pour qu’elle retourne la liste des \(H_\ell\) (approximation incluse), puis comparer les trois signaux : quelle échelle du signal Stationnaire porte le plus d’information alors qu’elle porte peu d’énergie ?
Indices. - # Etape 1 : normaliser par l’énergie totale des coefficients (toutes échelles) ; - # Etape 2 : attention au logarithme de zéro — ne sommer que sur les coefficients non nuls ; - # Etape 3 : comparer au profil d’énergie du §3 — l’ordre des échelles est-il le même ?
def entropie_par_echelle(coeffs):
"""Entropie de Shannon (bits) de la distribution d'energie de CHAQUE echelle.
coeffs : [a_L, d_L, ..., d_1] (sortie de wavedec_scratch)
retour : liste de len(coeffs) entropies, une par echelle
"""
# Etape 1 : energie totale sur toutes les echelles
# TODO etudiant
# Etape 2 : pour chaque echelle, p_k = c_k^2 / energie_totale, H = -somme p log2 p (p > 0)
# TODO etudiant
# Etape 3 : retourner la liste des H
# TODO etudiant
return None # TODO etudiant
# verification attendue : sur le signal Stationnaire, l'entropie totale
# (somme des H) doit etre FINIE et l'echelle dominante doit diffeRer du profil d'energie
# H_stationnaire = entropie_par_echelle(wavedec_scratch(signals['Stationnaire'], DB4, L))
# print(H_stationnaire)
print("Exercice a completer")Exercice a completerContexte. Le seuil universel \(\lambda = \hat\sigma\sqrt{2\ln N}\) utilise le même seuil pour toutes les échelles, avec un \(\hat\sigma\) estimé sur le seul détail le plus fin. C’est conservateur : au §4, il ne gardait que 0,4 à 2,8 % des coefficients de détail. Un raffinement standard (par exemple dans les estimateurs par échelle de la littérature) estime \(\hat\sigma_\ell\) séparément à chaque échelle \(\ell\) (le bruit blanc reste blanc à toutes les échelles, mais l’estimation locale gagne en robustesse quand le signal dépose beaucoup d’énergie à certaines échelles) et applique \(\lambda_\ell = \hat\sigma_\ell \sqrt{2 \ln N_\ell}\) où \(N_\ell\) est le nombre de coefficients de l’échelle.
Objectif. Implémenter shrink_par_echelle, débruitager le banc Doppler avec, et comparer au seuil global du §4 : le SNR global monte-t-il ? Le SNR local (autour des hautes fréquences du début du chirp) change-t-il dans le même sens ?
Indices. - # Etape 1 : la longueur de l’échelle \(\ell\) est len(coeffs[i]) — le facteur \(\sqrt{2\ln N_\ell}\) diminue aux échelles grossières ; - # Etape 2 : estimer \(\hat\sigma_\ell\) par MAD sur cette échelle ; - # Etape 3 : réutiliser snr et le protocole du §4 (même bruit : default_rng(101)).
def shrink_par_echelle(coeffs, soft=False):
"""Seuillage avec UN seuil par echelle, estime sur l'echelle elle-meme.
coeffs : [a_L, d_L, ..., d_1] ; l'approximation (coeffs[0]) est toujours gardee.
retour : coefficients seuilles, meme structure
"""
out = [coeffs[0]]
for c in coeffs[1:]:
# Etape 1 : N_ech = len(c), facteur sqrt(2 ln N_ech)
# TODO etudiant
# Etape 2 : sigma_ech par MAD local, lam_ech = sigma_ech * facteur
# TODO etudiant
# Etape 3 : seuiller (dur ou doux) a lam_ech
# TODO etudiant
out.append(c) # ligne provisoire : a remplacer par le coefficient seuille
return out
# protocole : meme bruit que le banc Doppler du §4
# rng = np.random.default_rng(101)
# y = signals["Doppler"] + SIGMA * rng.standard_normal(N)
# coeffs = wavedec_scratch(y, DB4, L)
# xh = waverec_scratch(shrink_par_echelle(coeffs), DB4)
# print("SNR seuil global :", snr(signals["Doppler"], R["Doppler"]["xh"]))
# print("SNR seuil/echelle :", snr(signals["Doppler"], xh))
print("Exercice a completer")Exercice a completerContexte. Au §4, le burst du banc Stationnaire était connu (position 1000). En pratique on ne le connaît pas : on veut détecter où il est. Les détails fins \(d_1, d_2\) d’une DWT répondent fort et localement aux discontinuités — c’est la contrepartie temps-fréquence que Fourier n’offre pas (une détection par transformée de Fourier nécessite une fenêtre glissante, c’est-à-dire… une décomposition temps-fréquence).
Objectif. Écrire localise_burst(y) qui retourne la position estimée du transitoire d’un signal \(y\) bruité, en n’utilisant que l’énergie locale des détails fins : reconstruire l’énergie glissante de \(d_1\) (et éventuellement \(d_2\)), chercher son maximum. Valider sur le banc Stationnaire : l’erreur de localisation doit être de l’ordre de la dizaine d’échantillons.
Indices. - # Etape 1 : un coefficient \(d_1\) se calcule sur deux échantillons adjacents — la position temporelle du coefficient \(k\) est environ \(2k\) ; - # Etape 2 : sommer \(|d_1|^2\) dans une fenêtre glissante (par exemple np.convolve sur les carrés) pour lisser le bruit ; - # Etape 3 : np.argmax de l’énergie glissante, convertir l’indice coefficient → indice échantillon, comparer à 1000.
def localise_burst(y):
"""Position estimee du transitoire d'un signal bruite, par energie locale des details fins.
y : signal observe (longueur N) ; retour : indice echantillon estime.
"""
# Etape 1 : d1 (et d2 au besoin) par notre moteur, a la bonne profondeur
# TODO etudiant
# Etape 2 : energie glissante des carres (fenetre ~ 32 coefficients)
# TODO etudiant
# Etape 3 : argmax puis conversion coefficient -> echantillon
# TODO etudiant
return None # TODO etudiant
# validation attendue (bruit du banc Stationnaire) :
# rng = np.random.default_rng(303)
# b = np.zeros(N); b[1000:1016] = 1.2
# y = signals["Stationnaire"] + b + SIGMA * rng.standard_normal(N)
# pos = localise_burst(y)
# print("position estimee :", pos, "| erreur :", abs(pos - 1000))
print("Exercice a completer")Exercice a completerCe que ce notebook a établi, par la mesure. 1. La DWT orthonormale se réimplémente en une trentaine de lignes lisibles (gather/scatter adjoints), reconstruit à l’erreur machine, et s’accorde coefficient par coefficient avec la bibliothèque de référence — à condition d’épingler deux conventions (l’ordre des coefficients publiés, le signe du miroir alterné) que le test de reconstruction seul ne détecte pas. 2. Le profil d’énergie par échelle est une empreinte lisible : tons concentrés, chirp étalé, sauts dilués sur toutes les échelles. 3. Le débruitage par seuillage dur réalise sans aucun oracle ce qu’un passe-bas Fourier ne réalise qu’avec oracle sur le terrain non stationnaire (Doppler), et perd sur le terrain stationnaire — il n’existe pas de base universellement parcimonieuse, et le choix de base encode une hypothèse sur la structure du signal. 4. Le seuillage doux paie son biais (−3,4 à −5,6 dB mesurés) ; le seuil universel est conservateur — les raffinements (par échelle : exercice 2 ; SURE, Bayesian shrinkage : littérature) sont le pont vers le bloc B.
La suite de l’épic #16055. - A.2 — Ondelettes 2D : extension séparable, damiers de détails, compression d’image (le codage entropique de l’exercice 1 devient central) ; - A.3 — Scattering : cascade ondelette → module → moyenne, invariance par translation garantie par construction ; - B — SOTA : PyWavelets pour le débruitage industriel (BayesShrink, SureShrink), kymatio pour le Scattering GPU, et le tableau comparatif from scratch vs SOTA.
Voir README.md de la série pour l’état d’avancement.