Échantillonnage et reconstruction d'un signal

2.a. Transformée de Fourier et spectre d'un signal

La transformée de Fourier du signal u(t) est définie par :

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

Pour que cette intégrale soit définie (quelle que soit la fréquence f), la fonction u(t) doit bien sûr vérifier certaines conditions. Si le signal u(t) est limité dans le temps (fonction à support borné) les bornes de l'intégrale définissant la transformée de Fourier sont finies et la fonction U(f) est donc définie pour toute fréquence f réelle. La transformée de Fourier pose des difficultés lorsque u(t) est périodique, c'est-à-dire lorsqu'il existe un réel P tel que u(t+P)=u(t) pour tout t. Une fonction périodique n'est évidemment pas à support borné. Sous certaines conditions de régularité (classe $C^{\infty }$ par morceaux), elle peut s'écrire comme la somme de sa série de Fourier. Chaque terme de cette série est une fonction périodique sinusoïdale. Considérons donc un signal sinusoïdal de fréquence f₁ :

$$u\left(t\right)=A_{1}cos\left(2\pi f_{1}t\right)\text{(2)}$$

et l'intégrale suivante :

$$\begin{array}{ccc}\hfill {\int }_{-\frac{T}{2}}^{\frac{T}{2}}u\left(t\right)exp(-j2\pi ft)dt& =\frac{A_1}{2}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}(exp(-j2\pi (f_1-f)t)+exp(-j2\pi (-f_1-f)t))dt\hfill & \hfill \text{(3)}\\ \hfill & =\frac{A_1}{2}T(sinc\left((f_1-f)T\right)+sinc\left((f_1+f)T\right))\hfill & \hfill \text{(4)}\end{array}$$

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

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

Lorsque $T\to \infty$, l'intégrale précédente tend vers zéro si $f\ne f_1$ et $f\ne -f_1$. Si $f=f_1$, $sinc\left((f-f_1)T\right)=1$ et l'intégrale tend donc vers l'infini. Il en est de même si $f=-f_1$. La transformée de Fourier d'une fonction sinusoïdale (et d'une fonction périodique) n'est donc pas définie en tant que fonction. Dans ce cas, U(f) désigne une distribution. Remarquons que les signaux qu'on considère en pratique (c'est-à-dire dans la réalité) sont toujours à support borné. Par conséquent, leur transformée de Fourier est une fonction.

Le signal u(t) peut être reconstitué à partir de sa transformée de Fourier, par la transformée de Fourier inverse :

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

Cette expression montre que U(t) s'écrit comme une somme continue de fonctions sinusoïdales et que U(f) est l'amplitude complexe de la composante de fréquence f dans cette somme. En conséquence, la transformée de Fourier du signal définit son spectre en fréquences, ou plus simplement son spectre. On dit aussi que U(f) est la représentation fréquentielle du signal. Si u est à valeurs réelles, on a $U\left(f\right)=U^{*}\left(f\right)$ (complexe conjugué).

La représentation graphique (partielle) du spectre consiste à tracer le module de U(f) en fonction de la fréquence.

Un signal est dit à bande de fréquence limitée si sa transformée de Fourier est à support borné, c'est-à-dire s'il existe une fréquence maximale f_m telle que :

$$U\left(f\right)=0si\left|f\right|>f_m\text{(7)}$$

Certains signaux théoriques ne sont pas à bande de fréquence limitée. Par exemple la fenêtre rectangulaire de largeur T, 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}$$

a pour transformée de Fourier :

$$W_r\left(f\right)={\int }_{-\frac{T}{2}}^{\frac{T}{2}}exp(-j2\pi ft)dt=Tsinc\left(fT\right)\text{(8)}$$

Cette fonction n'est pas à support borné mais elle tend vers zéro lorsque $f\to \infty$, donc son module peut être considéré comme négligeable à partir d'une certaine fréquence. D'un point de vue pratique, tout signal peut être considéré comme limité en fréquence mais la question est de savoir comment la fréquence maximale f_m se situe par rapport à la fréquence d'échantillonnage.

2.b. Échantillonnage

Soit u(t) une fonction du temps (ou de tout autre variable réelle) à valeurs complexes. Cette fonction représente un signal analogique, c'est-à-dire un signal à temps continu. Une grandeur physique est représentée par un nombre réel mais tout ce qui est exposé ci-après s'applique aussi bien si u(t) est complexe.

L'échantillonnage de u(t) consiste à considérer les valeurs qu'elle prend pour des instants discrets régulièrement espacés :

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

