Accueil > Unclassified > Discrete Fourier transform: Fourier transform
Unclassified

Discrete Fourier transform: Fourier transform

1. Fourier transform

Let be a signal u (t). Its Fourier transform (TF) is: S (f) = ∫-∞∞u (t) exp (-j2πft) dt

If u (t) is real, its Fourier transform has the following parity: S (-f) = S (f) *

The signal is expressed with its TF by the inverse Fourier transform: u (t) = ∫-∞∞S (f) exp (j2πft) df

During the digital processing of a signal, we have u (t) over a duration T, for example over the interval [-T / 2, T / 2]. An approximation of the TF is calculated in the form: Sa (f) = ∫-T / 2T / 2u (t) exp (-j2πft) dt

Let be a sampling of N points, obtained for 0≤k≤N-1: tk = -T2 + kTNuk = u (tk)

An approximation is obtained by the method of rectangles: Sa (f) ≃TNexp (jπfT) ∑k = 0N-1ukexp (-j2πkfTN)

We look for the TF for the following frequencies, with 0≤n≤N-1: fn = nT

that is to say: Sa (fn) = Texp (jπn) 1N∑k = 0N-1ukexp-j2πknN

By denoting Sn the discrete Fourier transform (DFT) of uk, we therefore have: Sa (fn) ≃Texp (jπn) Sn

In a spectral analysis, we are generally interested in the modulus of S (f), which allows to ignore the term exp (jπ n)

The spectrum obtained is by nature discrete, with lines spaced 1 / T apart. It is therefore the spectrum of a periodic signal of period T. To simulate a continuous spectrum, T must be chosen large enough.

The spectrum obtained is periodic, of periodicity fe = N / T, the sampling frequency.

2. Signal with bounded support

2.a. Example: Gaussian

We choose T such that u (t) = 0 for | t |> T / 2. Consider for example a Gaussian centered at t = 0: u (t) = exp-t2a2

whose Fourier transform is S (f) = aπexp (-π2a2f2)

By choosing for example T = 10a, we have | u (t) | <10-10 for 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

The following function calculates the TFD for a duration T and a number of 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];
]
            

For comparison, we will also plot the theoretical TF of the Gaussian:

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

Example :

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

The spectrum obtained must be extended by periodicity of period N / T (5 in the example). To increase the resolution, you need to increase T. It is interesting to keep the N / T sampling frequency constant:

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

The sampling frequency must verify Shannon’s criterion: N / T must be greater than 2 times the greatest frequency present in the spectrum. Here is an example where this criterion is not verified:

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. Example: sinusoid modulated by a Gaussian

We consider the following signal: 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

In this case, it is necessary to choose a sampling frequency higher than twice the frequency of the sine wave, ie 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. Rectangular window

Consider a rectangular window of width a:

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

Let’s see its spectrum with T = 100 and 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

Unbounded medium signal

In this case, the window [-T / 2, T / 2] is arbitrarily imposed by the measurement system. For example on a digital oscilloscope, T can be adjusted by setting the time base. Consider for example a periodic signal comprising 3 harmonics:

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

The sampling frequency should check 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

We recognize in the spectrum the lines corresponding to the 3 harmonics. The Fourier coefficients are found by dividing by T. Let us modify the frequency of the 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

The base of the lines is notably wider. The sampled signal is in fact the product of the periodic signal defined above by a rectangular window h (t) of width T. The TF is therefore the convolution product of S with the TF of h: H (f) = Tsin ( πTf) πTf

which exhibits slowly decreasing oscillations which can lead to interference between the different harmonics of the signal. To remedy this problem, we replace the rectangular window by a window whose spectrum has weaker side lobes, for example the Hamming window:

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

It may interest you

Leave a Reply

Your email address will not be published.

Solve : *
19 − 12 =