1. TFD of a sampled signal
The spectrum of an analog signal is obtained with its Fourier transform. The spectrum of a digital signal is obtained by its discrete Fourier transform, according to the method described in TFD applied to the Fourier transform.
Consider as an example a periodic signal (period 1):
signal[t_]:=Cos[2*Pi*t]+0.5*Cos[2*2*Pi*t]+0.25*Sin[3*2*Pi*t]+0.125*Sin[4*2*Pi*t];
To sample this signal, we set a number of points and a total duration which does not coincide with a period of the signal, as is the case in the digitizations of real signals:
tmax=50.32563;
n=1000;
te=tmax/n;
fe=1/te (*frequence d'echantillonnage*)
19.870590790418323
echant=Table[signal[(k-1)*te],{k,n}];
We calculate the TFD and we construct the corresponding frequency scale:
tfd=Fourier[echant,FourierParameters->{-1,-1}];
spectre=Table[{(k-1)/tmax,Abs[tfd[[k]]]},{k,n}];
ListPlot[spectre,Filling->Axis,PlotRange->{{0,5},Full},AxesLabel->{"f","A"}]
If the signal were sampled over an infinite time, the lines would be zero width (since the signal is periodic). The broadening of the base of the lines is an effect of the finite duration of the sample, that is to say of the application of a rectangular window to the signal. Errors are also observed in the relative heights of the lines, which come from the fact that the frequency resolution (inverse of the duration of the sample) is insufficient to capture the maximum of the lines.
2. Truncation windows
Let u (t) be the signal analyzed over the interval [0, T]. The truncation window w (t) is a null function outside this interval. The spectral analysis then consists in calculating the following Fourier transform: S (f) = ∫0Tu (t) w (t) exp (-i2πft) dt (1)
In the previous example, the truncation window is rectangular (w (t) = 1).
If we denote by U (f) the Fourier transform of u (t) and W (f) that of the window w (t), the calculated Fourier transform is the convolution product: S (f) = ∫- ∞ + ∞U (fx) W (x) dx (2)
The Fourier transform of the rectangular window is: W (f) = Texp (-iπfT) sinc (πfT) (3)
To characterize the effect of the window, we plot the modulus of its Fourier transform in decibel (by convention T = 1):
Wr[f_]:=Sinc[Pi*f];
Plot[20*Log[10,Abs[Wr[f]]],{f,-10,10},AxesLabel->{"f","dB"},
PlotRange->{Full,{-100,0}}]
The Hamming window is defined by: w (t) = 0.54-0.5cos2πTt (4)
hamming[t_]:=Piecewise[{{0,t<0},{0,t>1}},0.54-0.5*Cos[2*Pi*t]];
Plot[hamming[t],{t,0,1}]
We calculate its Fourier transform then we plot its modulus in decibel:
Wh=FourierTransform[hamming[t],t,f,FourierParameters->{0,-2*Pi}];
Plot[20*Log[10,Abs[Wh]],{f,-10,10},AxesLabel->{"f","dB"},PlotRange->{Full,{-100,0}}]
The first secondary lobe on either side of the principal maximum is much lower than for the rectangular window. In contrast, the main lobe has a width at -3dB greater. Let’s see the effect of this window on the DFT of the signal:
echantHamming = Table[echant[[k]]*hamming[k/n],{k,n}];
tfd=Fourier[echantHamming,FourierParameters->{-1,-1}];
spectre=Table[{(k-1)/tmax,Abs[tfd[[k]]]},{k,n}];
ListPlot[spectre,Filling->Axis,PlotRange->{{0,5},Full},AxesLabel->{"f","A"}]
The lines no longer show the widening of their base observed with a rectangular window. This is a consequence of the much lower value of the sidelobes. This property is important for obtaining continuous spectra. The height of the lines is better restored because the central lobe is wider: for a given total duration, there are therefore more points in the central lobe. When possible, the duration T can of course be increased to refine the precision of the lines.
By comparison, the rectangular window has a good frequency resolution (narrower central lobe) but does not correctly reproduce the continuous spectra because of the too intense secondary lobes.
3. Mathematica calculation function
The following function calculates the spectrum of a discrete signal. The sampling frequency is provided as an argument. The type of window can be specified as an option (hamming by default):
Options[Spectre]={Fenetre->"hamming"};
Spectre[signal_,fe_,opts___]:=Module[{ne,fenetre,w,echant,tfd,spectre},
fenetre = Fenetre/.{opts}/.Options[Spectre];
ne=Length[signal];
w["rectangulaire"][k_]:=1.0;
w["hamming"][k_]:=0.54-0.46*Cos[2*Pi/ne*k];
w["hanning"][k_]:=0.5-0.5*Cos[2*Pi/ne*k];
w["blackman"][k_]:=0.42-0.5*Cos[2*Pi/ne*k]+0.08*Cos[4*Pi/ne*k];
echant=Table[signal[[k+1]]*w[fenetre][k],{k,0,ne-1}];
tfd=Fourier[echant,FourierParameters->{-1,-1}];
spectre=Table[{fe*k/ne,Abs[tfd[[k]]]},{k,ne}];
Return[spectre];
];
Example with the signal defined above:
spectre=Spectre[echant,fe,Fenetre->"hanning"];
ListPlot[spectre,Filling->Axis,PlotRange->{{0,5},Full},AxesLabel->{"f","A"}]
