The complex mathematical formules inside this text are there for reference. They are not essential to understand the text."
What is a transformation into a Fourier series ?
Suppose an audio signal sampled at 8kHz. This means that every succesive eighth of a milliseconds one makes a measurement of the intensity of the signal.
For the remainder of this text, we will assume we're working on a sample of EIGHT measurements.
Here is an exemple of such a sample :
50 mV 206 mV -100 mV -65 mV -50 mV -6 mV 100 mV -135 mV
We could represent it by this graph :
The intersection of the first vertical bar with the drawing of the signal is at 50 mV, the second at 206 mV, etc...
Our goal is to break up this signal into a sum of sinusoids. In the case of our example, that yields :
(If you sum these sinusoids, you get our signal as a result.)
Or, written with words :
- A sinusoid of 1 kHz with an amplitude of 50 mV, and a phase shift of pi/2. (or a time shift of 2 eight-thousandth of a second)
- A sinusoid of 2 kHz with an amplitude of 100 mV, and a phase shift of 0.
- A sinusoid of 3 kHz with an amplitude of 100 mV, and a phase shift of 0.
S(t) = 50mV . sin (2 pi 1000 t + pi/2) +
100mV . sin (2 pi 2000 t + 0 ) +
100mV . sin (2 pi 3000 t + 0 )
or :
S(t) = 50mV . sin (2 pi 1000 (t + 2/8000)) +
100mV . sin (2 pi 2000 (t + 0 )) +
100mV . sin (2 pi 3000 (t + 0 ))
The way this decomposition is carried out is traditionally the following one :
One DECIDES the signal is composed of sinusoids of 0 kHz, 1 kHz, 2 kHz, 3 kHz, and 4 kHz.
Then remains to determine the amplitude and the phase shift of each sinusoid.
(If a sinusoid is not at all present in the signal, its amplitude will be null, and its phase shift will be of no importance.)
(The sinusoid of 0 kHz is in fact simply the DC component of the signal, a constant. Its phase does not have any importance, as long it is not infinite.)
(Determining the phase and the amplitude of the sinusoid at 4 kHz is somehow wobbly : several sinusoids, with different amplitudes and phase shifts, can do the job and yield a correct result.)
(One considers the signal to be repetitive : it is supposed that a ninth measure would have yielded the same result as the first, a tenth the same as the second...)
How to do it
Our goal is thus to determine the amplitude and the phase shift of the sinusoids. To understand how it works, the reader must considerate and accept the three following mathematical facts :
1. A signal can be represented as being the sum of several pikes.
If we state the equation of a pike is the following :
y(t) = pike (a, dc, t)
with a the amplitude, dc the time shift, and t the time (the X-coordinate).
Then our S(t) signal can be written the following way :
S(t) = pike ( 50mV, 0, t) +
pike ( 206mV, 1/8000, t) +
pike (-100mV, 2/8000, t) +
pike ( -65mV, 3/8000, t) +
pike ( -50mV, 4/8000, t) +
pike ( -6mV, 5/8000, t) +
pike ( 100mV, 6/8000, t) +
pike (-135mV, 7/8000, t)
2. A pike shaped signal can be broken up into a precise Fourier transform series.
Say Fe the sampling rate (8000 Hz in our case)
Say ne the number of measures (8 in our case)
(In human language : one produces a pike by making the sum of cossinusoids.)
The same, in complex exponential notation (see the appendix) :
The complete graph of a pike with amplitude 200 mV and shift 0 is :
The five sinusoids it is the sum of are :
pique(200mV, 0, t) = 200 . 1/8 . sin (2 pi 0 (t - 0) + pi/2) +
200 . 1/4 . sin (2 pi 1000 (t - 0) + pi/2) +
200 . 1/4 . sin (2 pi 2000 (t - 0) + pi/2) +
200 . 1/4 . sin (2 pi 3000 (t - 0) + pi/2) +
200 . 1/8 . sin (2 pi 4000 (t - 0) + pi/2)
That is to say :
25 . sin (2 pi 0 t + pi/2) + 50 . sin (2 pi 1000 t + pi/2) + 50 . sin (2 pi 2000 t + pi/2) + 50 . sin (2 pi 3000 t + pi/2) + 25 . sin (2 pi 4000 t + pi/2)
And with complex exponentials :
pique(200mV, 0, t) = RE( 200 . 1/8 . 1 . e 0,2 pi 0 t +
200 . 1/4 . 1 . e 0,2 pi 1000 t +
200 . 1/4 . 1 . e 0,2 pi 2000 t +
200 . 1/4 . 1 . e 0,2 pi 3000 t +
200 . 1/8 . 1 . e 0,2 pi 4000 t )
That is to say :
RE( 25 e 0,2 pi 0 t +
50 e 0,2 pi 1000 t +
50 e 0,2 pi 2000 t +
50 e 0,2 pi 3000 t +
25 e 0,2 pi 4000 t )
You will notice this graph is not null everywhere away from the pike. What matters is the fact it is null at the exact places where the other measures are situated.
(Note that this graph presents another pike, identical, in the series of measures that follows.)
3. The sum of two sinusoids having the same frequencies (but with different phase shifts and amplitudes), is a sinusoid.
(sgn (x) yields 1 if x is positive, and -1 if x is negative.)
Using complex exponential notation :
Synthesis
If the signal is the sum of eight pikes,
S(t) = pique (a1, dc1, t) +
pique (a2, dc2, t) +
pique (a3, dc3, t) +
pique (a4, dc4, t) +
pique (a5, dc5, t) +
pique (a6, dc6, t) +
pique (a7, dc7, t) +
pique (a8, dc8, t)
and if each pike is the sum of five sinusoids,
then the signal is the sum of 8 . 5 = 40 sinusoids.
However, all the pikes are made up of five sinusoids having the same frequencies (0 kHz, 1 kHz, 2 kHz, 3 kHz, 4 kHz). (but with different amplitudes and phase shifts)
In other words; amongst the 40 sinusoids, there are 8 sinusoids of 0 kHz, 8 sinusoids of 1 kHz, 8 sinusoids of 2 kHz, 8 sinusoids of 3 kHz, and 8 sinusoids of 4 kHz.
The sum of 8 sinusoids having the same frequency, yields one sinusoid, having that same frequency.
One can thus simplify the sum of the 40 sinusoids, downto a sum of just 5 sinusoids :
S(t) = a0 sin (2 pi 0 t + d0) +
a1 sin (2 pi 1000 t + d1) +
a2 sin (2 pi 2000 t + d2) +
a3 sin (2 pi 3000 t + d3) +
a4 sin (2 pi 4000 t + d4)
Using the complex notation :
S(t) = RE( (a0, b0) e 0,2 pi 0 t +
(a1, b1) e 0,2 pi 1000 t +
(a2, b2) e 0,2 pi 2000 t +
(a3, b3) e 0,2 pi 3000 t +
(a4, b4) e 0,2 pi 4000 t )
Our sample of eight measures has thus be broken down, a mechanical and imparable way, into the sum of five sinusoids.
(The initial signal is made out of 8 numbers (50, 206, -100, -65, -50, -6, 100, -135). One breaks it up into 5 sinusoids, characterized each one by 2 numbers (amplitude, phase shift), that is to say 10 numbers. Does the result of a transformation into Fourier series thus use more place than the initial signal ? The answer is no : one can drop away the phase shifts of the sinusoids of 0 kHz and 4 kHz. That way 8 numbers remain.)
Application
Here a program written in BASIC which implements the pikes system to carry out a decomposition in Fourier series.
As far as possible, the author tried to use variables names identical to those contained in the text, and wrote the program so it could run on the majority of BASIC interpreters or compilers (it has been tested on BBC BASIC and QBASIC (this later is provided as a standard with MS-DOS, thus available on old PC computers)). (This program uses the arctangent function (ATN). If you wish a version that uses the arccosine function (ACOS), clic here.)
If the following parameters are given to the program :
SAMPLING RATE ? 8000 NUMBER OF SAMPLES ? 8
Then the eight following numbers :
50 206 -100 -65 -50 -6 100 -135
The program will answer this, with always some little round-off errors :
0Hz A=0 D=1.57079633 1000Hz A=50 D=1.57079633 2000HZ A=100 D=0 3000HZ A=100 D=0 4000HZ A=0 D=3.14159265
Which commes over well with the example given at the beginning of this text. (Phase shifts do not have any importance for the sinusoids whose amplitude or frequency is null.)
The program :
PY = 4 * ATN(1) : REM 3.14159265, YOU KNOW
INPUT "SAMPLING RATE ", FE
INPUT "AMOUNT OF MEASURES ", NE
NS = NE / 2 + 1: REM AMOUNT OF SINUSOIDS
REM CREATE A MATRIX C WITH (COS(DN) , -SIN(DN)) :
DIM C(NE, NS, 2)
FOR N = 0 TO NE / 2
FOR E = 0 TO NE - 1
DC = E / FE
FEN = N * FE / NE
DN = 2 * PY * FEN * DC
C(E, N, 1) = COS(DN)
C(E, N, 2) = -SIN(DN)
NEXT E
NEXT N
REM FILL THE MATRIX B WITH
REM THE AMPLITUDES :
DIM B(NS)
B(0) = 1 / NE
FOR N = 1 TO NE / 2 - 1
B(N) = 2 / NE
NEXT N
B(NE / 2) = 1 / NE
REM S IS THE MATRIX CONTAINING THE MEASURES :
DIM A(NE)
REM R IS THE MATRIX CONTAINING THE RESULTS :
DIM R(NS, 2)
REM CONSTANTS HAVE BEEN CALCULATED
REM HERE IS THE "LIVING" PART :
REM ENTER THE MEASURES :
FOR E = 0 TO NE - 1
PRINT "MEASURE "; E; " ";
INPUT A(E)
NEXT E
REM THE CALCULATION :
FOR N = 0 TO NE / 2
R(N, 1) = 0
R(N, 2) = 0
FOR E = 0 TO NE - 1
A = A(E) * B(N) * C(E, N, 1)
B = A(E) * B(N) * C(E, N, 2)
R(N, 1) = R(N, 1) + A
R(N, 2) = R(N, 2) + B
NEXT E
NEXT N
REM OUTPUT THE RESULTS :
PRINT "COMPLEXS COEFFICIENTS :"
FOR N = 0 TO NE / 2
PRINT N * FE / NE; "HZ ("; R(N, 1); ","; R(N, 2); ")"
NEXT N
PRINT "AMPLITUDE AND PHASE SHIFTS OF THE SINUSOIDS :"
FOR N = 0 TO NE / 2
A = SQR(R(N, 1) ^ 2 + R(N, 2) ^ 2)
IF -R(N, 2) = 0 THEN
D = PY / 2
IF R(N, 1) < 0 THEN D = -D
ELSE
D = ATN(R(N, 1) / -R(N, 2))
REM NO ARCCOSSINUS IN QBASIC
IF -R(N, 2) < 0 THEN D = D + PY
ENDIF
PRINT N * FE / NE; "HZ A="; A; " D="; D
NEXT N
Appendix :
Reminder on exponentials and complex numbers
Thanks
I wish to thank my friends Frédéric Cloth and Didier Bizzarri for their help and data.
Eric Brasseur
-
1994
[ En français ]