où q est un nombre entier. T_e est la période d'échantillonnage et on utilisera aussi la fréquence d'échantillonnage f_e=1/T_e. La suite de nombres complexes u_q est le signal échantillonné. Il s'agit d'un signal à temps discret. On emploie aussi la dénomination signal numérique mais dans ce cas on fait aussi référence à la représentation de u_q (nombre réel) sous la forme d'un nombre entier. Dans ce document, on ne s'intéresse pas au problème de la représentation d'un nombre réel par un nombre entier (problème de la quantification) mais il s'agit évidemment d'un problème de grande importance pratique.

2.c. Transformée de Fourier à temps discret

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

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

Multipliée par $T_e$, elle constitue une approximation de la transformée de Fourier, d'autant meilleure que la période d'échantillonnage est petite. En effet, on a :

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

La TFTD définit le spectre du signal échantillonné, comme la TF définit celui du signal analogique.

La TFTD est périodique :

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

Le spectre d'un signal échantillonné est donc périodique, de plus petite période égale à la fréquence d'échantillonnage.

Remarquons que la fréquence intervenant dans la TFTD est une variable réelle continue (comme pour la TF).

Le calcul numérique de la TFTD n'est possible que si la somme comporte un nombre fini de termes. C'est le cas si u(t) est à support borné. Notons [0,T] un intervalle en dehors duquel u(t)=0 et échantillonnons cet intervalle par N instants :

$$t_n=n\frac{T}{N}pourn=0,1,\dots N-1\text{(14)}$$

avec une période d'échantillonnage $T_e=\frac{T}{N}$.

Cela nous amène à la définition de la transformée de Fourier discrète (TFD) :

$$U_k=\sum _{n=0}^{N-1}{u}_{n}exp(-j2\pi \frac{kn}{N})pourk=0,1\dots N-1\text{(15)}$$

Soient les fréquences :

$$f_k=k\frac{1}{T}pourk=0,1,\dots N-1\text{(16)}$$

On a :

$$U_k=U_{T_e}\left(f_k\right)\text{(17)}$$

Autrement dit, la TFD permet de calculer la TFTD mais seulement pour les fréquences ci-dessus et si le signal est nul en dehors de l'intervalle [0,T]. On peut donc dire que la TFD réalise un échantillonnage en fréquence de la TFTD. Plus T est grand, plus cet échantillonnage est fin et plus on s'approche de la TFTD. En pratique, on peut augmenter la durée T en ajoutant des zéros au signal (méthode de remplissage par des zéros).

Dans le cas (théorique) où le signal u(t) n'est pas à support borné (par exemple un signal périodique), la TFD permet de calculer le spectre d'un signal échantillonné égal à u_q sur l'intervalle [0,T] et nul en dehors de cet intervalle. Autrement dit, le calcul d'une TFD suppose implicitement que le signal est prolongé en dehors de l'intervalle [0,T] par des zéros. Étrangement, les N nombres U_k sont aussi les coefficients de Fourier d'une fonction de période T, donc d'une fonction prolongée par périodicité et non pas par des zéros (cette propriété est vraie si la condition de Nyquist-Shannon, énoncée plus loin, est satisfaite).

2.d. Théorème du repliement

Le théorème du repliement (aliasing theorem [1]) donne la relation entre le spectre du signal échantillonné (sa TFTD) et le spectre du signal analogique (sa TF) :

$$U_{T_e}\left(f\right)=\frac{1}{T_e}\sum _{p=-\infty }^{\infty }U(f+\frac{p}{T_e})\text{(18)}$$

Lorsque $T_e\to \infty$, on a $\mathrm{(12)}$ c'est-à-dire que la TF se confond avec la TFTD (au facteur T_e près).

Pour visualiser la relation entre la TFTD et la TF, supposons que le signal soit à bande de fréquence limitée (fréquence maximale f_m) et que $f_e>2f_m$.

La figure suivante montre un exemple schématique de spectres (tracé des parties réelles) :

Figure pleine page

On définit aussi la fréquence de Nyquist :

$$f_n=\frac{f_e}{2}\text{(19)}$$

Dans le cas où $f_e>2f_m$, la TFTD (multipliée par T_e) est simplement un prolongement de la TF par périodicité car les différents termes de la somme $\mathrm{(18)}$ ne se recouvrent pas. Si on s'intéresse à la TFTD (multipliée par T_e) dans l'intervalle [-f_n,f_n], elle correspond exactement à la TF. Rappelons que la TFD permet de calculer la TFTD pour certaines fréquences (multiples de 1/T) mais seulement si le signal est nul en dehors de l'intervalle [0,T].

