Analyse spectrale d'un signal numérique

1. Introduction

Ce document montre comment faire l'analyse spectrale des signaux numériques, c'est-à-dire des signaux échantillonnés.

L'analyse spectrale d'un signal repose sur la transformée de Fourier.

Les fonctions de fenêtrage utilisées pour l'analyse spectrale sont étudiées dans Fenêtres pour l'analyse spectrale.

2. Transformée de Fourier

Soit u(t) un signal, c'est-à-dire une fonction du temps à valeurs réelles. La transformée de Fourier de u(t) est une fonction de la fréquence définie par :

$$U\left(f\right)={\int }_{-\infty }^{+\infty }u\left(t\right)exp(-j2\pi ft)dt\text{(1)}$$

La fonction u(t) étant à valeurs réelles, on a :

$$U(-f)=U{\left(f\right)}^{*}\text{(2)}$$

où ${\left[z\right]}^{*}$ désigne le complexe conjugé de $z$.

Soit s(t) un signal dont on souhaite faire l'analyse spectrale. Ce signal a une durée T. On considère par convention qu'il commence à t=-T/2 et qu'il se termine à t=T/2. L'analyse spectrale consiste à calculer la transformée de Fourier de ce signal prolongé sur R. Il faudra éventuellement modifier le signal par une fonction de fenêtrage. La fonction dont on considère la transformée de Fourier est :

$$u\left(t\right)=s\left(t\right)w\left(t\right)\text{(3)}$$

w(t) est la fonction de fenêtrage (ou plus simplement fenêtre), qui est nulle en dehors de l'intervalle [-T/2,T/2].

La fenêtre la plus simple est la fenêtre rectangulaire, définie par :

$$\begin{array}{cc}\hfill & w_r\left(t\right)=1pourt\in [-\frac{T}{2},\frac{T}{2}]\hfill \\ \hfill & w_r\left(t\right)=0pour\left|t\right|>\frac{T}{2}\hfill \end{array}$$

Considérons l'exemple important suivant : un signal s(t) sinusoïdal et une fenêtre rectangulaire. On pose :

$$s\left(t\right)=S_{0}cos(2\pi f_{1}t+\varphi )\text{(4)}$$

La transformée de Fourier de ce signal multiplié par la fenêtre rectangulaire de durée T s'écrit :

$$U\left(f\right)={\int }_{-\frac{T}{2}}^{\frac{T}{2}}{S}_{0}cos(2\pi f_{1}t+\varphi )exp(-j2\pi ft)dt\text{(5)}$$

L'utilisation de la formule d'Euler permet d'écrire :

$$U\left(f\right)=\frac{S_0}{2}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}(exp\left(j(2\pi (f_1-f)t+\varphi )\right)+exp\left(j(-2\pi (f_1+f)t-\varphi )\right))dt\text{(6)}$$

Le calcul conduit à :

$$U\left(f\right)=\frac{S_0}{2}T{e}^{j\varphi }(sinc\left((f_1-f)T\right)+sinc\left((f_1+f)T\right))\text{(7)}$$

où sinc est la fonction sinus cardinale, définie par :

$$sinc\left(u\right)=\frac{sin\left(\pi u\right)}{\pi u}\text{(8)}$$

La transformée de Fourier de la fonction sinusoïdale contient deux termes similaires. Le premier est un sinus cardinal centré en f=f₁, le second un sinus cardinal centré en f=-f₁ (identique au premier).

Nous pouvons faire une représentation graphique du module de U(f). Par convention, on prend la période de s(t) égale à 1, ce qui revient à considérer que cette période est l'unité de temps. La durée T est donc en fait le rapport de la durée du signal par sa période. Remarque : la notation T est souvent utilisée pour la période d'un signal périodique alors qu'ici il s'agit de la largeur de la fenêtre; ce choix est volontaire et sera justifié plus loin.

import numpy as np
from matplotlib.pyplot import *
f1 = 1
T = 10
f = np.linspace(-2,2,100000)
U = 0.5*T*(np.sinc((f-f1)*T)+np.sinc((f+f1)*T))
figure(figsize=(16,6))
plot(f,np.absolute(U))
grid()
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U(f)|/S_0$",fontsize=16)

              

fig1.pdf

Voici le module de la transformée de Fourier pour une fenêtre rectangulaire 10 fois plus longue (soit 100 fois la période) :

T = 100
U = 0.5*T*(np.sinc((f-f1)*T)+np.sinc((f+f1)*T))
figure(figsize=(16,6))
plot(f,np.absolute(U))
grid()
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U(f)|/S_0$",fontsize=16)
              

fig2.pdf

Si l'on s'intéresse au lobe principal de la fonction sinus cardinal, celui-ci peut être défini par $u\in [-1,1]$. La largeur du lobe principal (pour une fenêtre rectangulaire) est donc :

$$\Delta f=\frac{2}{T}\text{(9)}$$

Comme il est démontré dans Fenêtres pour l'analyse spectrale, d'autres fonctions de fenêtrage peuvent avoir une largeur de lobe principale légèrement différente mais elle est toujours inversement proportionnelle à la durée T. Le lobe principal du sinus cardinal centré en f₁ peut être considéré comme la raie associée à la fréquence f₁. Alors que sa largeur est inversement proportionnelle à T, sa hauteur est proportionnelle à T. La limite $T\to \infty$ correspond au cas théorique de la transformée de Fourier d'une fonction périodique sinusoïdale. La largeur de la raie tend vers zéro et sa hauteur tend vers l'infini. La fonction limite est une distribution, la distribution de Dirac centrée en f=f₁.

En pratique, ce qui nous intéresse dans l'analyse spectrale du signal s(t) est sa fréquence et son amplitude (égale à 1 dans le calcul ci-dessus). On représente donc plutôt le module de la transformée de Fourier divisé par T (transformée de Fourier normalisée) :

