Download Fourier Transformation - Module and more Assignments Mathematics in PDF only on Docsity!
The Fourier Transform
The Fourier transform is crucial to any discussion of time series analysis, and this chapter discusses the
definition of the transform and begins introducing some of the ways it is useful.
We will use a Mathematica -esque notation. This includes using the symbol I for the square root of
minus one. Also, what is conventionally written as sin(t) in Mathematica is Sin[t] ; similarly the cosine
is Cos[t]. Finally, the irrational number 2.71828... is represented by the symbol E.
The contents of this chapter are:
Fourier Series
Fourier Transform
Example and Interpretation
Oddness and Evenness
The Convolution Theorem
Discrete Fourier Transforms
Definitions
Example
Implementation
Author
√ Fourier Series
Recall the Fourier series, in which a function f[t] is written as a sum of sine and cosine terms:
f#t'
a
0
cccccc
n 1
a
n
Cos#nt' ≈
n 1
b
n
Sin#nt'
or equivalently:
f#t' ≈
n
c
n
E
Int
n
c
n
+Cos#nt' I Sin#nt'/
The coefficients are found from the fact that the sine and cosine terms are orthogonal, from which:
a
n
cccc
S
t 0
2 S
f#t'Cos#nt'≈ t
b
n
cccc
S
t 0
2 S
f#t'Sin#nt'≈ t
Fourier series are used, for example, to discuss the harmonic structure of the tonic and overtones of a
vibrating string.
Note that the series represents either f[t] over a limited range of 0 < t < 2S, or we assume that the
function is periodic with a period equal to 2S.
Also note that, as opposed to the Taylor series, the Fourier series can represent a discontinuous func-
tion:
S S 2 S 3 S
t
1
f+t/
f#t'
cccc
cccc
S
Sin#t'
cccccccccccccccccc
Sin# 3 t'
cccccccccccccccccccccc
Sin# 5 t'
cccccccccccccccccccccc
√ Fourier Transform
In the previous section we defined the series over the interval (0, 2S). Say instead we are interested in
the interval (-L, L). Then the coefficients in the Fourier series are:
a
n
cccc
L
L
L
f#t' Cos$
nSt
cccccccccc
L
(≈ t
b
n
cccc
L
L
L
f#t' Sin$
nSt
cccccccccc
L
(≈ t
Thus we write the series of f as a function of a dummy variable x as:
f#x'
cccc
S
0
≈ Z
≈ t f#t'Cos#Z+t x/'
cccccccc
2 S
≈ Z
≈ t f#t' Cos#Z+t x/'
Since the sine is odd:
Sin#T' Sin#T'
we can write:
0 I
L
N
M
M
cccccccc
2 S
≈ Z
≈ t f#t'Sin#Z+t x/'
\
^
]
]
Adding this to the expression for f[x] gives:
f#x'
cccccccc
2 S
≈ Z e
IZx
t f
t '
e
IZt
The Fourier transform F 1
[Z ] of f[t] is:
F
1
#Z'
f#t' e
IZt
≈ t
Note that it is a function of Z. If we interpret t as the time, then Z is the angular frequency. Thus we
have replaced a function of time with a spectrum in frequency.
The inverse Fourier transform takes F[ Z ] and, as we have just proved, reproduces f[t]:
f#t'
cccccccc
2 S
F
1
#Z' e
I
Z t
≈ Z
You should be aware that there are other common conventions for the Fourier transform (which is why
we labelled the above transforms with a subscript). For example, some texts use a different normalisa-
tion:
F
2
#Z'
cccccccccccccc
r
2 S
f#t' e
IZt
≈ t
f#t'
cccccccccccccc
r
2 S
F
2
#Z' e
IZt
≈ t
Still others reverse the transform and its inverse:
F
3
#Z'
cccccccccccccc
r
2 S
f#t' e
I
Z t
≈ t
f#t'
cccccccccccccc
r
2 S
F
3
#Z' e
IZt
≈ t
The only difference between the "type-2" definition and the "type-3" one is the relative signs of the real
and imaginary parts of the transforms.
By default, Mathematica uses this "type-3" definition of the Fourier transform. In this class we will
almost always be using the "type-1" convention.
Say we have a function of the position x : g[x]. Then the type-1 Fourier transform and inverse transform
are:
G
1
#k'
g#x' e
Ikx
≈ x
and:
g#x'
cccccccc
2 S
G
1
#k' e
Ikx
≈ k
In this case the transform is a function of the wavenumber k = 2 S / O.
√ Example and Interpretation
Say we have a function:
I F#Z'
cccccccc
2 S
Sin#+Z
0
Z/nS s Z
0
cccccccccccccccccccccccccccccccccccccccccccccccccccc
2 +Z
0
Z/
For n = 5 and Z 0
= 100, the right hand side of the above looks like:
60 80 100 120 140 160
The zeroes in the above occur at:
Z
0
Z
cccccccccccccccc
Z
0
n
GZ
ccccccc
Z
0
n î 1, î 2, î 3, ...
Since the contributions outside the central maximum are small, we may take:
GZ
Z
0
cccccc
n
to be a measure of the width of the peak.
For n = 10 and the same value of Z 0
the plot looks like:
60 80 100 120 140 160
Clearly, the width of the curve is now decreased.
Curves such as the above will occur sufficiently often that we will give the function that generates them
a name: the sinc :
Sinc#x' ù
Sin#Sx'
ccccccccccccccccccccc
Sx
One interpretation of the above Fourier transform is that F[ Z ] is the frequency spectrum of a sine wave
signal f[t] which is varying in time; thus Z is the angular frequency. The main frequency component
occurs at the frequency of the sine wave, Z 0
, but there are other frequency components that cancel out
the signal for values of the time whose absolute value is greater than n S / Z 0
. If we think about letting n
go to infinity, then the sine wave is non-zero for all values of the time from -à to à; in this case the
width of the Fourier transform goes to zero and become a Dirac delta function centered at Z 0
.
Thus if we have an infinite sine wave but only measure it for a finite period of time, the measurement
will introduce "sidebands" in the frequency spectrum.
Another interpretation of the transform is that the symbol t is the finite width of a slit; the Fourier
transform of f[t] is then the amplitude of the diffraction pattern of the slit. The fact that a wider slit
produces a narrower transform means that to get, say, good dispersion of the high tones from a loud-
speaker requires that the speaker be small.
Yet another interpretation is that f[t] is the amplitude of an electromagnetic wave that is passing by us.
The period of the sine wave itself is
T
2 S
cccccccc
Z
0
and there are n cycles of the sine wave in f[t] , so it takes a time:
't n
2 S
cccccccc
Z
0
for the wave to pass us.
The width in the peak of the Fourier transform is a way of saying there is an uncertainty in the "true"
value of the frequency. But we know from Planck that the frequency is related to the energy E accord-
ing to:
E
h
cccccccc
2 S
Z
Thus, the uncertainty in the frequency corresponds to an uncertainty in the energy:
'E
h
cccccccc
2 S
GZ
Above we estimated GZ to be Z 0
/ n so:
'E
h Z
0
cccccccccccc
2 S n
√ The Convolution Theorem
I hope that after going through some of the interpretations of the Fourier transform above, you are
already convinced that it is one of the "keys to the universe." Here we present one of the most impor-
tant keys in the context of time series analysis.
Imagine we have a function f[t] whose Fourier transform is F[ Z ] , and another function g[t] whose
transform is G[ Z ]. Then the convolution is:
f#t' g#t'
√
f#u' g#t u'≈ u
We write g[t - u] in terms of the inverse Fourier transform:
g#t u '
cccccccc
2 S
G#Z'E
IZ +t u/
≈ Z
Thus:
f#t' g#t'
√
f#u'
cccccccc
2 S
G#Z'E
I
Z + t
u
/
≈ Z ≈ u
cccccccc
2 S
G#Z'E
IZt
f#u'E
IZu
≈ u≈ Z
But the right hand integral above is just the Fourier transform of f[u] , so:
f#t' g#t'
cccccccc
2 S
F#Z' G#Z'E
IZt
≈ Z
In words:
The inverse Fourier transform of a product of Fourier transforms is the convolution of the original
functions.
Here is an example. The Fourier transform of a Sinc function is just the rectangle function that in the
Convolution chapter we gave the symbol æ:
t
f#t' Sin#Z 0
t' s +Z 0
t/
Z 0
Z 0
Z
F#Z' ¾
Thus F[ Z ] only passes frequencies with an absolute value less than Z 0
. So if we convolve a Sinc func-
tion with some signal g[t] , the Sinc is performing as an ideal "low pass filter."
The example shows that we can consider the design of a filter in two separate domains:
- In the frequency domain we can just multiply the Fourier transform of the original time series by
some desired filter, a low pass filter in the above example.
- In the time domain, we can convolve the time series with the inverse Fourier transform of the desired
filter.
The fact that the amplitude of the Sinc function approaches zero asymptotically as t ë ±à means that
doing the full convolution would require measuring f[t] from t = -à to t = +à. Since this is physically
impossible, we have proved that the ideal low pass filter can not be built.
√ Discrete Fourier Transforms
So far in this chapter we have only considered continuous functions f[t]. Here we extend to a time
series that is a sample of f[t]. This section is divided into three subsections: Definitions, Examples ,
and Implementation.
Ç Definitions
If f[t] is a times series of length n , then we replace the continous Fourier transform:
F#Z'
f#t' e
I
Zt
≈ t
with a sum:
F#Z
j
L
N
M
M
M
M
k 0
n 1
f#t
k
' E
I Z j
t k
\
^
]
]
]
]
So now we can write the discrete Fourier transform as:
F#Z
j
L
N
M
M
M
M
k 0
n 1
f
k
E
I j +
2 S
cccccc
T
/ k'
\
^
]
]
]
]
L
N
M
M
M
M
k 0
n 1
f
k
E
I 2 S j k s n
\
^
]
]
]
]
To get the inverse Fourier transform:
f
t
k
cccccccc
S
F
#Z'
E
IZt
w
L
N
M
M
M
M
M
cccccccc
S
j 0
n 1
F
j
E
I Z
j
t
k
\
^
]
]
]
]
]
GZ
The value of GZ is just how much Z j
changes with each change from j to j + 1. We just saw that it is:
GZ
2 S
cccccccc
T
2 S
cccccccc
n'
So the discrete inverse Fourier transform is:
f#t
k
L
N
M
M
M
M
M
cccc
n
j
n
1
F
j
E
I 2
S j k
s n
\
^
]
]
]
]
]
cccc
Note that in both the transform and its inverse, by making the usual choice that the sampling interval '
is one simplifies the definition. We shall make that choice for the remainder of this chapter. Thus:
F
j
k 0
n
1
f
k
E
I 2
S j k
s n
f
k
cccc
n
j 0
n
1
F
j
E
I 2
S j k
s n
We recall that F j
is a shorthand for F[ Z
j
], where:
Z
j
j-
2 S
cccccccc
n'
1 j-
2 S
cccccccc
n
and that f k
is a shorthand for f[ t
k
] , where:
t
k
k' k
Ç Example
As already mentioned, the built in discrete Fourier transform routines in Mathematica use a different
convention than we use by default in these notes. Earlier, we called Mathematica 's choice a "type-3"
definition, so we write:
F
3
##j''
cccccccccc
r
n
k 1
n
f##k'' E
I 2 S +j 1 / + k 1 / s n
Here the sign of the exponential is different than the type-1 definition. Also note that k goes from 1
through n , because that is the way the Mathematica accesses members of a list. Recall that we usually
write the elements of a time series as:
timeSeries f
0
, f
1
, f
2
, ... , f
n
1
Then Mathematica 's access of elements of the list is:
timeSeries##k'' f
k 1
The Mathematica function implementing this definition is called Fourier.
The inverse Fourier transform in Mathematica , InverseFourier , is:
f
3
##k''
cccccccccc
r
n
j 1
n
F
3
##j''E
I 2
S + j
1
/ + k
1
/ s n
It is fairly simple to use Mathematica 's functions to implement the "type-1" conventions that we have
been using, but in this subsection we will not bother. The only real difference between the two conven-
tions is the signs of the real and imaginary parts and the normalisation.
Say we have a n = 200 element time series of five cycles of a sine wave:
convince yourself that the frequency can just as well be a negative number as a positive one. The alias
is a representation of the negative frequency component of the Fourier transform. Ideally we would like
it to look something like:
Z
Ideal plot of the Fourier transform
Note : the above plot also shows the same bug as the previous one. There is an extra dot high on the
vertical axis. Ignore it.
But since the Fourier transform is just a list of numbers, not {frequency, number} pairs, Mathematica
couldn't do that. Instead it joined the negative frequency values to the end of the positive frequency
ones.
We will have a lot more to say about aliasing in the next chapter.
Ç Implementation
One "brute force" way of calculating the sums of a Fourier transform is to define:
W ù E
I2Ssn
Then:
F
j
k 0
n
1
f
k
E
I 2
S j k
s n
k 0
n
1
f
k
W
jk
We can think of f as a vector of length n, and W as a matrix of dimension n l k. This multiplication
requires n
2
calculations, and evaluating the sum requires a smaller number of operations to generate the
powers of W. Thus calculating the Fourier transform this way is a O( n
2
) process. Doubling the number
of points in a time sample quadruples the time necessary to calculate the transform; tripling the number
of points requires nine times as much time.
However, there is a Fast Fourier Transform algorithm that can give an immense improvement. The
basic idea is that we split the sum into two parts:
F
j
k 0
n
s 2
1
f
2 k
E
I 2
S j 2 k
s n
k 0
ns 2 1
f
2 k 1
E
I 2 S j +2 k 1 /s n
The first sum involves the even terms in f , and the second sum the odd ones.
Using the same W as before, we can write:
F
j
k 0
ns 2 1
f
2 k
E
I 2 S j k s +ns 2 /
W
k
k 0
ns 2 1
f
2 k 1
E
I 2 S j k s +ns 2 /
F
k
e
W
k
F
k
o
We can apply the same procedure recursively to the sums represented by F
k
e
and F
k
o
. Eventually, if n is
a power of two, we end up with no summations at all, just a product of terms.
It turns out that this procedure is O(n log
2
n). So, doubling the number of points in the time series only
doubles the time necessary to calculate the transform using this algorithm, and tripling the number of
points increases the time by about 4.75.
Many people treat this Fast Fourier Transform as "magic", and it is does seem magical in its properties.
But notice that it all hinges on n being a power of two. For example, if one constructs a time series:
timeSeries ≈
t 1
32700
E
ts+ 32700 s 6 /
the values range from 0.999817 to 0.00247. Mathematica 's Fourier in one computing environment
took 13.12 seconds to transform this series. Since the last terms in the series are so close to zero, it is
quite reasonable to add zeroes to the end of the time series. We added 68 zeroes, so the total length of
the series becomes 32,768 = 2
15
. Now it took Fourier only 5.16 seconds to do the transform.
The lesson, then, is that if one is taking data in the lab and later the data will be Fourier transformed, if
there is a lot of data the total number of data points should always be a factor of two.
