( ) ( ) t , tN = N t = T

AMS 147 Computational Methods and Applications
Lecture 14
Copyright by Hongyun Wang, UCSC
Recap:
*) Application #2 of cubic spline:
Find the distance from a point to a trajectory represented by N data points
*) Least squares fitting
*) Fourier transform
Suppose ƒ(t) is periodic with period = T. It can be written as
f (t ) =
c
k = k
kt exp i2 T
where the Fourier coefficients {ck} has the expression
1
ck =
T
T
kt f (t ) exp i2 T dt
0
Fourier transform is a mapping
{ f (t ) ,
0 t < T} { ck ,
< k < }
Next, we discuss the discrete version of Fourier transform.
Discrete version of Fourier transform:
Question: Why do we discuss discrete version of Fourier transform?
Answer:
•
Data is discrete
•
Computation is always done in discrete form
Consider a grid.
t =
T
N
t j = j t ,
j = 0, 1, …, ( N 1)
t 0 = 0 , t1 = t , … , t N 1 = ( N 1) t , t N = N t = T
{}
(Draw the real axis to show the grid t j )
-1-
AMS 147 Computational Methods and Applications
{t
The grid:
j
}
, 0 j N 1
{
}
Suppose function f ( t) is known only on the grid t j , 0 j N 1 .
That is, ƒ(t) is represented only by N data points:
{ f , j = 0, 1, …, N 1}
j
where
( )
fj = f t j
The exact Fourier coefficients require the full information about function f(t).
1
ck =
T
T
kt
f (t ) exp i 2 T dt
0
{ f , j = 0, 1, …, N 1} , we use the composite trapezoidal rule to approximate
With only the data
j
the Fourier coefficients.
The composite trapezoidal rule applied to a periodic function:
g0 + gN N 1 g
t
dt
t
+ gj (
)
0
2
j =1
T
Function g(t) is periodic
==>
==>
gN = g (T ) = g ( 0 ) = g0
T
N 1
0
j=0
g (t ) dt t g j
Applying the composite trapezoidal rule to approximate the Fourier coefficients, we have
ck kt j 1 N 1
t f j exp i 2
T
T j=0
1 T
1
1
t = =
T N N
T
Using
==>
1
ck N
N 1
f
j=0
j
and
k T kj
k
tj = j =
T N N
T
k j
exp i 2 N
Let
c˜k =
1
N
N 1
f
j =0
j
k j
exp i 2 N
-2-
AMS 147 Computational Methods and Applications
{ck }
are called discrete Fourier coefficients.
It appears that the discrete version of Fourier transform is a mapping
{ f , j = 0, 1, …, N 1}
j
{ck ,
< k < }
As we will see later, the discrete Fourier coefficients satisfy
ck + N = ck
and thus, the discrete version of Fourier transform is a mapping
{ f , j = 0, 1, …, N 1}
j
N
N
ck , < k 2
2
We study the discrete Fourier coefficients ck and their connection to the exact Fourier
coefficients ck.
kt
We start by examining the behavior of Fourier mode exp i2 on the grid t j
T
{}
kt
Fourier mode exp i2 has k oscillations in [0, T ] .
T
{}
The grid t j has N points in [0, T ]
==>
Each oscillation is represented by
{}
N
points on the grid t j .
k
We need at least 2 points to resolve each oscillation.
{}
The grid t j can only possibly resolve
==>
N
kt exp i 2 , for k 2
T
For k <<
{}
kt
N
, Fourier mode exp i2 is well resolved on the grid t j and we have
T
2
ck ck for k <<
For k >
N
2
{}
kt
N
, the grid t j cannot resolve Fourier mode exp i2 .
T
2
kt
N
Question: What happens to exp i2 and ck when k > ?
T
2
Answer:
Aliasing
-3-
AMS 147 Computational Methods and Applications
Aliasing:
( k + N )t = exp i 2 kt .
On the grid t j , exp i 2
T
T t t j j
{}
Claim:
Verify:
(k + N ) t exp i2
T
Using
t =t j
N tj
T
=
k tj
N tj
= exp i2
i2
T
T N T
j = j
T N
ktj
= exp i 2
i 2 j
T
ktj = exp i 2
exp ( i 2 j )
T k t
= exp i 2 T ==>
t =t j
(k + N ) t exp i2
T
t =t j
kt
= exp i2 T
t= t j
( k + N ) t and exp i 2 k t appear
That is, on the grid t j , Fourier modes exp i 2
T
T exactly the same. This is called aliasing.
{}
Therefore, we have
ck + N = ck
and
ck + n N = ck
for all integer values of n.
Summary of the aliasing:
{}
The grid t j can only resolve
kt exp i 2 ,
T
{}
On the grid t j , for k >
N
N
<k
2
2
kt
N
, the Fourier mode exp i2 appears as
T
2
-4-
AMS 147 Computational Methods and Applications
(k mN )t exp i2
T
where m is an integer satisfying N
N
< k mN . And we have
2
2
ck = ck m N
N
N
So we only need to keep discrete Fourier coefficients ck , < k .
2
2
A short digression (frequency, sampling rate and Nyquist frequency)
Definition:
Frequency =
# of waves ( # of oscillations )
time
kt k
Frequency of exp i 2 =
T T
Definition:
Sampling rate =
# of data points N
=
time
T
Definition:
Nyquist frequency = Largest frequency that can be resolved
=
N
2 = N
2T
T
Note that
1
Nyquist frequency= Sampling rate
2
Aliasing (in the language of frequency, sampling rate and Nyquist frequency):
With sampling rate =
N
k
, a mode of frequency appears as a mode of frequency
T
T
N
k
m
T
T
where m is an integer satisfying
-5-
AMS 147 Computational Methods and Applications
N
k
<
T
2T
Nyquist
frequency
N
N m T
2T
Sampling Nyquist
rate
frequency
Example:
T = 1,
N = 256
N
= 256 Hz
T
N
Nyquist frequency =
= 128 Hz
2T
Sampling rate =
A mode of frequency 60 Hz appears as a mode of frequency 60 Hz.
A mode of frequency 160 Hz appears as a mode of frequency 160 256 = 96 Hz .
A mode of frequency 260 Hz appears as a mode of frequency 260 256 = 4 Hz .
A mode of frequency 460 Hz appears as a mode of frequency 460 2 256 = 52 Hz .
Back to …
The discrete version of Fourier transform
{ f , j = 0, 1, …, N 1} j
c˜k =
1
N
N 1
f
j =0
j
N
N
ck , < k 2
2
k j
exp i 2 N
Note: the index range for ck is selected based on
ck ck , for N
N
<k
2
2
However, this index range is not the most convenient one for computation.
For the convenience of computation, we use a different index range, which leads us to discrete
Fourier transform.
DFT (discrete Fourier transform)
{
}
DFT: F = f j , j = 0,1,…, N 1 Y = {y k , k = 0,1,…, N 1}
N 1
kj
y k = f j exp i 2 , k = 0,1,…, N 1
N
j =0
Note: In Matlab, DFT is implemented by the built-in function fft.
-6-
AMS 147 Computational Methods and Applications
N
N
We study the connection between {y k , k = 0,1,…, N 1} and c˜k , < k .
2
2
{y k , k = 0,1,…, N 1} is what we get (the output of DFT).
N
N
c˜k , < k is what we want (it contains approximations to the true Fourier
2
2
coefficients)
We need to relate these two to each other.
N
2
y k = N c˜k
For 0 k N
< k N 1
2
For
yk = N ck = N ck N
==>
y0 , y1 , …, y N , y N +1 , …, yN 2 , yN 1 = N c0 , c1 , …, c N , c N +1 , …, c2 , c1 2
2
2
2
If we put the indices into a vector, we have
N
N
ind = 0, 1, …, , + 1, …, 2, 1
2
2
In Matlab,
ind = [0:N/2,-N/2+1:-1];
This vector is very useful in Matlab programming involving FFT.
IDFT (inverse discrete Fourier transform)
{
}
IDFT: Y = {y k , k = 0,1,…, N 1} F˜ = ˜f j , j = 0,1,…, N 1
˜f = 1
j
N
N 1
y
k=0
k
kj
exp i2 ,
N
j = 0,1,…, N 1
Note: In Matlab, IDFT is implemented by the built-in function ifft.
To justify the name “inverse discrete Fourier transform”, we need the theorem below.
DFT :
IDFT :
{f }
j
{ yk }
{ yk }
{ f }
j
-7-
AMS 147 Computational Methods and Applications
{f }
IDFT DFT :
j
{ f }
j
Theorem:
IDFT DFT = Identity
That is, fj = f j , for j = 0, 1, …, N 1
(See Appendix for the proof of this theorem)
p
When N = 2 , both DFT and IDFT can be implemented very efficiently.
We are going to discuss FFT (fast Fourier transform) after we discuss applications of FFT.
(Go through sample codes in assignment #4)
Appendix:
Proof of the theorem
Theorem:
IDFT DFT = Identity
That is, fj = f j , for j = 0, 1, …, N 1
Proof of the theorem:
˜f = 1
l
N
=
=
=
N 1
y
1
N
1
N
1
N
k =0
k
kl
exp i 2 N
N 1 N 1
f
k = 0 j =0
j
N 1 N 1
f
k = 0 j =0
j
kj
kl
exp i 2 exp i2 N
N
k ( l j )
exp i 2
N N 1
k (l j) f
j exp i 2 N j=0
k =0
N 1
(E01)
For j = l , we have
N 1
expi 2
k=0
k (l j )
=N
N For j l , we have
-8-
AMS 147 Computational Methods and Applications
k ( l j ) N 1 k
= ,
N k = 0
N 1
exp i 2
k=0
(l j ) = exp i 2
N Using the sum of geometric series
N 1
= 1 ,
k=0
N 1
k
we obtain
N 1
expi 2
k=0
k (l j ) exp (i2 (l j )) 1
=
=0
N (
)
l
j
1
exp i2
N Substituting this result back into (E01), we obtain
˜f = 1
l
N
N 1
N 1
f exp i2
j =0
j
k=0
k ( l j )
= fl
N -9-

Download Report