T = 100
Un = 0.5*(np.sinc((f-f1)*T)+np.sinc((f+f1)*T))
figure(figsize=(16,6))
plot(f,np.absolute(Un))
grid()
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U(f)|/(TS_0)$",fontsize=16)
              

fig3.pdf

La fonction de fenêtrage rectangulaire comporte des lobes secondaires particulièrement intenses. Cela n'est pas nécessairement gênant, pourvu qu'on soit bien conscient que ces lobes viennent de la fenêtre et non pas du phénomène physique étudié. Dans certains cas, il peut être nécessaire d'atténuer ces lobes. Cette opération est nommée apodisation et se fait en changeant la fonction de fenêtrage w(t). La fenêtre de Blackman, étudiée dans Fenêtres pour l'analyse spectrale, a une TF dont des lobes secondaires sont très faibles. Voici le spectre obtenu avec cette fenêtre :

def Wb(f1,T,f):
    x = (f-f1)*T
    return 0.42*np.sinc(x)+0.5/2*(np.sinc(1-x)+np.sinc(1+x))+0.08/2*(np.sinc((2-x))+np.sinc((2+x)))
Un = 0.5*(Wb(f1,T,f) + Wb(-f1,T,f))
figure(figsize=(16,6))
plot(f,np.absolute(Un))
grid()
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U(f)|/(TS_0)$",fontsize=16)
                 

fig4.pdf

La transformée de Fourier (TF) est le produit de convolution de la TF du signal s(t) par la TF de la fenêtre w(t). Or le maximum de cette dernière est inférieur à 1. La valeur du maximum pour la fenêtre de Blackman est 0,42. Il faut donc diviser le module par cette valeur (en plus de la division par T) :

figure(figsize=(16,6))
plot(f,np.absolute(Un)/0.42)
grid()
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U(f)|/(TS_0)$",fontsize=16)
                 

fig5.pdf

L'inconvénient de la fenêtre de Blackman est la grande largeur de son lobe principal (plus de 2 fois la largeur du lobe de la fenêtre rectangulaire). En revanche, les lobes secondaires sont négligeables. L'intérêt pratique de cette fenêtre dépend du signal étudié. Si les raies sont assez espacées, c'est certainement une très bonne fenêtre. Si l'espacement des raies est inférieur à la largeur du lobe, il suffit d'augmenter la durée T (ce qui bien sûr n'est pas toujours possible). Voici le résultat avec une durée 2 fois plus grande :

T = 200
Un = 0.5*(Wb(f1,T,f) + Wb(-f1,T,f))
figure(figsize=(16,6))
plot(f,np.absolute(Un)/0.42)
grid()
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U(f)|/(TS_0)$",fontsize=16)
                 

fig6.pdf

3. Transformée de Fourier à temps discret

Un signal échantillonné avec la période d'échantillonnage T_e est défini par :

$$\begin{array}{ccc}\hfill & t_q=q{T}_e\hfill & \hfill \text{(10)}\\ \hfill & u_q=u\left(t_q\right)\hfill & \hfill \text{(11)}\end{array}$$

où q est un nombre entier (relatif). Le signal u_q est qualifié de signal numérique ou bien de signal à temps discret. L'approximation de l'intégrale $\mathrm{(1)}$ par la formule des rectangles s'écrit :

$$U\left(f\right)\approx \sum _{q=-\infty }^{\infty }{T}_e{u}_{q}exp(-j2\pi fq{T}_e)\text{(12)}$$

La transformée de Fourier à temps discret (TFTD) [1] est définie par :

$$U_{T_e}\left(f\right)=\sum _{q=-\infty }^{\infty }{u}_{q}exp(-j2\pi fq{T}_e)\text{(13)}$$

Dans le domaine du traitement numérique du signal, il est d'usage de poser $Z=exp\left(j2\pi f{T}_e\right)$. La TFTD s'écrit alors :

$$U_{T_e}\left(Z\right)=\sum _{q=-\infty }^{\infty }{u}_q{Z}^{-q}\text{(14)}$$

Cette relation définit la transformée en Z du signal échantillonné.

La TFTD constitue une approximation de la TF, d'autant meilleure que la période d'échantillonnage est petite :

$$\underset{T_e\to 0}{lim}{T}_e{U}_{T_e}\left(f\right)=U\left(f\right)\text{(15)}$$

Notons f_e=1/T_e la fréquence d'échantillonnage. Le théorème de Shannon nous permet d'énoncer deux conditions pour que la TFTD corresponde à la TF :

• La transformée de Fourier U(f) doit être à support borné : il doit exister une fréquence f_m telle que pour |f|>f_m on ait U(f)=0.
• La fréquence d'échantillonnage doit vérifier f_e>2f_m (condition de Nyquist-Shannon).

En pratique, ces deux conditions ne sont pas toujours vérifiées, mais idéalement elles doivent l'être. Une conséquence de ces conditions est que la TFTD ne peut être utilisée pour calculer la TF que pour des fréquences appartenant à l'intervalle [-f_e/2,f_e/2].