Si $f_e\le 2f_m$, la somme $\mathrm{(18)}$ se fait sur des fonctions qui se recouvrent :

Figure pleine page

Dans ce cas, la TFTD dans l'intervalle [-f_n,f_n] (multipliée par T_e) n'est pas égale à la TF à cause des zones de recouvrement des différents termes de la somme. Ce phénomène est nommé repliement de bande car la partie de la TF dans l'intervalle [f_e/2,f_e] qui se retrouve dans l'intervalle [0,f_e] peur se voir comme un repliement du spectre de l'intervalle [0,f_e/2], c'est-à-dire une symétrie par rapport à la fréquence f_e/2. Le terme anglais pour désigner ce phénomène est aliasing, qui signifie que le spectre de l'intervalle [f_e/2,f_e] se réplique (alias) dans l'intervalle [0,f_e/2]. Il ne faut toutefois pas perdre de vue que le repliement est en fait une somme (relation $\mathrm{(18)}$) et que le nombre de termes sommés peut être supérieur à 2 (auquel cas la notion de repliement perd son sens).

2.e. Théorème de l'échantillonnage

La question importante qui se pose est la suivante : peut-on reconstuire le signal u(t) (fonction du temps) à partir des échantillons u_q ?

Cette question en contient une autre : comment peut-on reconstruire une fonction du temps à partir des échantillons ? La transformée de Fourier inverse permet de réaliser cette reconstruction. Pour obtenir un signal (fonction du temps), il faut appliquer la transformée de Fourier inverse à la fonction de la fréquence égale à $T_e{U}_{T_e}\left(f\right)$ sur l'intervalle [-f_n,f_n], égale à zéro en dehors de cet intervalle. Le signal recontruit s'écrit donc :

$$\tilde{u}\left(t\right)=T_e{\int }_{-f_n}^{f_n}{U}_{T_e}\left(f\right)exp\left(j2\pi ft\right)df\text{(20)}$$

Si $f_e\le 2f_m$, la fonction $U_{T_e}\left(f\right)$ n'est pas égale à $U\left(f\right)$ sur tout l'intervalle [-f_n,f_n]. En conséquence $\tilde{u}\left(t\right)$ n'est pas identique à $u\left(t\right)$ et il est donc impossible de recontruire le signal à partir des échantillons.

Si $f_e>2f_m$, on a $U_{T_e}\left(f\right)=U\left(f\right)$ sur l'intervalle [-f_n,f_n] et l'expression $\mathrm{(20)}$ est identique à celle de la transformée de Fourier inverse de U(f), ce qui implique que $\tilde{u}\left(t\right)=u\left(t\right)$ pour tout t. En conséquence, le signal recontruit est identique au signal analogique u(t).

Explicitons $\mathrm{(20)}$ au moyen de l'expression de la TFTD $\mathrm{(11)}$ :

$$\begin{array}{ccc}\hfill \tilde{u}\left(t\right)& =T_e{\int }_{-f_n}^{f_n}\sum _{q=-\infty }^{\infty }{u}_{q}exp(-j2\pi fq{T}_e)exp\left(j2\pi ft\right)df\hfill & \hfill \text{(21)}\\ \hfill & =\sum _{q=-\infty }^{\infty }{u}_q{T}_e{\int }_{-f_n}^{f_n}exp\left(j2\pi f(t-q{T}_e)\right)df\hfill & \hfill \text{(22)}\\ \hfill & =\sum _{q=-\infty }^{\infty }{u}_{q}sinc(\frac{t}{T_e}-q)\hfill & \hfill \text{(23)}\end{array}$$

La fonction sinus cardinale est définie par $\mathrm{(5)}$.

