Accueil > Unclassified > Discrete Fourier transform: Fourier series
Unclassified

Discrete Fourier transform: Fourier series

1. Fourier coefficients

Let u (t) be a signal of period T developable in Fourier series: u (t) = ∑n = -∞∞cnexpj2πnTt

The Fourier coefficients are, for n∈Z: cn = 1T∫0Tu (t) exp-j2πnTtdt

We consider a sampling of u (t) of N points, with 0≤k≤N-1: tk = kTNuk = u (tk)

An approximation of the Fourier coefficients can be obtained by the method of rectangles: cn≃1T∑k = 0N-1ukexp-j2πnkNTN

The formula obtained defines a sequence of period N. It is therefore sufficient to calculate, for 0≤n≤N-1: Sn = 1N∑k = 0N-1ukexp-j2πnkN

The application which associates the sequence of N numbers uk with the sequence Sn is the discrete Fourier transform (DFT). Sn is an approximation of the Fourier coefficient cn, corresponding to the frequency harmonic: fn = nT

2. TFD with Mathematica

The function allowing to calculate the DFT (by the fast Fourier transform algorithm) is: Fourier [{s0, s1, …, sN-1}, FourierParameters -> {a, b}]

which calculates the following sum (index n from 1 to N): 1N (1-a) / 2∑k = 1Nukexpj2πb (n-1) k-1N

To obtain the DFT defined above, it is therefore necessary to set a = -1 and b = -1.

3. Applications

3.a. Trigonometric polynomial

We define a signal in the form of a trigonometric polynomial of period T:

T = 1;
w0 = 2*Pi*T;
signal[t_] := 1 + 2*Cos[w0*t] + 0.5*Cos[2*w0*t] + 1*Sin[3*w0*t] + 0.2*Cos[10*w0*t];
                
Plot[signal[t], {t, 0, 1},AxesLabel->{'t','s'}]
plotA.png

Signal sampling:

npoints = 50; (* number of points *)
te = T/npoints; (* sampling time *)
time[k_] := (k - 1)*te;
samples = Table[signal[time[k]], {k, npoints}];
                

Calculation of the DFT and plot of the spectrum (modulus of | Sn |):

coefFourier = Fourier[echantillons, FourierParameters -> {-1, -1}];
frequences[k_] := (k - 1)/T;
raies = Table[Line[{{frequences[k], 0}, {frequences[k], Abs[coefFourier[[k]]]}}], {k, npoints}];           
                
Show[Graphics[{RGBColor[1,0,0],raies},Frame->True,PlotRange->{{0,npoints/T},{0,1}},
			AspectRatio->0.6,FrameLabel->{'f','|Sn|'}]]
plotB.png

We recognize on the first half of the spectrum the Fourier coefficients c0 = 1, c1 = 2/2, c2 = 0.5 / 2, c3 = -j / 2, c10 = 0.2 / 2

We notice the following property: | Sn | = | SN-n |

Indeed, if the function u (t) has real values, we have: SN-n = Sn *

The second part of the spectrum (also called the image of the spectrum) therefore corresponds to the Fourier coefficients of negative indices: c-n≃SN-n

There is also the periodicity of the sequence Sn: Sn = SN + n

We therefore obtain a spectrum whose periodicity is the sampling frequency fe = NT

It is effectively the spectrum of a signal sampled at this frequency.

To simplify, let’s define a function that calculates the spectrum for a given sample rate:

raiesSpectre[npoints_]:=Module[{te,temps,echantillons,frequences},
	te = T/npoints;
	temps[k_] := (k - 1)*te;
	echantillons = Table[signal[temps[k]], {k, npoints}];
	coefFourier = Fourier[echantillons, FourierParameters -> {-1, -1}];
	frequences[k_] := (k - 1)/T;
	Return[Table[Line[{{frequences[k], 0}, {frequences[k], Abs[coefFourier[[k]]]}}],
						{k, npoints}]];
]
                