La TFTD possède une propriété remarquable (et loin d'être a priori évidente) : si les deux conditions ci-dessus sont vérifiées, la TFTD est exactement égale à la TF sur l'intervalle [-f_e/2,f_e/2]. Cette propriété est une conséquence du théorème de Shannon, qui établit qu'un signal peut être entièrement reconstruit à partir de ses échantillons, pourvu que les deux conditions précédentes soient vérifiées.

La relation $\mathrm{(13)}$ permet en fait de définir une fonction sur R (fréquence réelle quelconque). La TFTD ainsi définie est périodique :

$$U_{T_e}(f+\frac{1}{T_e})=U_{T_e}\left(f\right)\text{(16)}$$

Autrement dit, la période de la TFTD est égale à la fréquence d'échantillonnage.

La figure suivante montre le module de la TF d'un signal à bande de fréquence limitée (dont la TF est à support borné) et le spectre de la TFTD correspondante, dans le cas où la fréquence d'échantillonnage vérifie la condition de Nyquist-Shannon.

Figure pleine page

La définition de la TFTD fait intervenir une somme infinie mais, puisque le signal u_q obtenu par numérisation est non nul seulement sur l'intervalle [-T/2,T/2], cette somme contient en fait un nombre fini de termes. Soit N le nombre d'échantillons sur cet intervalle de temps, indexés par $k=0,1,\dots N-1$ :

$$t_k=-\frac{T}{2}+k\frac{T}{N}\text{(17)}$$ $$U_{Te}\left(f\right)=\sum _{k=0}^{N-1}{T}_e{u}_{k}exp(-j2\pi fk{T}_e)exp\left(j2\pi f\frac{N{T}_e}{2}\right)\text{(18)}$$

Sachant que $T=N{T}_e$, on obtient :

$$U_{Te}\left(f\right)=U_{T,N}\left(f\right)=\frac{T}{N}exp\left(j\pi fT\right)\sum _{k=0}^{N-1}{u}_{k}exp(-j2\pi fk\frac{T}{N})\text{(19)}$$

Nous changeons au passage la notation de la TFTD car cette nouvelle expression fait intervenir T et N à la place de T_e.

4. Transformée de Fourier discrète

Pour calculer numériquement la TFTD, il faut échantillonner l'intervalle de fréquence [-f_e/2,f_e/2]. En raison de la périodicité de la TFTD, on peut aussi considérer que la fréquence appartient à l'intervalle [0,f_e] et échantillonner cet intervalle. Dans un premier temps, on choisit les N fréquences suivantes :

$$f_n=\frac{n}{T}(n=0,1,\dots N-1)\text{(20)}$$

Le fait de choisir un nombre d'échantillons pour la fréquence égal au nombre d'échantillons pour le temps peut sembler a priori une simple convention. Pour ces fréquences, la TFTD s'écrit :

$$U_{T,N}\left(f_n\right)=\frac{T}{N}exp\left(j\pi n\right)\sum _{k=0}^{N-1}{u}_{k}exp(-j2\pi k\frac{n}{N})\text{(21)}$$

La transformée de Fourier discrète est définie par :

$$U_n^{tfd}=\sum _{k=0}^{N-1}{u}_{k}exp(-j2\pi k\frac{n}{N})\text{(22)}$$

Il est important de faire la distinction entre la TFTD et la TFD :

• La transformée de Fourier à temps discret (TFTD), définie par $\mathrm{(13)}$, est définie pour toute fréquence réelle, ce qui correspond au fait que le signal échantillonné n'est pas limité dans le temps.
• La transformée de Fourier discrète (TFD), définie par $\mathrm{(22)}$, est définie seulement pour les fréquences $f_n=\frac{n}{T}$. Le signal est limité dans le temps (de durée T).

Le signal dont on cherche à calculer la TFTD est bien limité dans le temps à cause du fenêtrage (dans le sens où son support est borné) donc la somme qui définie la TFTD a un nombre de termes fini. Dans ce cas, la distinction importante entre dans la TFTD est celle-ci : dans la TFTD seul le temps est discret, alors que dans la TFD, le temps et la fréquence sont discrets. Numériquement, on calcule toujours une TFD puisqu'il est évidemment nécessaire de discrétiser la fréquence. La figure suivante montre la TF d'un signal à bande de fréquence limitée, sa TFTD et sa TFD pour N = 8.

Figure pleine page

La figure ci-dessus peut laisser croire que la TFTD peut être reconstruite à partir de la TFD par interpolation linéaire, ce qui n'est pas vrai en général. Nous verrons plus loin comment la TFTD peut être recontruite à partir de la TFD part la méthode du remplissage par des zéros.

La transformée de Fourier discrète permet aussi de calculer les coefficients de Fourier d'une fonction périodique de période T (la durée du signal) qui serait prolongée par périodicité. Ce lien entre les coefficients de Fourier et la TFD est développé plus loin et dans Introduction à l'analyse spectrale. La fréquence f_n est alors la fréquence de l'harmonique de rang n de ce signal périodique. La TFD est d'ailleurs le plus souvent introduite comme un moyen de calculer les coefficients de Fourier [2]. La TFD donne même les valeurs exactes des coefficients de Fourier si la série s'arrête à un rang fini et si N est assez grand pour que la condition de Nyquist-Shannon soit vérifiée.

Pour les fréquences f_n, nous obtenons la valeur de la TFTD :

$$U_{T,N}\left(f_n\right)=\frac{T}{N}exp\left(j\pi n\right)U_n^{tfd}\text{(23)}$$

qui finalement nous donne la TF pour ces fréquences (si la condition de Nyquist-Shannon est vérifiée) :

$$U\left(f_n\right)=\frac{T}{N}exp\left(j\pi n\right)U_n^{tfd}\text{(24)}$$

La TFTD est définie pour toute fréquence dans l'intervalle [-f_e/2,f_e/2] mais la TFD ne permet de la calculer que pour les fréquences f_n (soit N fréquences dans cet intervalle). Nous verrons plus loin qu'il est nécessaire d'obtenir la TFTD pour un plus grande nombre de fréquences dans cet intervalle afin d'obtenir une meilleure représentation de la TF (c.a.d. un meilleur échantillonnage). Le calcul de la TFTD pour ces fréquences supplémentaires est obtenu par la méthode du remplissage par des zéros, qui sera développée plus loin.

La TFTD définit le spectre d'un signal échantillonné de durée infinie (même si son support est borné à cause du fenêtrage). La TFD définit le spectre d'un signal échantillonné de durée T. Les fréquences de ce spectre sont $f_n=\frac{n}{T}$ où n est un nombre entier.

Tout comme la TFTD, la TFD est périodique :

$$U_{n+N}^{tfd}=U_n^{tfd}\text{(25)}$$

Il en résulte que :

$$U_{-n}^{tfd}=U_{N-n}^{tfd}\text{(26)}$$

Autrement dit, la valeur de la TFD pour la fréquence -f_n est la même que pour la fréquence f_N-n. Supposons pour simplifier l'écriture que N soit pair. Les termes de la TFD de rangs $n=0,1,\dots N/2$ permettent de calculer la TF pour les fréquences $0,\frac{1}{T},\frac{2}{T}\dots \frac{NT}{2}$ (cette dernière fréquence est la moitié de la fréquence d'échantillonnage). Les termes $n=N/2+1,\dots N-1$ permettent de calculer la TF pour les fréquences $-(\frac{N}{2}+1)\frac{1}{T},\dots -\frac{1}{T}$.

