Site icon Lab4Sys.com

Transformée de Fourier discrète : transformée de Fourier

1. Transformée de Fourier

Soit un signal u(t). Sa transformée de Fourier(TF) est :S(f)=∫-∞∞u(t)exp(-j2πft)dt

Si u(t) est réel, sa transformée de Fourier possède la parité suivante :S(-f)=S(f)*

Le signal s’exprime avec sa TF par la transformée de Fourier inverse :u(t)=∫-∞∞S(f)exp(j2πft)df

Lors du traitement numérique d’un signal, on dispose de u(t) sur une durée T, par exemple sur l’intervalle [-T/2,T/2]. Une approximation de la TF est calculée sous la forme :Sa(f)=∫-T/2T/2u(t)exp(-j2πft)dt

Soit un échantillonnage de N points, obtenu pour 0≤k≤N-1 :tk=-T2+kTNuk=u(tk)

Une approximation est obtenue par la méthode des rectangles :Sa(f)≃TNexp(jπfT)∑k=0N-1ukexp(-j2πkfTN)

On recherche la TF pour les fréquences suivantes, avec 0≤n≤N-1 :fn=nT

c’est-à-dire :Sa(fn)=Texp(jπn)1N∑k=0N-1ukexp-j2πknN

En notant Sn la transformée de Fourier discrète (TFD) de uk, on a donc :Sa(fn)≃Texp(jπn)Sn

Dans une analyse spectrale, on s’intéresse généralement au module de S(f), ce qui permet d’ignorer le terme exp(jπ n)

Le spectre obtenu est par nature discret, avec des raies espacées de 1/T. C’est donc le spectre d’un signal périodique de période T. Pour simuler un spectre continu, T devra être choisi assez grand.

Le spectre obtenu est périodique, de périodicité fe=N/T, la fréquence d’échantillonnage.

2. Signal à support borné

2.a. Exemple : gaussienne

On choisit T tel que u(t)=0 pour |t|>T/2. Considérons par exemple une gaussienne centrée en t=0 :u(t)=exp-t2a2

dont la transformée de Fourier estS(f)=aπexp(-π2a2f2)

En choisissant par exemple T=10a, on a |u(t)|<10-10 pour t>T/2

a=1;
signal[t_]:=Exp[-t^2/a^2];
            
Plot[signal[t],{t,-5,5},PlotRange->{{-5,5},{0,1}},AxesLabel->{"t","s"}]
plotA.png

La fonction suivante calcule la TFD pour une durée T et un nombre points N

spectre[T_,npoints_]:=Module[{te,temps,echantillons,transFourier,frequences,spectre},
	te=T/npoints;
	temps[k_]:=-T/2+(k-1)*te;
	echantillons:=Table[signal[temps[k]],{k,npoints}];
	transFourier=Fourier[echantillons,FourierParameters->{-1,-1}]*T;
	frequences[k_]:=(k-1)/T;
	spectre=Table[{frequences[k],Abs[transFourier[[k]]]},{k,npoints}];
	Return[spectre];
]
            

Pour comparaison, on tracera aussi la TF théorique de la gaussienne :

tfSignal[f_]:=a*Sqrt[Pi]*Exp[-(Pi*a*f)^2];
            

Exemple :

T=20; npoints=100;
sp=spectre[T,npoints];
p1=ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{-npoints/T,npoints/T},{0,2}},
									PlotStyle->{RGBColor[1,0,0]}];
p2=Plot[tfSignal[f],{f,-npoints/T,npoints/T},PlotRange->{{-npoints/T,npoints/T},{0,2}}];
            
Show[p1,p2]
plotB.png

Le spectre obtenu doit être prolongé par périodicité de période N/T (5 sur l’exemple). Pour augmenter la résolution, il faut augmenter T. Il est intéressant de maintenir constante la fréquence d’échantillonnage N/T :

T=100; npoints=500;
sp=spectre[T,npoints];
p1=ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{-npoints/T,npoints/T},{0,2}},
									PlotStyle->{RGBColor[1,0,0]}];
p2=Plot[tfSignal[f],{f,-npoints/T,npoints/T},PlotRange->{{-npoints/T,npoints/T},{0,2}}];
            
Show[p1,p2]
plotC.png

