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"}]
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]
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]
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]
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"}]
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]}]
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]}]
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]}]
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]}]
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]}]