Par ailleurs, si la fonction u(t) est à valeurs réelles, on a :

$$U_{N-n}^{tfd}={\left[U_n^{tfd}\right]}^{*}\text{(27)}$$

ce qui est en accord avec le fait que :

$$U(-f)=U^{*}\left(f\right)\text{(28)}$$

Il peut arriver qu'on ait à calculer la TF d'une fonction à variables complexes (notamment dans les calculs de conception des filtres numériques) mais dans ce document u(t) désigne un signal, donc une fonction à valeurs réelles.

Voyons comment la TFTD peut être calculée pour des fréquences autres que les fréquences multiples de 1/T. Voici à nouveau l'expression de la TFD pour N échantillons :

$$U_{N,n}^{tfd}=\sum _{k=0}^{N-1}{u}_{k}exp(-j2\pi k\frac{n}{N})\text{(29)}$$

Complétons les N nombres $u_0,u_1,\dots u_{N-1}$ par N nombres $u_N,u_{N+1}\dots 2N-1$ tous nuls. La TFD de ces 2N nombres s'écrit :

$$U_{2N,n}^{tfd}=\sum _{k=0}^{2N-1}{u}_{k}exp(-j2\pi k\frac{n}{2N})=\sum _{k=0}^{N-1}{u}_{k}exp(-j2\pi k\frac{n}{2N})\text{(30)}$$

On obtient la TFD des N nombres $u_0,u_1,\dots u_{N-1}$ mais pour les 2N fréquences $f_n=\frac{T}{2N}$ avec $n=0,1,\dots 2N-1$. En conséquence, ajouter N zéros au signal échantillonné puis calculer la TFD des 2N valeurs revient à calculer la TFTD de $u_0,u_1,\dots u_{N-1}$ mais pour un échantillonnage de fréquence comportant 2 fois plus de fréquences dans l'intervalle [0,f_e]. La généralisation de cette méthode est immédiate : on peut doubler la résolution en fréquence de la TFTD en multipliant par 2 le nombre d'échantillons, les échantillons ajoutés étant tous nuls. Cette méthode de calcul de la TFTD est la méthode du remplissage par des zéros.

La fonction numpy.fft.fft calcule la TFD. Comme son nom l'indique, cette fonction fait le calcul au moyen de l'algorithme de Transformée de Fourier rapide (Fast Fourier Transform). Cet algorithme nécessite cependant que N soit une puissance de 2. En réalité, numpy.fft.fft fonctionne pour une taille de tableau quelconque. L'algorithme implémenté par cette fonction consiste à appliquer la FFT tant que la taille du tableau est divisible par 2 puis à terminer par un calcul direct de la TFD. Par exemple, si le tableau contient 100 éléments, on se ramène au calcul de deux TFD à 50 éléments, puis à 4 TFD à 25 éléments. Puisque 25 est impair, le calcul d'une TFD à 25 éléments doit se faire de manière directe, c'est-à-dire par calcul direct des 25 sommes définissant la TFD (chacune comportant 25 termes). Pour que le calcul soit optimisé, il est donc préférable que N soit une puissance de 2.

Remarque : dans la présentation théorique ci-dessus, l'intervalle de temps de la fenêtre est [-T/2,T/2]. En pratique, on utilise plutôt [0,T]. Les tableaux sont indexés à partir de zéro et on pose $t_k=k\frac{T}{N}$ avec $k=0,1,\dots N-1$.

Considérons comme exemple le cas de la fonction périodique $s\left(t\right)=S_{0}cos(2\pi f_1+\varphi )$ dont nous avons étudié la TF précédemment pour différentes fonctions de fenêtrage. Le but est d'utiliser la TFD pour obtenir la TF. Plus précisément, on cherche à calculer $U\left(f\right)/T$ (TF normalisée).

On reprend les conventions du calcul précédent : la fréquence du signal est fixée à f₁=1 (ce qui revient à utiliser sa période comme unité de temps). La durée T de la fenêtre doit être assez grande car 1/T est la résolution fréquentielle du spectre obtenu par la TFD. Le nombre d'échantillons N pour cette durée doit être assez grand pour satisfaire la condition de Nyquist-Shannon. Dans le cas présent, il faut que N>2T. Par ailleurs, on évitera de choisir une valeur de T entière car la durée du signal n'a aucune raison d'être multiple de sa période.

from numpy.fft import fft,fftfreq
f1 = 1
T = 100.48
N = 500
Te = T/N
fe = 1/Te
t = np.arange(N)*Te 
s = np.cos(2*np.pi*f1*t+0.274)
u = s # fenêtre rectangulaire
tfd = fft(u)/N
freq = np.arange(N)*1/T
figure(figsize=(16,6))
plot(freq,np.absolute(tfd),'r.')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
grid()
  
                  

