Lab IV : Identification des paramètres d'une sinusoïde de fréquence connue

• Cours "Physique du Numérique" - Portail René Descartes - AMU

Préparé par :

• Jean-Marc Themlin (v. 2021-09), Aix-Marseille Université © Contenus à diffusion restreinte, dans le cadre de ce cours.

Remarque

A toutes fins utiles, on rappelle ici qu'une solution de secours (moins pratique, et probablement moins stable) permet d'accéder aux énoncés sur des notebooks fonctionnels sur un serveur distant (Binder) uniquement depuis votre navigateur. Vous pourrez exécuter, modifier et compléter les cellules de code et remplir les cellules de texte ; à l'issue de votre travail, vous pourrez importer (commande Download dans le menu File) le notebook complété pour le stocker localement sur votre PC et le déposer sur AMeTICE. Pour accéder à Binder pour travailler ce premier TP, il suffit de cliquer sur l'icône ci-dessous et de sélectionner dans le répertoire des fichiers accessibles, le notebook Enonce_Lab_4_Identification_Sinusoides.ipynb ; soyez patient, ça peut prendre quelques minutes pour lancer l'environnement :

Introduction

Le produit scalaire, qui permet de quantifier la ressemblance entre deux signaux numériques $x[n]$ et $y[n]$ définis sur le même domaine, est défini par : $$<s,r>=\sum_{n=1}^{N} s[n] \ r^*[n] \ \Delta t$$ où $\Delta t$ est la période d'échantillonnage commune aux deux signaux numériques.

Dans ce TP, vous allez utiliser le produit scalaire pour déterminer (identifier) l'amplitude $A$ et la phase à l'origine $\phi$ d'une (co)sinusoïde $s(t)=A \ \cos (2\pi f_0 t + \phi)$ de fréquence $f_0$ connue, qu'on supposera échantillonnée sur l'intervalle $T_0$ entre les temps $t_1$ et $t_2$ par la (co)sinusoïde discrète $s[n] = A \ \cos (2\pi f_0 n \Delta t +\phi)$, avec une période d'échantillonnage $\Delta t$. Le signal (co)sinusoïdal d'entrée à analyser est donc représenté par N échantillons de la sinusoïde numérique $s[n]$, avec $N = \frac{t_2-t_1}{\Delta t}=\frac{T_0}{\Delta t}$.