Le théorème de l'échantillonnage de Nyquist-Shannon s'en déduit : pour que le signal u(t) puisse être reconstruit à partir du signal échantillonné, il faut et il suffit que le signal soit à bande de fréquence limitée (existence d'une fréquence maximale f_m dans sa TF) et que la fréquence d'échantillonnage soit strictement supérieure au double de f_m. La reconstruction se fait de la manière suivante :

$$u\left(t\right)=\sum _{q=-\infty }^{\infty }{u}_{q}sinc(\frac{t}{T_e}-q)\text{(24)}$$

Cette relation est nommée formule de Shannon. Dans le cas où $f_e\le 2f_m$ (sous-échantillonnage), il est bien sûr possible d'appliquer cette formule mais elle ne donne pas u(t).

Pour un indice q donné, la fonction sinc qui apparaît dans la somme est un sinus cardinal dont le maximum est à l'instant t=qT_e, qui est l'instant correspondant à l'échantillon u_q. Par ailleurs, ce sinus cardinal s'annule pour tous les autres instants du signal échantillonné. La formule peut s'interpréter comme un produit de convolution. Posons :

$$h_q\left(t\right)=sinc(\frac{t}{T_e}-q)\text{(25)}$$

On peut alors écrire :

$$u\left(t\right)=\sum _{q=-\infty }^{\infty }{u}_q{h}_q\left(t\right)\text{(26)}$$

Pour illustrer la formule de Shannon, voici la reconstruction à partir d'un signal échantillonné très simple :

import numpy as np
from matplotlib.pyplot import *

uq = [1,1,1,1,1]
Te = 1
N = 100000
t = np.linspace(-10*Te,10*Te,N)
u = np.zeros(N)
figure(figsize=(16,6))
for q in range(len(uq)):
    h = np.sinc(t/Te-q)
    u += uq[q]*h
    plot(t,uq[q]*h,'k--')
plot(t,u,'r-')
grid()
xlabel('t',fontsize=16)
ylabel('u(t)',fontsize=16)
                    

fig1.pdf

Les différents termes de la somme sont tracés en traits interrompus.

2.f. Reconstruction

La reconstruction d'un signal à partir des échantillons peut se faire quelle que soit la fréquence d'échantillonnage mais il faut que f_e>2f_m pour que la recontruction donne le signal u(t) analogique.

Il faut distinguer la reconstruction numérique de la reconstruction analogique (après conversion numérique-analogique).

La reconstruction numérique consiste à appliquer un algorithme qui permet d'obtenir la fonction u(t) à partir des échantillons. Obtenir numériquement u(t) signifie obtenir un échantillonnage de u(t) à une fréquence d'échantillonnage plus grande que f_e. Si f_e est à peine supérieure à 2f_m, il faudra augmenter la fréquence d'échantillonnage d'un facteur assez grand (au moins 10). Si f_e est grand devant 2f_m, un facteur plus petit sera suffisant; dans ce cas une reconstruction par interpolation linéaire peut d'ailleurs suffire. Si on fait la reconstruction avec la formule de Shannon $\mathrm{(24)}$, il faut échantillonner chaque terme de la somme (un sinus cardinal) avec une période d'échantillonnage beaucoup plus petite que T_e. La simulation ci-dessus montre un exemple. Le nombre de points de l'échantillonnage du sinus cardinal est N. Si P désigne le nombre d'échantillons u_q (5 dans l'exemple ci-dessus), il faut que N soit supérieur à P. La complexité temporelle de ce calcul est donc P², ce qui en fait une méthode très inefficace, seulement utilisable en pratique lorsque le nombre d'échantillons est faible (comme sur l'exemple précédent). On verra cependant plus loin une méthode de reconstruction par filtrage RIF, qui revient à appliquer la formule de Shannon mais avec une fonction h_q(t) limitée dans le temps.

Une autre manière de faire la reconstruction numérique est d'utiliser la relation $\mathrm{(20)}$, que l'on peut écrire de la manière suivante si $f_e>2f_m$ :

$$u\left(t\right)=T_e{\int }_{-f_n}^{f_n}{U}_{T_e}\left(f\right)exp\left(j2\pi ft\right)df\text{(27)}$$

Pour y parvenir, il faut calculer numériquement la TFTD $U_{T_e}\left(f\right)$. Cela est bien sûr impossible si le signal n'est pas à support borné mais dans la réalité les signaux étudiés sont à support borné (leur durée est celle de l'expérience). Si le signal u(t) est à support borné, la TFTD peut être calculée pour certaines fréquences (multiples de l'inverse de la durée du signal) au moyen de la TFD, qui elle-même se calcule très efficacement avec l'algorithme de transformée de Fourier rapide.

Pour reconstruire le signal, il faut qu'il soit à bande de fréquences limitée. Théoriquement, un signal à support borné ne peut être à bande de fréquences limitée. En pratique, même si le signal n'est pas à bande de fréquences limitée, il existe bien une fréquence f_m telle que sa TF est négligeable pour f>f_m. On pourra donc reconstruire le signal (avec une très faible erreur) si $f_e>2f_m$.

Supposons que le signal prenne des valeurs non nulles seulement dans l'intervalle [0,T₁]. La durée T₁ et la durée de l'expérience qui a permis de faire l'acquisition du signal (par conversion analogique-numérique). Le calcul de la transformée de Fourier discrète pour N₁ échantillons répartis uniformément sur l'intervalle [0,T₁] conduit aux valeurs de la TFTD pour N₁ fréquences multiples de 1/T₁. Lorsqu'on calcule la TFD afin de faire l'analyse spectrale du signal, on cherche à affiner l'échantillonnage en fréquence de la TFD en augmentant la durée et en remplissant par des zéros (on prolonge le signal échantillonné par des zéros). Soit T la durée finale, typiquement 5 fois plus grande que T₁ et telle que le nombre d'échantillons N soit une puissance de 2 (condition nécessaire pour l'algorithme FFT). Dans le cas présent, l'ajout de zéros est nécessaire pour obtenir un nombre d'échantillons N puissance ce 2. On a bien sûr $T=N{T}_e$. La TFD appliquée à ces échantillons donne la TFTD (de la fonction à support bornée) pour les fréquences suivantes :

$$f_k=\frac{k}{T},k=0,1\dots N-1\text{(28)}$$

On dispose ainsi d'un échantillonnage en fréquence de la fonction $U_{T_e}\left(f\right)$, qui doit permettre d'appliquer la formule $\mathrm{(27)}$ au moyen de la TFD inverse. Cependant, comme déjà remarqué, le nombre d'échantillons de u(t) qu'on souhaite obtenir doit être souvent beaucoup plus grand que le nombre d'échantillons N. La reconstruction est en effet une sorte d'interpolation. Pour obtenir cette augmentation du nombre d'échantillons, il faut étendre le spectre (la TFD) en le complétant par des zéros. En verra plus loin comment faire cet ajout de zéros. Le calcul de la TFD inverse de ce spectre étendu permet de reconstuire le signal, c'est-à-dire d'obtenir un échantillonnage du signal à une fréquence plus élevée que f_e. Le nombre de zéros à ajouter au spectre est d'autant plus grand que f_e est proche de 2f_m. L'intérêt de cette méthode, comparée à la précédente, est qu'elle repose sur l'algorithme de transformée de Fourier rapide.

La reconstruction analogique se fait au moyen d'un convertisseur numérique/analogique (N/A). La fréquence de la conversion N/A peut être égale à f_e (la fréquence d'échantillonnage du signal analogique) mais elle peut aussi être supérieure. Par exemple, dans le domaine audio, la fréquence d'échantillonnage est de 44 kHz mais il est fréquent que la conversion N/A se fasse à une fréquence 4 fois plus grande. Dans ce cas, il faut commencer par reconstruire les échantillons supplémentaires dûs à l'augmentation de fréquence d'échantillonnage. Il y a donc une première étape de reconstruction numérique, suivie de la conversion N/A. La sortie d'un convertisseur N/A délivre une tension égale à la valeur de l'échantillon (à une constante multiplicative près) et qui se maintient pendant une durée pratiquement égale à la période de conversion. On obtient donc un signal en marches d'escalier qu'il convient de filtrer pour recontruire le signal (ce filtrage n'est pas toujours nécessaire). Le filtre passe-bas idéal qui convient est un filtre passe-bas idéal de fréquence de coupure f_e/2, où f_e désigne la fréquence d'échantillonnage initiale (avant l'augmentation). En pratique, on doit évidemment se contenter d'un filtre passe-bas analogique loin d'être idéal. Étant donné la difficulté et le coût de réalisation de filtres passe-bas très sélectifs, la méthode généralement préférée (lorsqu'elle est possible) consiste à se passer de filtre analogique et à augmenter la fréquence de conversion en faisant une reconstruction numérique. On verra plus loin comment une reconstruction numérique peut être réalisée dans un traitement en temps réel, au moyen d'un filtre RIF.

3.a. Échantillonnage du signal

L'exemple développé ci-après consiste à échantillonner un signal défini par une fonction. Il s'agit de simuler l'acquisition d'un signal sinusoïdal pendant une durée finie. On doit donc définir une fonction sinusoïdale multipliée par une fonction de fenêtrage (qui limite la durée de la sinusoïde). Dans la plupart des cas, la fonction de fenêtrage qui est par défaut réalisée lors d'une expérience est une fenêtre rectangulaire. Cela correspond à un convertisseur A/N qui est mis en fonctionnement à l'instant t=0 puis est stoppé à l'instant t=T. On nommera période celle de la sinusoïde bien que le signal considéré ne soit pas périodique (on peut le qualifier de quasi périodique). Si la période est très petite devant T, une fenêtre rectangulaire peut convenir pour l'analyse spectrale. Sinon, il peut être intéressant d'appliquer une fonction de fenêtrage différente, afin de réduire les lobes secondaires de la raie du spectre.

f1 = 1.12 # fréquence de la sinusoïde
fe = 10
Te = 1/fe
T = 100 # durée du signal 
N = int(T/Te)
t = np.arange(N)*Te 
u = np.sin(2*np.pi*f1*t)
figure(figsize=(16,6))
plot(t,u,'.r')
xlabel('t',fontsize=16)
ylabel('u')
xlim(0,5)
                  

fig2.pdf

3.b. Transformée de Fourier à temps discret

L'évaluation numérique de la transformée de Fourier à temps discret (TFTD) consiste à calculer la TFTD pour un échantillonnage de la fréquence dans l'intervalle [0,f_e]. Rappelons que la TFTD est périodique, de période f_e. Pour obtenir un échantillonnage assez fin, il faut généralement étendre le signal échantillonné par des zéros (remplissage par des zéros, en anglais zero padding), en faisant en sorte que le nombre final d'échantillons soit une puissance de 2, afin de pouvoir appliquer pleinement l'algorithme de transformée de Fourier discrète. Lorsque le but est de faire l'analyse spectral du signal, on cherche en général le spectre du signal avant fenêtrage, ce qui est en principe impossible puisque justement on ne dispose pas de ce signal. Cependant, on peut s'en approcher si la fenêtre est assez longue (la durée T assez grande) et si une fonction de fenêtrage est appliquée (voir Analyse spectrale d'un signal numérique).

Dans le cas présent, on s'intéresse à la TFTD du signal échantillonné avec le fenêtrage rectangulaire (c'est-à-dire à la TFTD du signal sinusoïdal multiplié par la fenêtre rectangulaire), donc on n'applique pas de fonction de fenêtrage.

from numpy.fft import fft,ifft

nz = 10 # nombre de blocs de zéros ajoutés
p = int(np.log(nz*N)/np.log(2))+1
Np = 2**p
u1 = np.concatenate((u,np.zeros(Np-N)))
t1 = np.arange(Np)*Te
tfd = fft(u1)/N # transformée de Fourier discrète
freq = np.arange(Np)*1/(Np*Te)
figure(figsize=(16,6))
plot(freq,np.absolute(tfd))
xlabel('f',fontsize=16)
ylabel(r"$U_{T_e}$",fontsize=16)
grid()
                 

fig3.pdf

On obtient ainsi une représentation graphique du module de la TFTD sur l'intervalle [0,f_e], qui se prolonge par périodicité en dehors de cet intervalle. La fenêtre rectangulaire se manifeste par une raie en forme de sinus cardinal. En effet, la transformée de Fourier du signal fenêtré est égale au produit de convolution de la transformée de Fourier du signal sinusoïdal (distribution de Dirac formée d'un raie à f₁ et d'une raie à -f₁) et de celle de la fenêtre rectangulaire (donnée plus haut). Voici un détail de la raie :

xlim(f1-0.3,f1+0.3)                  
                   

fig4.pdf

La courbe que l'on voit sur cette figure correspond à la valeur absolue d'un sinus cardinal centré en f=f₁. La largeur de son lobe principal est 2/T.

Il faut remarquer que la fonction sinus cardinal n'est évidemment pas à support borné (aucune transformée de Fourier de fonction de fenêtrage ne l'est). En conséquence, la TF du signal fenêtré n'est pas à bande de fréquence limitée (contrairement au signal sinusoïdal non fenêtré). Il s'en suit que le respect de la condition de Nyquist-Shannon ($f_e>2f_m$) est en principe impossible pour le signal fenêtré. Cependant, la décroissance de la fonction sinus cardinal (relativement lente) fait qu'un compromis est possible : la fréquence d'échantillonnage doit être assez grande pour que l'intensité des lobes secondaires du sinus cardinal à la fréquence f_e/2 soit négligeable. Cette condition est beaucoup plus facile à satisfaire avec les fenêtres de Hamming, Hann et Blackman car leur transformée de Fourier décroît beaucoup plus vite que celle de la fenêtre rectangulaire. Lorsqu'on emploie une fenêtre rectangulaire, il est préférable d'opter pour une fréquence d'échantillonnage beaucoup plus grande que 2f (ici, nous avons un facteur 5). Malgré cela, il y a un léger repliement de bande puisque le support de la TF s'étend au delà de f_e/2. En théorie, il sera donc impossible de faire une reconstruction parfaite du signal fenêtré (mais cela est vrai quelle que soit la fonction de fenêtrage). En pratique, l'effet du repliement devrait être négligeable.

3.c. Reconstruction par interpolation linéaire

Cette reconstruction consiste à faire une interpolation linéaire entre deux échantillons consécutifs. La fonction matplotlib.plot réalise graphiquement cette interpolation :

figure(figsize=(16,6))
plot(t,u,'.r')
plot(t,u,'k-')
xlabel('t',fontsize=16)
ylabel('u')
xlim(0,5)
                 

fig5.pdf

Il est évident que cette reconstruction n'est satisfaisante que si la fréquence d'échantillonnage est beaucoup plus grande que la fréquence de la sinusoïde, ce qui est loin d'être le cas sur cet exemple. On peut cependant remarquer que si l'on compte les cycles sur une durée donnée, on peut déterminer la fréquence f. Même si la forme d'onde est fausse, la fréquence est correcte (car f_e>2f).

3.d. Reconstruction par la transformée de Fourier discrète inverse

Nous disposons de la TFTD du signal fenêtré pour un échantillonnage de fréquences. Voici ce que donne la TFD inverse (il faut multiplier par N) :

u1 = np.real(ifft(tfd))*N # partie imaginaire nulle en théorie mais pas tout à fait en réalité
figure(figsize=(16,6))
plot(t1,u1,'k-')
xlabel('t',fontsize=16)
ylabel('u')
                

fig6.pdf

Sans surprise, on obtient le signal complété par des zéros. Voici un détail :

xlim(0,5)           
                

fig7.pdf

On retrouve évidemment les échantillons du signal échantillonné. Pour faire une reconstruction plus fine, il faut étendre la TFTD par des zéros. Le bloc de zéros doit être centré en f_e/2. Afin que le nombre total d'éléments reste une puissance de 2, ajoutons Np zéros :

def ajout_zero_tfd(tfd):
    Np= len(tfd)
    return np.concatenate((tfd[0:Np//2],np.zeros(Np),tfd[Np//2:Np]))
    
tfd2 = ajout_zero_tfd(tfd)
freq2 = np.arange(len(tfd2))*1/(Np*Te)

figure(figsize=(16,6))
plot(freq2,np.absolute(tfd2))
xlabel('f',fontsize=16)
ylabel(r"$U_{T_e}$",fontsize=16)
grid()
                 

fig8.pdf

Comme prévu, la fréquence d'échantillonnage est multipliée par deux.

Voici le signal recontruit par transformée discrète inverse de cette TFD :

u2 = np.real(ifft(tfd2))*N 
t2 = np.arange(len(u2))*Te/2
figure(figsize=(16,6))
plot(t2,u2,'k-')
xlabel('t',fontsize=16)
ylabel('u')               
                 

fig9.pdf

xlim(0,5)           
                

fig10.pdf

On peut encore multiplier par deux la fréquence d'échantillonnage :

tfd3 = ajout_zero_tfd(tfd2)
freq3 = np.arange(len(tfd3))*1/(Np*Te)
u3 = np.real(ifft(tfd3))*N 
t3 = np.arange(len(u3))*Te/4
figure(figsize=(16,6))
plot(t3,u3,'k-')
xlabel('t',fontsize=16)
ylabel('u')
                

fig11.pdf

xlim(0,5)           
                

fig12.pdf

xlim(T-5,T)           
                

fig13.pdf

On obtient ainsi une reconstruction satisfaisante. Il faut cependant remarquer que la reconstruction n'est pas parfaite au début et à la fin du signal. C'est une conséquence du fait que la TF de la fenêtre rectangulaire décroît relativement lentement (fonction sinus cardinal) et que la condion de Nyquist-Shannon ne peut être satisfaite parfaitement pour cette fenêtre.

La reconstruction ne nécessite pas de compléter le signal initial par des zéros (ce qui était nécessaire pour obtenir une bonne représentation de la TFTD) mais il est tout de même intéressant d'ajouter le nombre de zéros minimale nécessaire pour obtenir une taille puissance de 2 :

nz = 1 # nombre de blocs de zéros ajoutés, un seul suffit pour la reconstruction
p = int(np.log(nz*N)/np.log(2))+1
Np = 2**p
u1 = np.concatenate((u,np.zeros(Np-N)))
t1 = np.arange(Np)*Te
tfd = fft(u1)/N # transformée de Fourier discrète  
tfd2 = ajout_zero_tfd(tfd)
tfd3 = ajout_zero_tfd(tfd2)
freq3 = np.arange(len(tfd3))*1/(Np*Te)
u3 = np.real(ifft(tfd3))*N 
t3 = np.arange(len(u3))*Te/4
figure(figsize=(16,6))
plot(t3,u3,'k-')
xlabel('t',fontsize=16)
ylabel('u')
xlim(0,5)
                        

fig14.pdf

3.e. Reconstruction par un filtre RIF

Lorsque la reconstruction doit se faire en temps réel (typiquement en amont d'un convertisseur N/A), la méthode par transformée de Fourier inverse est difficile à mettre en œuvre. Dans ce cas, il est beaucoup plus simple de mettre en place un filtrage par convolution pour faire la reconstruction. Un filtre de convolution est un filtre à réponse impulsionnelle finie (RIF). Si l'on souhaite une très forte sélectivité du filtre, donc une réponse impulsionnelle très longue, une méthode par transformée de Fourier discrète est plus efficace (voir Filtres RIF : filtrage par convolution et par transformée de Fourier rapide) mais elle est beaucoup plus difficle à implémenter.

Pour faire une reconstruction par filtrage RIF, on doit commencer par augmenter la fréquence d'échantillonnage en dupliquant les échantillons déjà existants. Voici par exemple une augmentation d'un facteur 4 de la fréquence d'échantillonnage (facteur utilisé pour la conversion audio) :

u4 = np.zeros(len(u)*4)
for k in range(4):
    u4[k::4] = u
t4 = np.arange(len(u4))*(Te/4)
figure(figsize=(16,6))
plot(t4,u4,'r.')
xlabel('t',fontsize=16)
ylabel('u')
xlim(0,5)
                    

fig15.pdf

Le signal échantillonné qu'on obtient ainsi correspond au signal échantillonné à la sortie d'un convertisseur N/A (avec une fréquence d'échantillonnage 4 fois plus grande que la fréquence de conversion). Ce qui suit constitue donc aussi une simulation d'un filtrage passe-bas analogique en aval du convertisseur N/A.

Le filtre RIF doit avoir en principe une fréquence de coupure égale à f_e/2 (il s'agit de la fréquence d'échantillonnage initiale). Ces coefficients peuvent être calculée par la méthode du fenêtrage :

from scipy.signal import firwin,freqz
P = 100
M = 2*P+1
fc = fe/2 # fréquence de coupure
h = firwin(numtaps=M,cutoff=[fc/(4*fe)],window='hann',nyq=0.5)
                      

Il faut remarquer que la fréquence d'échantillonnage du signal traité est 4 fois la fréquence initiale.

Voici la réponse fréquentielle de ce filtre :

w,H=freqz(h,1000)
figure(figsize=(12,8))
subplot(211)
plot(w/(2*np.pi),20*np.log10(np.absolute(H)))
xlabel("f/fe")
ylabel("GdB")
grid()
subplot(212)      
plot(w/(2*np.pi),np.unwrap(np.angle(H)))
xlabel("f/fe")
ylabel("phase")
grid()
                      

fig16.pdf

Voici le résultat du filtrage RIF par convolution :

from scipy.signal import convolve
u5 = convolve(u4,h,mode='same')
figure(figsize=(16,6))
plot(t4,u5,'r.')
xlabel('t',fontsize=16)
ylabel('u')
xlim(0,5)

                      

fig17.pdf

Le résultat est très semblable à celui obtenu précédemment par transformée de Fourier discrète inverse (pour le même facteur d'augmentation de la fréquence d'échantillonnage). La méthode par TFD inverse est plus efficace, donc elle est préférable si on doit traiter un grand nombre de signaux. Elle l'est aussi si on doit traiter un seul signal très long mais son implémentation dans ce cas est beaucoup plus complexe (voir Filtres RIF : filtrage par convolution et par transformée de Fourier rapide).

La reconstruction par un filtre RIF revient à appliquer la formule de Shannon mais avec une fonction h_q(t) tronquée et multipliée par une fonction de fenêtrage. Le fait de tronquer la réponse impulsionnelle réduit considérablement le temps de calcul. La multiplication par une fonction de fenêtrage (Hann dans l'exemple ci-dessus) est nécessaire pour éviter les ondulations du gain dans la bande passante.