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'}]
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|'}]]
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|'}]]
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|'}]]
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|'}]]
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|) |
| 1 | 1 |
| 3 | 0.998684 |
| 5 | 0.996057 |
| 7 | 0.992122 |
| 9 | 0.986892 |
| 11 | 0.980376 |
| 13 | 0.972592 |
| 15 | 0.963556 |
| 17 | 0.953292 |
| 19 | 0.941822 |
| 21 | 0.929174 |
| 23 | 0.915377 |
| 25 | 0.900464 |
| 27 | 0.884471 |
| 29 | 0.867433 |
| 31 | 0.849391 |
| 33 | 0.830387 |
| 35 | 0.810465 |
| 37 | 0.789671 |
| 39 | 0.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|) |
| 1 | 1 |
| 3 | 0.999947 |
| 5 | 0.999842 |
| 7 | 0.999684 |
| 9 | 0.999474 |
| 11 | 0.999211 |
| 13 | 0.998895 |
| 15 | 0.998527 |
| 17 | 0.998106 |
| 19 | 0.997633 |
| 21 | 0.997107 |
| 23 | 0.99653 |
| 25 | 0.995899 |
| 27 | 0.995217 |
| 29 | 0.994482 |
| 31 | 0.993695 |
| 33 | 0.992857 |
| 35 | 0.991966 |
| 37 | 0.991023 |
| 39 | 0.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
- Telemecanic programmable logic controllers
- Best site for healthy recipes, cooking and beauty for women
- Download Vijeo Citect Version 7.5 software
- Download proteus software version 8
- Download the student version of the EPLAN Electric software