La fréquence d’échantillonnage doit vérifier le critère de Shannon : N/T doit être supérieur à 2 fois la plus grande fréquence présente dans le spectre. Voici un exemple où ce critère n’est pas vérifié :

T=100; npoints=100;
sp=spectre[T,npoints];
p1=ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{-npoints/T,npoints/T},{0,2}},
									PlotStyle->{RGBColor[1,0,0]}];
p2=Plot[tfSignal[f],{f,-npoints/T,npoints/T},PlotRange->{{-npoints/T,npoints/T},{0,2}}];
            
Show[p1,p2]
plotD.png

2.b. Exemple : sinusoïde modulée par une gaussienne

On considère le signal suivant :u(t)=exp(-t2/a2)cos(2πtb)

avec b≪a .

a=1;
b=0.1; 
signal[t_]:=Exp[-t^2/a^2]*Cos[2*Pi*t/b];
                
Plot[signal[t],{t,-5,5},PlotRange->{{-2,2},{-1,1}},
                        AxesLabel->{"t","s"}]
plotE.png

Dans ce cas, il faut choisir une fréquence d’échantillonnage supérieure à 2 fois la fréquence de la sinusoïde, c.a.d. N/T>2/b.

T=100; npoints=3000;
sp=spectre[T,npoints]; 
                
ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{0,npoints/T},{0,2}},
                        PlotStyle->{RGBColor[1,0,0]}]
plotF.png

2.c. Fenêtre rectangulaire

Soit une fenêtre rectangulaire de largeur a :

a=1;
signal[t_]:=If[Abs[t]<(a/2),1,0];
                

Voyons son spectre avec T=100 et N=1000 :

T=100; npoints=3000;
sp=spectre[T,npoints]; 
                
ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{0,npoints/T},{0,2}},
                        PlotStyle->{RGBColor[1,0,0]}]
plotG.png

3. Signal à support non borné

Dans ce cas, la fenêtre [-T/2,T/2] est arbitrairement imposée par le système de mesure. Par exemple sur un oscilloscope numérique, T peut être ajusté par le réglage de la base de temps. Considérons par exemple un signal périodique comportant 3 harmoniques :

b=1; (*periode*)
w0=2*Pi/b;
signal[t_]:=Cos[w0*t]+0.5*Cos[2*w0*t]+0.1*Cos[3*w0*t];
            

La fréquence d’échantillonnage doit vérifier N/T>6/b.

T=100; npoints=1000;
sp=spectre[T,npoints];
            
ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{0,npoints/T},{0,T}},
                    PlotStyle->{RGBColor[1,0,0]}]
plotH.png

On reconnait dans le spectre les raies correspondant aux 3 harmoniques. Les coefficients de Fourier se retrouvent en divisant par T. Modifions la fréquence du signal :

b=0.9;
w0=2*Pi/b;
signal[t_]:=Cos[w0*t]+0.5*Cos[2*w0*t]+0.1*Cos[3*w0*t];
T=100; npoints=1000;
sp=spectre[T,npoints];
            
ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{0,npoints/T},{0,T}},PlotStyle->{RGBColor[1,0,0]}]
plotI.png

La base des raies est notablement plus large. Le signal échantillonné est en fait le produit du signal périodique défini ci-dessus par une fenêtre h(t) rectangulaire de largeur T. La TF est donc le produit de convolution de S avec la TF de h :H(f)=Tsin(πTf)πTf

qui présente des oscillations lentement décroissantes qui peuvent conduire à des interférences entre les différentes harmoniques du signal. Pour remedier à ce problème, on remplace la fenêtre rectangulaire par une fenêtre dont le spectre présente des lobes secondaires plus faibles, par exemple la fenêtre de Hamming :

b=0.9;
w0=2*Pi/b;
hamming[t_]:=0.54+0.46*Cos[2*Pi*t/T];
signal[t_]:=(Cos[w0*t]+0.5*Cos[2*w0*t]+0.1*Cos[3*w0*t])*hamming[t];
T=100; npoints=1000;
sp=spectre[T,npoints];
            
ListPlot[sp,AxesLabel->{"f","S"},PlotRange->{{0,npoints/T},{0,T}},
                    PlotStyle->{RGBColor[1,0,0]}]
plotJ.png
Quitter la version mobile