Calculer une TFD avec rfft : fréquences des cases, amplitudes, nombre entier de périodes

1. Le programme :

import numpy as np
fe = 1000.0
for N in (1000, 1024):
    t = np.arange(N) / fe
    x = 1.0 + 2.0*np.cos(2*np.pi*50*t) + 0.5*np.cos(2*np.pi*120*t + 0.7)
    X = np.fft.rfft(x)
    f = np.fft.rfftfreq(N, 1/fe)
    A = 2*np.abs(X)/N
    A[0] = np.abs(X[0])/N          # case 0 : valeur moyenne
    for k in np.argsort(A)[::-1][:5]:
        print(N, f[k], round(A[k], 3), round(np.angle(X[k]), 3))

Pour $N=1000$, rfft renvoie $501$ valeurs ; la case $k$ est à la fréquence $k{f}_e/N=k$ hertz, de $0$ à $500Hz=f_e/2$. La durée d'acquisition est $T_{\text{acq}}=N/f_e=1s$ et la résolution $\delta f=1/T_{\text{acq}}=1Hz$.

2. Le programme affiche $A_0=1,000$, $A_{50}=2,000$, $A_{120}=0,500$, avec $\arg X_{120}=0,700$ ; toutes les autres cases valent moins de $10^{-14}$ (erreurs d'arrondi). La fenêtre d'une seconde contient exactement $50$ périodes de la vibration et $120$ périodes de la composante d'engrènement : chaque sinusoïde tombe sur une case. Écrite $\frac{c}{2}({\mathrm{e}}^{j\theta }+{\mathrm{e}}^{-j\theta })$, une sinusoïde d'amplitude $c$ partage son poids entre la case $k$ et la case $N-k$ (non renvoyée par rfft) ; la somme sur les $N$ échantillons donne $X_k=\frac{N}{2}c{\mathrm{e}}^{j\phi }$, d'où le facteur $2/N$. La composante continue n'a pas de double : $X_0=N{c}_0$, d'où $1/N$.

3. Avec $N=1024$, l'écart entre cases devient $\delta f=1000/1024=0,977Hz$, et $50Hz=51,2\delta f$ ne tombe plus sur une case. Le programme lit $1,87V$ à $49,80Hz$ (case $51$), soit $6,4\%$ de moins que la vraie amplitude, et des valeurs non nulles sur les cases voisines : $0,47V$ à $50,78Hz$, $0,31V$ à $48,83Hz$, puis $0,21$ et $0,17V$ (figure). De même, $120Hz$ est lu $0,486V$ à $120,12Hz$, et la valeur moyenne $1,006V$. C'est la fuite spectrale : la fenêtre de $1,024s$ contient $51,2$ périodes du $50Hz$ ; la TFD, qui traite la fenêtre comme une période, voit un signal périodisé présentant un raccord brutal, dont le spectre s'étale. L'erreur sur la case la plus proche se prévoit (résultat admis) : pour un décalage de $0,2$ case, l'amplitude lue est multipliée par $\frac{\sin (0,2\pi )}{0,2\pi }=0,936$, soit $1,87V$.

4. Pour une lecture exacte, la fenêtre doit contenir un nombre entier de périodes de chaque composante : $50T_{\text{acq}}$ et $120T_{\text{acq}}$ entiers, soit $T_{\text{acq}}$ multiple de $0,1s$, et $N$ multiple de $100$. Quand les fréquences ne sont pas connues d'avance, on allonge l'acquisition (la fuite reste, mais les cases se resserrent) ; un $N$ « rond » pour l'algorithme n'a aucune raison d'être rond pour le signal.

Ce que l'exercice installe : la case $k$ est à $k{f}_e/N$ ; on lit $X_0/N$ pour la moyenne et $2\left|X_k\right|/N$ pour une amplitude ; ces lectures ne sont exactes que si la fenêtre contient un nombre entier de périodes, faute de quoi la raie fuit sur les cases voisines.