Analyse spectrale

Une sinusoïde de pulsation, d'amplitude et de phase inconnues est observée dans du bruit :

$$x[n]=A\cos ({\omega }_{0}n+\phi )+w[n],w[n]\sim \mathcal{N}(0,{\sigma }^2),$$

avec $N=80$ échantillons, $A=1$, $\phi =\pi /4$ et $\sigma =0,5$. Estimer ${\omega }_0$ est le problème central de l'analyse spectrale, celui du radar Doppler, de l'accordeur de guitare et de l'analyse de vibrations. Le chapitre 2 a montré que sa loss est redoutable, multimodale avec un minimum global étroit ; ce projet fait d'abord découvrir la structure du modèle, qui ramène les trois inconnues à une seule, puis fait implémenter tous les algorithmes qui attaquent la loss restante, à la main d'abord, avec scipy ensuite, jusqu'à la stratégie professionnelle, la recherche en deux temps.

1. Lien avec le cours

Le projet réunit les deux pages du chapitre 2 : la base de Fourier du modèle linéaire et sa solution analytique des moindres carrés pour la partie linéaire du modèle, la recherche sur grille et la descente de gradient pour sa partie non linéaire, et la variance Monte-Carlo pour juger le résultat.

2. Structure du code

Le cadre commun des projets d'estimation s'applique : le modèle d'abord, une classe au gabarit de scipy.stats, puis la loss, une fonction séparée, et enfin les algorithmes qui la minimisent, à la main puis avec scipy.

class NoisySinusoid:
    """Sinusoid of unknown pulsation in white Gaussian noise."""

    def __init__(self, N=80, A=1.0, phi=np.pi / 4, sigma=0.5):
        ...

    def rvs(self, omega0, rng):
        """One realization x of length N."""
        ...


def fourier_basis(omega, N):
    """H(omega): two columns, cos(omega n) and sin(omega n)."""
    ...


def loss(omega, x):
    """Concentrated least squares loss J(omega)."""
    ...

Tests imposés. Sans bruit, $J({\omega }_0)\approx 0$ et l'estimateur en deux temps retrouve ${\omega }_0$, $A$ et $\phi$ ; et pour $N$ grand, les deux colonnes de fourier_basis sont quasi orthogonales, ${\mathbf{H}}^T\mathbf{H}\approx (N/2)\mathbf{I}$.

3. Travail demandé

1. La structure du modèle. Démontrer, par la formule d'addition du cosinus, la décomposition

$$A\cos (\omega n+\phi )=a\cos (\omega n)+b\sin (\omega n),$$

en explicitant $(a,b)$ en fonction de $(A,\phi )$, et réciproquement. Conclure sur la nature du modèle : il est non linéaire en $\omega$, mais linéaire en $(a,b)$, un modèle non linéaire à variables linéairement séparables. À $\omega$ fixé, le modèle s'écrit $\mathbf{x}=\mathbf{H}(\omega )\mathit{\theta }+\mathbf{w}$ avec $\mathit{\theta }=(a,b{)}^T$ : donner la matrice $\mathbf{H}(\omega )$, la base de Fourier à deux colonnes du cours. 2. La loss concentrée. Puisque $(a,b)$ se résout analytiquement à $\omega$ fixé, $\hat{\mathit{\theta }}(\omega )=({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathbf{x}$, la recherche non linéaire ne porte plus que sur $\omega$, à travers la loss concentrée

$$J(\omega )=\Vert \mathbf{x}-\mathbf{H}(\omega )\hat{\mathit{\theta }}(\omega ){\Vert }^2.$$

Implémenter $J$, simuler une réalisation avec ${\omega }_0=0,86$ rad/échantillon, et tracer $J(\omega )$ sur $]0,\pi [$. Décrire ce que l'optimiseur devra affronter : combien de minima locaux, quelle largeur pour le creux global ? Trois inconnues au départ, une seule à chercher numériquement : chiffrer ce que la séparation a fait gagner.

3. La grille, à la main. Implémenter la recherche sur grille : $K$ pulsations régulièrement espacées, évaluer $J$ partout, retenir la meilleure, puis en déduire $\hat{a}$, $\hat{b}$, et donc $\hat{A}$ et $\hat{\phi }$. Étudier l'erreur d'estimation en fonction de $K$ ($K=20$, $100$, $1000$, $10000$) et donner le coût en nombre d'évaluations de $J$. Quelle est la précision limite d'une grille à $K$ points ?
4. Le gradient, à la main. Implémenter la descente $\omega \leftarrow \omega -\eta \mathrm{\nabla }J(\omega )$, le gradient étant approché par différence finie centrée, $\mathrm{\nabla }J(\omega )\approx [J(\omega +\delta )-J(\omega -\delta )]/2\delta$. Étudier systématiquement les deux réglages : le pas $\eta$ (trop grand, divergence ; trop petit, lenteur) et l'initialisation (tracer $\hat{\omega }$ final en fonction de $\omega$ initial sur tout $]0,\pi [$ : la carte des bassins d'attraction). Conclure : de quoi la descente de gradient est-elle capable, et incapable ?
5. Les outils scipy. Résoudre le même problème avec scipy.optimize.minimize_scalar puis scipy.optimize.minimize (méthodes "BFGS" et "Nelder-Mead", avec la même initialisation que l'étape 4). Ces optimiseurs professionnels échappent-ils aux minima locaux ? Comparer précision et nombre d'évaluations avec vos implémentations.
6. La recherche en deux temps. Combiner : une grille grossière ($K=50$) pour localiser le bassin global, puis un raffinement local (votre gradient, ou minimize initialisé au meilleur point de grille). À budget d'évaluations égal, comparer la précision de la grille fine seule et de la stratégie en deux temps : chiffrer le gain, et expliquer en deux phrases pourquoi cette division du travail, un algorithme global grossier puis un algorithme local précis, est la stratégie standard.
7. Ouverture. Par Monte-Carlo ($500$ tirages, estimateur en deux temps), estimer la variance de ${\hat{\omega }}_0$ pour $N=20$, $40$, $80$, $160$. Vérifier empiriquement que la variance décroît environ en $1/N^3$, une décroissance bien plus rapide que le $1/N$ de la moyenne empirique : mesurer une fréquence est un problème remarquablement favorable. Remarquer enfin que maximiser ${\mathbf{x}}^T\mathbf{H}({\mathbf{H}}^T\mathbf{H}{)}^{-1}{\mathbf{H}}^T\mathbf{x}$, l'équivalent de minimiser $J$, redonne presque le périodogramme $|\sum _{n}x[n]e^{-\mathrm{i}\omega n}{|}^2$ : la loss concentrée de ce projet est l'outil historique de l'analyse spectrale.

4. Livrables

Le notebook reproductible des modalités communes, avec la démonstration de la séparation (étape 1) rédigée, la carte des bassins d'attraction de l'étape 4, le tableau comparatif précision/coût des étapes 3 à 6, et la courbe variance/$N$ de l'ouverture.