fig7.pdf

L'intervalle de fréquence représenté est [0,f_e[. Si l'on s'intéresse à la TFTD, c'est-à-dire au spectre du signal échantillonné, cette représentation est pertinente mais il faut se rappeler que la TFTD est périodique, de période f_e. Si l'on veut comparer à la TF du signal fenêtré, il est préférable de ramener la motié droite du spectre aux fréquences négatives. La fonction numpy.fft.fftfreq calcule les valeurs des fréquences avec des fréquences négatives pour la seconde partie du tableau :

figure(figsize=(16,6))
plot(fftfreq(N,d=Te),np.absolute(tfd),'r.')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
grid()
                  

fig8.pdf

Nous pouvons faire la comparaison entre cette TFD et la TF normalisée (divisée par T ) du signal pour la fenêtre rectangulaire (celle-ci doit être échantillonnée avec un très grande nombre de points). La TF est tracée sous la forme d'une ligne alors qu'on garde des points pour la TFD.

f = np.linspace(-fe/2,fe/2,100000)
U = 0.5*(np.sinc((f-f1)*T)+np.sinc((f+f1)*T)) # TF normalisée
figure(figsize=(16,6))
plot(f,np.absolute(U),label='TF')
plot(fftfreq(N,d=Te),np.absolute(tfd),'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
grid()
legend(loc='upper right',fontsize=16)
                  

fig9.pdf

Voici un détail :

figure(figsize=(16,6))
plot(f,np.absolute(U),label='TF')
plot(fftfreq(N,d=Te),np.absolute(tfd),'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
xlim(0.9,1.1)
grid()
legend(loc='upper right',fontsize=16)
                  

fig10.pdf

La TFD est une approximation très grossière de la TF dans le sens où l'échantillonnage de la fréquence n'est pas assez fin pour reproduire toutes les variations de la TF (avec une simple interpolation linéaire). L'inconvénient majeur est que la hauteur de la raie ne peut être obtenue (la hauteur théorique est 0,5). La raison de ce problème réside dans le fait que la TF de la fenêtre rectangulaire, de largeur T, ne peut être obtenue avec assez de détails puisque les échantillons sont répartis seulement sur la largeur de la fenêtre. La TFD contient pourtant toutes les informations permettant de reconstruire la TFTD. Pour y parvenir, il faut étendre la largeur de l'intervalle de temps échantillonné (Fenêtres pour l'analyse spectrale).

La TF de u(t) est le produit de convolution de la TF de s(t) par la TF de w(t). Il s'en suit qu'il faut étendre la largeur de l'intervalle de temps échantillonné pour obtenir une représentation convenable de la TF du signal fenêtré (c.a.d. un échantillonnage de la fréquence assez fin). Notons T' la durée de l'intervalle de temps étendu et posons T'=N'T_e. Sachant que w(t) est nulle en dehors de l'intervalle [0,T], cet extension de l'intervalle de temps consiste simplement à ajouter N'-N zéros. Si on s'intéresse seulement au module de la TFD, ces zéros peuvent être indifféremment ajoutés de part et d'autre du signal échantillonné (à gauche et à droite) ou simplement à droite, cette dernière solution étant plus simple à implémenter. Notons n_z le nombre de blocs de zéros que l'on veut ajouter (un bloc contient N zéros). Il est judicieux (mais pas obligatoire) de choisir N' égal à une puissance de 2 afin d'optimiser le calcul de la TFD par FFT. Cette méthode est appelée remplissage par des zéros (zero padding). Il faut remarquer qu'elle ne modifie ni la fréquence d'échantillonnage, ni la durée du signal analysée (qui reste T). Cette méthode revient en fait à augmenter la durée d'analyse de la fenêtre, qui est insuffisante si elle est égale à la durée T de la fenêtre. On peut aussi la voir comme une méthode de reconstruction de la TFTD à partir de la TFD : comme il a été démontré plus haut, ajouter à un signal de N échantillons N zéros avant de calculer la TFD des 2N valeurs revient à multiplier par deux le nombre de fréquences dans l'intervalle [-f_e/2,f₂/2]. Ainsi, si on ajoute 3 blocs de zéros, on multipie par 4 la résolution fréquentielle de la TFTD, si on ajoute 7 blocs, on la multiplie par 8.

nz = 5
p = int(np.log(nz*N)/np.log(2))+1
Np = 2**p # N' 
u_etendu=np.concatenate((u,np.zeros(Np-N)))
tfd = fft(u_etendu)/N
figure(figsize=(16,6))
plot(f,np.absolute(U),label='TF')
plot(fftfreq(Np,d=Te),np.absolute(tfd),'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
grid()
legend(loc='upper right',fontsize=16)
                      

fig11.pdf

et voici la raie en détail :

figure(figsize=(16,6))
plot(f,np.absolute(U),label='TF')
plot(fftfreq(Np,d=Te),np.absolute(tfd),'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
xlim(0.9,1.1)
grid()
legend(loc='upper right',fontsize=16)
                      

fig12.pdf

Le remplissage par 5 blocs de zéros donne donc une représentation satisfaisante de la TF (ou de la TFTD). Dans le cas de la fenêtre rectangulaire, les lobes secondaires sont très intenses. Le point important est que ce spectre permet une estimation correcte de la hauteur de la raie (ici 0,5). Il permet aussi d'améliorer la précision de la détermination de la fréquence correspondant au maximum, qui est la fréquence du signal s(t). Après remplissage par des zéros, la fréquence d'une raie peut être déterminée avec une précision bien meilleure que 1/T. Ce résultat peut sembler étonnant mais il faut noter que la résolution fréquentielle reste 1/T (et non pas 1/T') : deux raies dont la différence des fréquences est inférieure à 1/T ne seront pas résolues et l'ajout de zéro ne change rien à cette limite. Pour obtenir la fréquence du maximum, la méthode la plus simple consiste à prendre le point dont le module est le plus grand :

f1_exp = np.argmax(np.absolute(tfd[0:Np//2]))*1/(Te*Np)               
                           

print(f1_exp)
--> 0.9998394425507564

Voici les résultats pour 10 blocs de zéros ajoutés :

nz = 10
p = int(np.log(nz*N)/np.log(2))+1
Np = 2**p # N'
u_etendu=np.concatenate((u,np.zeros(Np-N)))
tfd = fft(u_etendu)/N
figure(figsize=(16,6))
plot(f,np.absolute(U),label='TF')
plot(fftfreq(Np,d=Te),np.absolute(tfd),'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
xlim(0.9,1.1)
grid()
legend(loc='upper right',fontsize=16)
                      

fig13.pdf

f1_exp = np.argmax(np.absolute(tfd[0:Np//2]))*1/(Te*Np)               
                           

print(f1_exp)
--> 0.9998394425507564

L'inconvénient de la fenêtre rectangulaire est la hauteur des lobes secondaires de sa TF. Voyons le calcul avec la fenêtre de Blackman, dont les lobes secondaires sont très faibles (voir Fenêtres pour l'analyse spectrale).

On commence par le calcul de la TFD avec N échantillons sur l'intervalle [0,T] (la largeur de la fenêtre) :

from scipy.signal.windows import blackman
                         
                             

u = s*blackman(N)
tfd = fft(u)/N
freq = np.arange(N)*1/T
f = np.linspace(-fe/2,fe/2,100000)
U = 0.5*(Wb(f1,T,f) + Wb(-f1,T,f)) # TF normalisée
figure(figsize=(16,6))
plot(f,np.absolute(U)/0.42,label='TF')
plot(fftfreq(N,d=Te),np.absolute(tfd)/0.42,'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
grid()
legend(loc='upper right',fontsize=16)
                  

fig14.pdf

Nous pouvons faire le même constat que pour la fenêtre rectangulaire : l'échantillonnage de la fenêtre sur sa durée T est insuffisant pour donner une approximation visuellement correcte de la TF. Il faut compléter par des zéros de manière à calculer la TFD pour un intervalle T'=N'T_e nettement plus grande que T :

nz = 5
p = int(np.log(nz*N)/np.log(2))+1
Np = 2**p # N'
u_etendu=np.concatenate((u,np.zeros(Np-N)))
tfd = fft(u_etendu)/N
freq = np.arange(N)*1/T
f = np.linspace(-fe/2,fe/2,100000)
U = 0.5*(Wb(f1,T,f) + Wb(-f1,T,f)) # TF normalisée
figure(figsize=(16,6))
plot(f,np.absolute(U)/0.42,label='TF')
plot(fftfreq(Np,d=Te),np.absolute(tfd)/0.42,'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
grid()
legend(loc='upper right',fontsize=16)
                  

fig15.pdf

figure(figsize=(16,6))
plot(f,np.absolute(U)/0.42,label='TF')
plot(fftfreq(Np,d=Te),np.absolute(tfd)/0.42,'r.',label='TFD')
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U^{tfd}|$",fontsize=16)
ylim(0,1)
xlim(0.9,1.1)
grid()
legend(loc='upper right',fontsize=16)
                  

fig16.pdf

5. Fonction de calcul du spectre

Nous résumons la méthode de calcul qui permet d'obtenir la TFTD d'un signal échantillonné sur une durée T. Dans la mesure où la condition de Nyquist-Shannon est vérifiée, cette TFTD est égale à la TF.

• Obtenir N échantillons du signal s(t) sur une durée T.
• Multiplier ces échantillons par une fenêtre de largeur T.
• Ajouter n_z bloc de N zéros (environ) de manière à obtenir N' échantillons (égal à une puissance de 2).
• Calculer la TFD de ces N' échantillons par l'algorithme de transformée de Fourier rapide.
• Normaliser la TFD en divisant par N si nécessaire et par la valeur maximale de la TF de la fenêtre (1 pour rectangulaire, 0,42 pour Blackman, 0,54 pour Hamming et 0,5 pour Hann).
• Calculer les fréquences sachant qu'elles sont espacées de $\frac{1}{T_{e}N\text{'}}$.
• Tracer le module de la TFD normalisée en fonction de la fréquence : elle correspond en théorie au module de la TF normalisée du signal fenêtré.

Remarque : la normalisation de la TF consiste à la diviser par T et aussi par la valeur maximale de la TF de la fenêtre. Cette normalisation permet de lire directement sur le spectre les amplitudes des composantes sinusoïdales du signal.

La fonction suivante effectue le calcul du spectre d'un signal échantillonné. Ses arguments sont :

• s : les échantillons du signal.
• Te : la période d'échantillonnage.
• win : le nom de la fenêtre (rect, hamming, hann ou blackman).
• nz : le nombre de blocs de zéro à ajouter.

Elle renvoie les fréquences de 0 à f_e, les fréquences de -f_e/2 à f_e/2 et la TFD normalisée.

from scipy.signal.windows import get_window

              
def spectre(s,Te,win,nz):
    N = len(s)
    if win=='rect':
        w = np.ones(N)
        a = 1
    elif win=='hamming': 
        w = get_window('hamming',N)
        a = 0.54
    elif win=='hann':
        w = get_window('hann',N)
        a = 0.5
    elif win=='blackman':
        w = get_window('blackman',N)
        a = 0.42
    else:
        return None
    p = int(np.log(nz*N)/np.log(2))+1
    Np = 2**p
    u = np.concatenate((s*w,np.zeros(Np-N)))
    tfd = fft(u)/N/a
    freq = np.arange(Np)*1/(Np*Te)
    return freq,fftfreq(Np,d=Te),tfd
              

La fonction numpy.fft.fft comporte un argument optionnel n qui permet de préciser le nombre d'échantillons utilisés pour le calcul de la TFD. Si ce nombre est supérieur à la taille du signal transmis, celui-ci est complété par des zéros pour parvenir à la taille demandée. On peut donc écrire :

tfd = fft(s*w,Np)/N/a

Voici par exemple, l'analyse d'un signal comportant deux sinusoïdes :

f1 = 1
f2 = 1.5
T = 103.58
N = 1000
Te = T/N
t = np.arange(N)*Te
s = 1.0*np.cos(2*np.pi*f1*t)+ 0.3*np.sin(2*np.pi*f2*t)
freq,freq1,tfd = spectre(s,Te,'blackman',10)
figure(figsize=(16,6))
plot(freq,np.absolute(tfd))
xlabel(r"$f/f_1$",fontsize=16)
ylabel(r"$|U|/T$",fontsize=16)
grid()
              

fig17.pdf

6. Analyse des signaux périodiques

7. Application : analyse d'un signal audio

L'analyse spectrale d'un signal audio a pour objectif de simuler la perception des sons. Nous sommes capable d'identifier un son, même très bref, comme étant approximativement périodique et de percevoir sa hauteur , c'est-à-dire (en première approximation) la fréquence de son harmonique fondamental. La hauteur d'un son musical ou vocal d'une durée de 1/10 seconde est perçue sans difficulté. Pour une fréquence de 400 Hz, une durée aussi courte correspond à 40 cycles. Pour simuler la perception des sons, il faut faire une analyse spectrale avec une fenêtre dont la durée est de l'ordre de grandeur de la résolution temporelle de l'audition (20 ms pour des sons de fréquence moyenne). Cette fenêtre doit être déplacée sur l'axe du temps (fenêtre glissante), par sauts plus ou moins grand. En général, deux positions successives de la fenêtre se recouvrent partiellement. Dans ce paragraphe, nous abordons le problème du calcul du spectre par fenêtre glissante et de sa représentation graphique sous forme de spectrogramme.

L'analyse spectrale avec une fenêtre ne comportant que 40 cycles ne permet pas de distinguer deux sons périodiques dont l'écart de fréquence est inférieur à un 1/40 ième de la fréquence, soit 10 Hz pour une fréquence de 400 Hz. Cependant, cela ne signifie pas que la hauteur d'un son analysé sur une durée de 1/10 secondes ne puisse pas être déterminée avec une plus grande précision. En effet, la méthode du remplissage par des zéros développée précédemment montre que la précision de la détermination de la fréquence d'un signal sinusoïdal peut être bien plus petite que l'inverse de la durée de la fenêtre d'analyse. Autrement dit, il faut bien distinguer la résolution fréquentielle qui permet de distinguer deux sons de fréquences très voisines de la résolution fréquentielle qui permet de déterminer la hauteur d'un son.

La fréquence d'un La₁ est 440 Hz (par convention pour un accord à 440 Hz). La fréquence du La₁# est alors de 466 Hz environ, soit un écart de 26 Hz. On voit donc que les différents sons quasi périodiques qui sont perçus simultanément ont rarement une différence de fréquence inférieure à 10 Hz (pour 400 Hz). Cela peut néanmoins se produire, par exemple le La₁ émis par un piano est en fait émis par trois cordes identiques qui en principe émettent exactement à la même fréquence lorsque le piano est parfaitement accordé. Lorsqu'un léger désaccord apparaît entre les trois cordes, nous ne percevons pas distinctement les trois fréquences mais le son a un timbre différent. Il est ainsi possible de constater qu'un piano est désaccordé simplement en jouant une seul note. Si la note est maintenue pendant plusieurs secondes, on entend par ailleurs des battements de l'amplitude (qui sont utilisés par l'accordeur pour accorder les trois cordes).

Comme exemple d'une analyse par fenêtre glissante, on considère un enregistrement (dans un ficher WAVE) de 30 s comportant un enchaînement mélodique (environ 6 notes par seconde).

Voici la lecture du fichier :

import scipy.io.wavfile as wave
fe,data = wave.read("ZOOM0004_TrLR.WAV")
Te = 1/fe
s1 = data[:,0]
s2 = data[:,1]
s1 = s1/s1.max()
s2 = s2/s2.max()
N = len(s1)
t = np.arange(N)*Te          
                  

print(fe)
--> 44100

print(N)
--> 1347456

La fréquence d'échantillonnage est 44,1 kHz. Voici la représentation temporelle du signal complet (voie 1 de l'enregistrement stéréo) :

figure(figsize=(16,6))
plot(t,s1)
grid()
xlabel('t (s)',fontsize=16)
                  

fig22.pdf

et un détail d'une seconde :

figure(figsize=(16,6))
plot(t,s1)
xlim(10,11)
grid()
xlabel('t (s)',fontsize=16)
                  

fig23.pdf

L'enchaînement des notes se devine mais il est impossible de les séparer par un algorithme simple. Pour faire l'analyse spectrale, il faut choisir la durée T de la fenêtre glissante. Dans le cas présent, les sons ont tous la même durée donc il semble à première vue pertinent de choisir T à peu près égale à cette durée. On choisit donc une fenêtre de longueur 8192, soit une durée d'environ 0,186 s. Il faut aussi choisir le déplacement de la fenêtre entre deux positions consécutives. On peut par exemple choisir un déplacement égal à 1/4 de la longueur de la fenêtre. Voici par exemple les spectres obtenus pour 8 positions consécutives de la fenêtre démarrant à l'instant t=10 s. On limite le tracé du spectre de 0 à 3000 Hz. Les spectres sont décalés verticalement pour plus de visibilité.

size = 2**13
step = size//4
ti = 10.0
i = int(ti/Te)
win = 'hann'
nz = 1
figure(figsize=(12,7))
for k in range(8):
    freq,freq1,U = spectre(s1[i:i+size],Te,win,nz)
    i += step
    plot(freq,k+np.absolute(U)*10,label='%d'%k)
xlim(0,3000)
ylim(0,10)
grid()
legend(loc='upper right')
xlabel('f (Hz)',fontsize=16)
                   

fig24.pdf

Évidemment, il y a peu de chance que la position d'une fenêtre coïncide avec ce que nous percevons comme un son de hauteur donnée (ici une note jouée), ni même en général que la durée de la fenêtre corresponde exactement à la durée d'un son.

L'enchaînement des spectres peut être obtenu par une animation :

from matplotlib.animation import FuncAnimation
size = 2**13
step = size//4
t1,t2 = 5,15 # instants initial et final
i1,i2 = int(t1/Te),int(t2/Te)
nwin = int((i2-i1)/step)
freq,freq1,U = spectre(s1[i1:i1+size],Te,'blackman',nz=5)
fig,ax = plt.subplots()
line0, = ax.plot(freq,np.absolute(U))
ax.grid()
ax.set_xlabel("f (Hz)")
ax.set_ylabel("U")
ax.axis([0,3000,0,0.5])

def animate(j):
    global s1,size
    i = i1+j*step
    freq,freq1,U = spectre(s1[i:i+size],Te,'hann',nz=5)
    line0.set_xdata(freq)
    line0.set_ydata(np.absolute(U))

interval = 100
n = nwin
anim = FuncAnimation(fig,animate,n,interval=interval,repeat=True)
anim.save('analyseSon.mp4')
plt.show()               
                    

La perception de l'amplitude des sons est plutôt logarithmique, du moins lorsque l'intensité sonore est relativement grande. On trace donc l'amplitude en échelle décibel :

Un spectrogramme est une représentation graphique des spectres obtenus par fenêtre glissante sous la forme d'une image. La fonction matplotlib.pyplot.specgram permet d'obtenir un spectrogramme. Voici comment on l'utilise pour obtenir un spectrogramme d'une partie de l'enregistrement :

NFFT = 8192 # largeur de la fenêtre
figure(figsize=(16,6))
t1,t2 = 5,15 # instants initial et final
i1,i2 = int(t1/Te),int(t2/Te)
spectrum,fresq,t,im = specgram(s1[i1:i2],NFFT=NFFT,noverlap=NFFT*3//4,Fs=fe,pad_to=NFFT*4,mode='magnitude',scale='dB',window=lambda data: data*np.blackman(len(data)))
ylim(0,3000) 
ylabel("f (Hz)",fontsize=16)
xlabel("t (s)",fontsize=16)
                      

fig25.pdf

Sur ce spectrogramme, l'amplitude est calculée en décibel afin de faciliter le répérage des harmoniques. noverlap est le recouvrement de deux fenêtres consécutives. Pour faire un déplacement de 1/4 de fenêtre il faut un recouvrement de 3/4 de fenêtre. pad_to est le nombre total de points de la FFT (plus grand que NFFT), les points ajoutés étant nuls (remplissage par des zéros). La fenêtre peut être donnée soit sous la forme d'une fonction soit sous la forme d'un tableau. La documentation de specgram ne donne pas de détail sur l'algoritme de génération de l'image mais on peut le deviner : l'axe du temps est divisé en un nombre de pixels suffisants pour que l'intervalle de temps correspondant au recouvrement soit bien visible. Pour chaque rangée verticale de pixels, le spectre affiché est un mélange des spectres obtenus par les fenêtres qui se superposent à l'instant correspondant (si le recouvrement est 3/4, il y en a 4). Pour chaque fréquence, l'amplitude est convertie en couleur au moyen d'une palette de couleurs.

En principe, le son correspondant à une note est caractérisé par un fondamental et des harmoniques dont les fréquences sont multiples de la fréquence fondamentale. Il faut donc repérer dans le spectrogramme ce type de configuration. L'inconvénient d'avoir choisi une durée de fenêtre égale à peu près à la durée d'une note est qu'une position aléatoire de la fenêtre a toutes les chances de contenir un mélange de sons de plusieurs hauteurs. Pour isoler ces sons, il est donc préférable de diminuer la durée de la fenêtre et d'annuler le recouvrement :

NFFT = 1024 # largeur de la fenêtre
figure(figsize=(16,6))
t1,t2 = 5,15 # instants initial et final
i1,i2 = int(t1/Te),int(t2/Te)
spectrum,fresq,t,im = specgram(s1[i1:i2],NFFT=NFFT,noverlap=0,Fs=fe,pad_to=NFFT*4,mode='magnitude',scale='dB',window=lambda data: data*np.blackman(len(data)))
ylim(0,3000) 
ylabel("f (Hz)",fontsize=16)
xlabel("t (s)",fontsize=16)
                      

fig26.pdf

L'inconvénient de la réduction de la durée de la fenêtre glissante est bien visible sur le spectrogramme : la résolution fréquentielle est moins bonne. Dans une analyse temps-fréquence, il faut trouver un compromis entre la résolution temporelle et la résolution fréquentielle. Ici, la largeur de la fenêtre est de 23 ms, ce qui donne une résolution fréquentielle de 43 Hz. Cette résolution fréquentielle est cependant largement suffisante pour distinguer les harmoniques. En observant attentivement le spectrogramme, on identifie bien pour la plupart des positions de la fenêtre un fondamental (jaune très lumineux) est des harmoniques.