Nous avons montré au cours l'intérêt d'utiliser le produit scalaire entre $s[n]$ et une exponentielle complexe de référence $r[n]=e^{j2\pi f_0t}$ (elle fixera notamment l'instant $t=0$ qui détermine le déphasage recherché) à la fréquence $f_0$, la même que celle de $s[n]$.

Avec l'hypothèse que la durée du signal est un multiple de $T_0=\frac{1}{f_0}$, le résultat du produit scalaire entre $s[n]$ et $r[n]$, qui est un nombre complexe, s'exprime en effet : $$\left\langle {s,r} \right\rangle \simeq {A \over 2}\;e^{j\phi } \;(t_2 - t_1 ) = {A \over 2}\;e^{j\phi } \;N\Delta t$$ Par conséquent, le calcul numérique du produit scalaire fournit directement les paramètres estimés : ParseError: KaTeX parse error: Undefined control sequence: \eqalign at position 2: \̲e̲q̲a̲l̲i̲g̲n̲{ & A_{Est} \…

Programmation : Quelques fonctions Python utiles

Exécutez Help(nom_de_la_fonction) ou nom_de_la_fonction? pour afficher l'aide des fonctions ci-dessous et comprendre ce qu'elles réalisent.

Dans la bibliothèque standard : len(x), abs(z) où z est un nombre complexe qui s'écrit z = a*1j+b ou z = c exp(1j*d).

Dans la bibliothèque math : floor(x)

Dans la bibliothèque numpy : np.size(t), np.zeros(N), np.angle(z) où z est un nombre complexe.

Dans la bibliothèque matplotlib.pyplot : plt.subplot(2,1,n=1,2)

Boucles for : Examinez l'exemple ci-dessous, qui contient deux exemples de boucles for. Dans le second, on remplit séquentiellement une liste avec un élément calculé dans la boucle, ce qui nous sera très utile par la suite.

Problème IV-1 : Paramètres d'une sinusoïde

Utilisez la cellule de code ci-dessous pour générer une (co)sinusoïde discrète $s[n]$

Q1: A partir du graphe, déterminez visuellement les paramètres $(A,f_0,\phi)$ de cette (co)sinusoïde d'expression générale $=A \ \cos (2\pi f_0 t + \phi)$.

Q1-R: Votre réponse ici : A = 25 f0 = 1/6 Phi = -2

Q2: En utilisant le vecteur $t$ comme support temporel, vérifiez la validité de votre estimation en superposant au graphe précédent la sinusoïde reconstituée à partir de vos estimations

Q2-R: Votre réponse ici : déterminer les paramètres graphiquement n'est jamais assez précis.

Problème IV-2 : Estimer l'amplitude et la phase d'un signal sinusoïdal de fréquence connue $f_0$

Utilisez la cellule de code ci-dessous pour construire une fonction Python apc pour calculer l'amplitue $A$ et la phase à l'origine $\phi$ d'un signal (co)sinusoïdal discret $s[n]$ de fréquence connue $f_0$.

Les paramètres d'entrée de la fonction seront :

• le vecteur support temporel du signal

• le signal à analyser (dont vous allez estimer les paramètres)

• la fréquence $f_0$ de ce signal

Les paramètres retournés par la fonction seront :

• $A_{Est}$, l'amplitude estimée de la sinusoïde

• $\phi_{Est}$, la phase à l'origine estimée de la sinusoïde

Test de la fonction $apc$

Test 1

Dans la cellule ci-dessous, utilisez le signal généré au IV-1 pour vérifier vos estimations initiales.

Q3 : Remplissez la table ci-dessous avec les résultats de l'estimation visuelle ($cf$ IV-1) et de l'estimation numérique.

Dans la cellule ci-dessous, pour vérifier la validité de l'estimation, superposez dans un graphe la sinusoïde initiale ($cf$ IV-1) et la (co)sinusoïde générée à partir des paramètres estimés numériquement (avec apc).

Test 2

Vous allez à présent tester le bon fonctionnement de votre fonction apc sur un signal de paramètres inconnus, hormis bien sûr sa fréquence $f_0$. Les premières instructions de la cellule ci-dessous vont télécharger dans l'environnement le fichier lab_data.mat, qui contient les variables suivantes :

• le vecteur support t_samp

• le vecteur signal s_samp, qui contient les échantillons d'une sinusoïde discrète d'une fréquence $f_0=100 \ Hz$.

Les valeurs retournées par votre fonction seront les suivants :

• la période d'échantillonnage $\Delta t$

• le nombre total $N$ de points que contiennent chacun de ces deux vecteurs

• la valeur de l'amplitude estimée $A_{Est}$

• la valeur de la phase à l'origine estimée $\phi_{Est}$

Pour vérifier la validité de l'estimation, complétez la cellule de code ci-dessous avec les instructions Python permettant de superposer dans un graphe la sinusoïde initiale (s_samp en fonction de t_samp) avec la (co)sinusoïde générée à partir des paramètres estimés numériquement.

Problème IV-3 : Tests de robustesse de la fonction d'estimation apc

Pour qu'un programme informatique soit opérationnel, il doit être testé de façon très approfondie lors de son développement. Ceci implique notamment d'explorer la plus large zone possible de l'espace des paramètres d'entrée du programme.

Cas non-idéal : Durée du signal différente d'un multiple de $T_0=\frac{1}{f_0}$

Les formules qui donnent les paramètres estimés ont été établies dans l'hypothèse que la durée du signal est un multiple (entier) de la période du signal $s[n]$, soit $T_0=\frac{1}{f_0}$.

Mais que se passe-t'il si vous appliquez la fonction apc sur un signal sinusoïdal qui ne contient pas un nombre entier de périodes $T_0=\frac{1}{f_0}$ ?

Dans la cellule de code ci-dessous, générez le signal indiqué :

Q4 : Que vaut la période de $s[n]$ ? Combien de périodes (chiffre décimal !) contient le signal $s[n]$ ?

Q4-R (votre réponse ici) : la période s[n] est 3. le signal s[n] contient D/T = 8/3 = 2,67périodes.

Q5 : Que valent les paramètres $(A,f_0,\phi)$ de cette (co)sinusoïde d'expression générale $=A \ \cos (2\pi f_0 t + \phi)$ ?

Q5-R (votre réponse ici) : d'après la fonction apc A = 7.802284479558461 phi = -1.5250792646745612 f0 = 1/3

Complétez la cellule de code ci-dessus pour estimer avec apc l'amplitude $A_{Est}$ et la phase à l'origine $\phi_{Est}$ de $s[n]$.

Complétez la table ci-dessous en y incluant les écarts absolus ($|A_{Est}-A|$ et $|\phi_{Est}-\phi|$) et relatifs ($\frac{|A_{Est}-A|}{A}$ et $|\frac{\phi_{Est}-\phi}{\phi}|$) :

8 | environ 0,2 | | | Phase à l'origine | -1.5250792646745612 | -pi/2 | | |

Visualiser la dépendance de l'estimation en fonction de la durée de la sinusoïde

Il est instructif de rechercher les tendances de l'évolution des valeurs estimées par apc en fonction de la durée D du signal sinusoïdal. Le but est de produire une figure montrant comment les résultats de l'estimation évoluent avec la longueur du signal.

A cet effet, dans la cellule de code ci-dessous, l'idée est d'initialiser un vecteur dur (un numpy.array de dimension 1) contenant des valeurs croissantes de la longueur D du signal et de son support temporel.

Ensuite, on utilisera une boucle for d'indice n pour construire des vecteurs t et s de longueur croissante, t variant entre 0 et une durée variable contenue dans Dur[n]. Dans la même boucle, on applique l'estimateur apc et on stocke les valeurs estimées dans deux vecteurs distincts, qui contiendront toutes les valeurs estimées de $A_{est}$ et de $\phi_{Est}$, qu'on pourra tracer ensuite dans une figure à deux graphes.

Pour fixer les idées, vous pouvez utiliser les paramètres suivants :

• Paramètres de la sinusoïde : $(A,f_0,\phi)=(3,1,-\pi/2)$

• Un vecteur D contenant des valeurs croissantes entre 0.25 et $D_{Max}=10 s$

Le résultat devrait faire apparaître des oscillations ainsi qu'une décroissance asymptotique de l'erreur des estimations fournies par apc. En vous basant sur l'examen de la figure, répondez aux questions suivantes :

Q6 : Pour quelles valeurs (relatives à $f_0$) les estimations de la phase et de l'amplitude égalent-elles les valeurs vraies des paramètres ? De ce fait, proposez une modification de la condition sur la durée (qui doit être égale à un multiple entier de la période $T_0=\frac{1}{f_0}$) pour obtenir une estimation exacte.

Q6-R (votre réponse ici) : Les multiples de 1 / f0 sont égalent au vrais valeursD = n*T0

Q7 : A partir de quelle durée l'estimation de l'amplitude diffère-t'elle de moins de 1% de la valeur vraie ?

Indication : L'instruction plt.ylim(A,A*1.01) vous permet de restreindre l'échelle en y entre les deux valeurs spécifiées.

Q7-R (votre réponse ici) : A partir de 9 secondes

Q8 : On considérant l'allure des courbes, pouvez-vous dire si on on a en moyenne plus de chances d'obtenir une valeur estimée en excès (supérieure à la valeur vraie) ou une valeur estimée par défaut (inférieure à la valeur vraie) ?

Q8-R (votre réponse ici) : Nous avons plus de chances d'obtenir une valeur estimée en excès

Effet d'une constante ajoutée à $s[n]$

• Dans la cellule ci-dessous, générez quelques périodes d'une sinusoïde discrète d'amplitude 1 à laquelle vous ajouterez une constante réelle $C$ (positive ou négative). Affichez le résultat dans une figure.

• Testez le comportement de votre fonction apc sur ce signal $x[n]=C+s[n]$. Vous testerez plusieurs valeurs de la constante pour valider vos observations et votre conclusion.

Q9 : Dans quelle mesure l'estimation réalisée par $apc$ est-elle affectée par l'addition d'une valeur constante ?

Q9-R (votre réponse ici) :

IV-4 : Compléments : Autres tests de robustesse de la fonction d'estimation apc

Vous n'aurez probablement pas le temps d'effectuer ces exercices de compléments lors de la séance de TP en présentiel. Nous vous conseillons pourtant vivement d'essayer de les résoudre par vous-même, car ils sont très instructifs, et qu'ils inspireront certainement les questions de l'examen de TP terminal.

Effet d'une incertitude sur la fréquence $f_0$

La fréquence $f_0$ de la sinusoïde est connue avec une certaine précision. Dans cette dernière partie, que vous travaillerez en devoir, le but est de tester la validité de l'estimation de l'amplitude et de la phase lorsque la fréquence supposée connue $f_t$, celle que l'on fournit à apc, est légèrement différente de la véritable $f_0$.

Dans la cellule ci-dessous, à l'aide d'une boucle for, générez un graphe qui permet de visualiser la valeur de l'erreur (absolue ou relative) sur les paramètres estimés $A_{est}$ et $\phi_{Est}$ en fonction de la valeur de la fréquence $f_t$ fournie à apc ; on maintiendra $f_t$ à $\pm 25\%$ de sa valeur vraie $f_0$.

Q10 : Dans quelle plage de fréquence l'erreur relative sur l'estimation de l'amplitude reste-t'elle inférieure à 40% ?

Q10-R (votre réponse ici) :

Effet d'un bruit additif gaussien

Si apc devait être utilisée pour estimer phase et amplitude de signaux réels (issus de capteurs ou d'un canal de communication), il y a de fortes chances que ces signaux soient affectés d'un bruit plus ou moins intense.

On peut très facilement simuler un bruit additif gaussien en ajoutant à chaque échantillon du signal des valeurs aléatoires réparties selon une loi normale, i.e. une gaussienne de moyenne $0$ et de variance $\sigma^2=1$. En Python, l'instruction nois=np.random.randn(N) génère un vecteur nois comprenant N valeurs (pseudo-)aléatoires distribuées selon la loi normale. On admettra que la puissance de ce bruit est donnée par la variance $\sigma^2$. Pour obtenir un bruit additif gaussien de puissance $P_N$, il suffira de multiplier nois par la valeur $P_N$.

Le rapport signal-sur-bruit (en abrégé SNR pour Signal-to-Noise ratio) quantifie en décibels l'amplitude relative du signal utile et du bruit. Il s'exprime ici : $$SNR_{dB}= 10 \log_{10}\frac{P_{Signal}}{P_{Noise}} = 10 \log_{10}\frac{A²}{σ²}$$ où $A$ est l'amplitude de la sinusoïde $s[n]$.

• Dans la cellule ci-dessous, générez quelques périodes d'une sinusoïde discrète d'amplitude 1 affectée d'un bruit additif gaussien de puissance $P_N=1$, et affichez le résultat dans une figure. Répétez le processus pour les puissances croissantes $P_N={2,5,20}$, et affichez le résultat dans une unique figure (Indication : utilisez plt.subplot).

• Testez le comportement de votre fonction apc sur une sinusoïde d'amplitude fixée à $1$ et une puissance de bruit prenant les valeurs successives $P_N={1,2,5,20}$. Stockez le résultat des estimations successives de l'amplitude et de la phase dans deux vecteurs, et montrez dans une seconde figure leur évolution en fonction de $P_N$.

Conclusions personnelles

Indiquez ci-dessous le temps approximatif que vous passé à travailler ce TP en-dehors de la séance.

Q11-R (votre réponse ici) : 5h

Ecrivez ci-dessous, en guise de conclusions, quelques phrases décrivant ce que ce TP vous a appris.

Q12 - R : Conclusions personnelles (votre réponse ici) : je n'ai pas bien compris