According to Shannon’s theorem, the sampling frequency must be greater (strictly) than twice the largest frequency present in the signal spectrum (10 times 2 on the previous example). Let’s see what happens if this condition is not true:

npoints=15;
raies = raiesSpectre[npoints];
                
Show[Graphics[{RGBColor[1,0,0],raies},Frame->True,PlotRange->{{0,npoints/T},{0,1}},
			AspectRatio->0.6,FrameLabel->{'f','|Sn|'}]]
plotC.png

We notice the appearance of a harmonic of order n = 5, a position which in fact corresponds to the image of the harmonic of order 10. This phenomenon is called band aliasing. It occurs when a significant part of the spectrum is superimposed on its image.

Let us consider the case of a periodic signal comprising a maximum frequency in its spectrum, that is to say that there exists p∈N such that cn = 0 for | n |> p. In this case, the signal is a trigonometric polynomial. The condition imposed by Shannon’s theorem is sufficient to obtain the Fourier coefficients of u (t), up to rounding errors. Let’s check it on the previous example by adopting N = 21, i.e. f_e = 21> 20.

npoints=21;
raies = raiesSpectre[npoints];
                
Show[Graphics[{RGBColor[1,0,0],raies},Frame->True,PlotRange->{{0,npoints/T},{0,1}},
			AspectRatio->0.6,FrameLabel->{'f','|Sn|'}]]
plotD.png

The 10th order harmonic is correctly calculated.

3.b. Square signal

This example highlights the problems associated with the existence of high frequencies in the spectrum. We define the signal on the interval [0, T]:

T=1
signal[t_]:=If[t<T/2,1,-1];     
                

In this case, there is no maximum frequency in the spectrum. The precision of the calculation will be all the greater as N is large.

npoints=100;
raies = raiesSpectre[npoints];
                
Show[Graphics[{RGBColor[1,0,0],raies},Frame->True,PlotRange->{{0,npoints/T},{0,1}},
			AspectRatio->0.6,FrameLabel->{'f','|Sn|'}]]
plotE.png

The relatively large value of the harmonics in the vicinity of fe / 2 suggests a non-negligible band aliasing phenomenon. In theory, the odd Fourier coefficients of the slot function verify: | c1 || cn | = n

The following table makes it possible to determine the precision of the calculation by the TFD:

indices=Table[2*p-1,{p,1,20}];
precision=Table[Abs[coefFourier[[2]]]/Abs[coefFourier[[2*p]]]/(2*p-1),{p,1,20}];
            
n|S1|/(n|Sn|)
11
30.998684
50.996057
70.992122
90.986892
110.980376
130.972592
150.963556
170.953292
190.941822
210.929174
230.915377
250.900464
270.884471
290.867433
310.849391
330.830387
350.810465
370.789671
390.768054

We see that the error exceeds 5 percent from the rank n = 19. This is due to the aliasing of harmonics of order higher than fe / 2 (50) on harmonics of lower order. To improve accuracy, let’s increase N:

npoints=500;
raies = raiesSpectre[npoints];
indices=Table[2*p-1,{p,1,20}];
precision=Table[Abs[coefFourier[[2]]]/Abs[coefFourier[[2*p]]]/(2*p-1),{p,1,20}];
            
n|S1|/(n|Sn|)
11
30.999947
50.999842
70.999684
90.999474
110.999211
130.998895
150.998527
170.998106
190.997633
210.997107
230.99653
250.995899
270.995217
290.994482
310.993695
330.992857
350.991966
370.991023
390.990029

The increase in the number of points is costly in calculation, especially if one seeks harmonics of high order. Another solution is to apply low pass filtering to the signal before calculating the DFT. Indeed, let us consider a low pass filter which eliminates, or strongly attenuates, the harmonics starting from the order n = 40. A good precision on the calculation of the spectrum for n <40 can be obtained with only N = 80. The filter used must have a response that is as constant as possible in the passband, so as not to alter the harmonics to be calculated.

It may interest you

Leave a Reply

Your email address will not be published.

Solve : *
30 + 24 =