58 Chapter 58 Fourier Analysis

58
Fourier Analysis
Kenneth Howell
University of Alabama in Huntsville
58.1
58.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58.2 The Function/Functional Theory . . . . . . . . . . . . . . . . . . .
58.3 The Discrete Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58.4 Relating the Functional and Discrete Theories . . . . . . .
58.5 The Fast Fourier Transform . . . . . . . . . . . . . . . . . . . . . . . .
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
58-1
58-2
58-8
58-12
58-17
58-21
Introduction
Fourier analysis has been employed with great success in a wide range of applications. The underlying theory is based on a small set of linear transforms on particular linear spaces. It may, in fact,
be best to refer to two parallel theories of Fourier analysis—the function/functional theory and the
discrete theory—according to whether these linear spaces are infinite or finite dimensional. The function/functional theory involves infinite dimensional spaces of functions on Rn and can be further divided
into two “subtheories”—one concerned with functions that are periodic on Rn (or can be treated as
periodic extensions of functions on finite subregions) and one concerned with more general functions
and functionals on Rn . This theory played the major role in applications up to half a century or so ago.
However, the tools from the discrete theory (involving vectors in C N instead of functions on Rn ) are
much more easily implemented on digital computers. This became of interest both because much of the
function/functional theory analysis can be closely approximated within the discrete theory and because
the discrete theory is a natural setting for functions known only by samplings. In addition, “fast” algorithms for computing the discrete transforms were developed, allowing discrete analysis to be done very
quickly even with very large data sets. For these reasons the discrete theory has become extremely important in modern applications, and its utility, in turn, has greatly extended the applicability of Fourier
analysis.
In the first part of this chapter, elements of the function/functional theory are presented and illustrated. For expediency, attention will be restricted to functions and functionals over subsets of R1 . A
corresponding review of the analogous elements of the discrete theory then follows, with a discussion of
the relations between the two theories following that. Finally, one of the fast algorithms is described and
the extent to which this algorithm improves speed and applicability is briefly discussed.
The development here is necessarily abbreviated and covers only a small fraction of the theory and
applications of Fourier analysis. The reader interested in more complete treatments of the subject is
encouraged to consult the references given throughout this chapter.
58-1
58-2
Handbook of Linear Algebra
58.2
The Function/Functional Theory
Definitions:
The following function spaces often arise in Fourier analysis. In each case I is a subinterval of R, and all
functions on I are assumed to be complex valued.
r C 0 (I), the normed linear space of bounded and continuous functions on I with the norm f =
sup{| f (x)| : x ∈ I}.
r L1 (I), the normed linear space of absolutely integrable functions on I with the norm f =
I
| f (x)| d x.
r L2 (I), the inner product space of square integrable functions on I with the inner product f, g =
I
f (x)g (x) d x.
If I is not specified, then I = R.
Let φ ∈ L1 (R). The Fourier transform F[φ] and the inverse Fourier transform F −1 [φ] of φ are the
functions
F[φ]|x =
∞
−∞
φ(y)e −i 2π xy d y
F −1 [φ]|x =
and
∞
−∞
φ(y)e i 2π xy d y .
Additionally, the terms “Fourier transform” and “inverse Fourier transform” can refer to the processes for
computing these functions from φ.
A function φ on R is said to be periodic if there is a positive value p, called a period for φ, such that
φ(x + p) = φ(x)
∀x ∈ R .
The smallest period, if it exists, is called the fundamental period.
The Fourier series for a suitably integrable, periodic function φ with period p > 0 is the infinite series
∞
c k e i 2π ωk x
k=−∞
where
ωk =
k
p
ck =
and
1
p
p
φ(y)e −i 2π ωk y d y.
0
The theory for Fourier transforms and series can be generalized so that the requirement of φ being
“suitably integrable” can be greatly relaxed (see [How01, Chap. 20, 30–34] or [Str94, Chap. 1–4]). Within
this generalized theory the delta function at a ∈ R, δa , is the functional limit
1
δa (x) = lim+ pulse (x − a)
→0 2
where
pulse (s ) =
1
if |s | ≤ 0
if |s | > That is, δa is the (generalized) function such that
∞
−∞
ψ(x)δa (x) d x = lim+
→0
∞
−∞
ψ(x)
whenever ψ is a function continuous at a.
An array of delta functions is any expression of the form
∞
k=−∞
c k δkx
1
pulse (x − a) d x
2
.
58-3
Fourier Analysis
where x > 0 is fixed (the spacing of the array) and c k ∈ C for each k ∈ Z. If the array is also periodic
with period p, then the corresponding index period is the positive integer N such that
p = Nx
φk+N = φk
and
∀ k ∈ Z.
The convolution φ ∗ ψ of a suitably integrable pair of functions φ and ψ on R is the function given by
φ ∗ ψ(x) =
∞
−∞
φ(x − y)ψ(y) d y .
Facts:
All the following facts except those with specific reference can be found in [How01] or [Str94].
1. Warning: Slight variations of the above integral formulas, e.g.,
1
√
2π
2.
3.
4.
5.
∞
−∞
φ(y)e −i xy d y
and
1
√
2π
∞
−∞
φ(y)e i xy d y,
are also often used to define, respectively, F[φ] and F −1 [φ] (or even F −1 [φ] and F[φ]). There is
little difference in the resulting “Fourier theories,” but the formulas resulting from using different
versions of the Fourier transforms do differ in details. Thus, when computing Fourier transforms
using different tables or software, it is important to take into account the specific integral formulas
on which the table or software is based.
Both F and F −1 are continuous, linear mappings from L1 (R) into C 0 (R). Moreover, if both φ and
F[φ] are absolutely integrable, then F −1 [F[φ]] = φ.
In the generalized theory (which will be assumed hereafter), F and F −1 are defined on a linear space
of (generalized) functions that contains all elements of L1 (R) and L2 (R), all Fourier transforms
of elements of L1 (R) and L2 (R), all piecewise continuous periodic functions, the delta functions,
and all periodic arrays of delta functions. F and F −1 are one-to-one linear mappings from this
space onto this space, and are inverses of each other.
Every nonconstant, piecewise continuous periodic function has a fundamental period, and every
period of such a function is an integral multiple of that fundamental period. Moreover, any two
Fourier series for a single periodic function computed using two different periods will be identical
after simplification (e.g.,
after eliminating zero-valued terms).
The Fourier series k∈Z c k e i 2π ωk x for a periodic function φ with period p can be written in
trigonometric form
c0 +
∞
[ak cos(2π ωk x) + bk sin(2π ωk x)]
k=1
where
ak = c k + c −k
2
=
p
and
bk = i c k − i c −k
2
=
p
p
φ(x) cos(2π ωk x) d x
0
p
φ(x) sin(2π ωk x) d x .
0
6. On I = (0, p) (or any other interval I of length p), the exponentials
e i 2π ωk x : ωk = k/ p, k ∈ Z
58-4
Handbook of Linear Algebra
form an orthogonal set in L2 (I), and the Fourier series k∈Z c k e i 2π ωk x for a periodic function φ
with period p is simply the expansion of φ with respect to this orthogonal set. That is,
ck =
φ(x), e i 2π ωk x e i 2π ωk x 2
∀ k ∈ Z.
Similar comments apply regarding the trigonometric form for the Fourier series and the set
{1, cos(2π ωk x), sin(2π ωk x) : ωk = k/ p, k ∈ N} .
7. A periodic function φ can be identified with its Fourier series
of criteria. In particular:
k∈Z c k e
i 2π ωk x
under a wide range
r If φ is smooth, then its Fourier series converges uniformly to φ.
r If φ is piecewise smooth, then its Fourier series converges pointwise to φ everywhere φ is con-
tinuous.
r If φ is square-integrable on (0, p), then
p
lim
(M,N)→(−∞,∞)
0
2
N
i 2π ωk x φ(x) −
c
e
k
d x = 0.
k=M
r Within the generalized theory,
∞
lim
(M,N)→(−∞,∞)
−∞
N
i 2π ωk x
ψ(x) φ(x) −
ck e
dx = 0
k=M
whenever ψ is a sufficiently smooth function vanishing sufficiently rapidly at infinity (e.g., any
Gaussian function).
8. [BH95, p. 186] Suppose φ is a periodic function with Fourier series k∈Z c k e i 2π ωk x . If φ is m-times
differentiable on R and φ (m) is piecewise continuous for some m ∈ N, then there is a finite constant
β such that
|c k | ≤
β
|k|m
∀ k ∈ Z.
9. For each a ∈ R and function φ continuous at a:
∞
−∞
φ(x)δa (x) d x = φ(a)
φδa = φ(a)δa .
and
10. For each a ∈ R:
r F[δ ]| = e −i 2πax
a x
r F[e i 2πax ]| = δ
x
a
and
and
F −1 [e −i 2πax ] = δa .
F −1 [δa ]|x = e i 2πax .
11. A function φ is periodic with period p and Fourier series
are the arrays
F[φ] =
∞
c k δkω
and
k∈Z c k e
F −1 [φ] =
k=−∞
i 2π ωk x
∞
if and only if its transforms
c −k δkω
k=−∞
with spacing ω = 1/ p.
12. If
φ=
∞
φk δkx
k=−∞
is a periodic array with spacing x, period p, and corresponding index period N, then:
58-5
Fourier Analysis
r The Fourier series for φ is given by
∞
n e i 2π ωn x
n=−∞
where ωn = n/ p = n/(Nx) and
N−1
1 φk e −i 2πkn/N .
Nx
n =
k=0
r The Fourier transforms of φ are also periodic arrays of delta functions with index period N. Both
have spacing ω = 1/ p and period P = 1/x. These transforms are given by
F[φ] =
∞
n δnω
F
and
−1
∞
[φ] =
n=−∞
−n δnω .
n=−∞
13. (Convolution identities) Let f , F , g , and G be functions with F = F[ f ] and G = F[g ]. Then,
provided the convolutions exist,
F[ f g ] = F ∗ G
F[ f ∗ g ] = F G .
and
Examples:
1. For α > 0,
F[pulseα ]|x =
∞
−∞
=
α
pulseα (y)e −i 2π xy d y
e −i 2π xy d y
−α
e i 2π αx − e −i 2π αx
sin(2π αx)
=
.
i 2π x
πx
=
2. Let f be the periodic function,
f (x) =
|x|
for − 1 ≤ x < 1
f (x − 2)
for all x ∈ R
The period of this function is p = 2, and its Fourier series
ωk =
k
2
k∈Z c k e
(so e i 2π ωk x = e i π kx )
and
1
ck =
2
=
2
0
|x|
if x < 1
|x − 2|
if 1 < x
−2
(kπ)
e −i πkx d x
1/2
if k = 0
.
k
(−1) − 1 if k = 0
.
i 2π ωk x
is given by
58-6
Handbook of Linear Algebra
Since f is continuous and piecewise smooth, f equals its Fourier series,
f (x) =
(−1)k − 1
1
e i π kx
+
2
k2π 2
∀x ∈ R.
k∈Z\{0}
Moreover,
F[ f ] =
(−1)k − 1
1
δ0 +
δk/2 .
2
k2π 2
k∈Z\{0}
3. Let
φ=
∞
φk δx
k=−∞
be the periodic array with spacing x = 1/3, index period N = 4, and coefficients
φ0 = 0, φ1 = 1, φ2 = 2,
and φ3 = 1 .
The period p of φ is then
p = Nx = 4 ·
1
4
= .
3
3
Its Fourier transform = F[φ] is also a periodic array,
∞
=
n δnω ,
n=−∞
with index period N = 4. The spacing ω and period P of are determined, respectively, from
the period p and spacing x of φ by
ω =
1
3
1
=
=
p
4/3
4
and
P =
1
1
=
= 3.
x
1/3
The coefficients are given by
n =
N−1
1 φk e −i 2πkn/N
Nx
k=0
1 −i 2π 0n/4
=
φ0 e
+ φ1 e −i 2π1n/4 + φ2 e −i 2π2n/4 + φ3 e −i 2π3n/4
Nx
1 0
=
0e + 1e −i nπ/2 + 2e −i nπ + 1e −i 3nπ/2
4(1/3)
=
3
[0 + (−i )n + 2 + i n ] .
4
Thus,
0 = 3,
1 =
3
,
2
2 = 0
and
3 =
3
.
2
58-7
Fourier Analysis
Applications:
1. [BC01, p. 155] or [Col88, pp. 159–160] (Partial differential equations) Using polar coordinates, the
steady-state temperature u at position (r, θ) on a uniform, insulated disk of radius 1 satisfies
r2
∂ 2u
∂ 2u
∂u
+ 2 =0
+r
2
∂r
∂r
∂θ
for 0 < r < 1 .
As a function of θ, u is periodic with period 2π and has (equivalent) Fourier series
∞
c k e i kx
c0 +
and
k=−∞
∞
[ak cos(kx) + bk sin(kx)]
k=1
where the coefficients are functions of r . If u satisfies the boundary condition
u(1, θ) = f (θ)
for 0 ≤ θ < 2π
for some function f on [0, 2π ), then the coefficients in the series can be determined, yielding
u(r, θ) =
∞
r k γk e i kθ = γ0 +
k=−∞
k=1
where
γk =
αk =
1
π
∞
k
r αk cos(kθ) + r k βk sin(kθ)
1
2π
2π
f (θ)e −i kθ dθ ,
0
2π
f (θ) cos(kθ) dθ,
and
0
βk =
1
π
2π
f (θ) sin(kθ) dθ .
0
2. (Systems analysis) In systems analysis, a “system” S transforms any “input” function f I to a corresponding “output” function f O = S[ f I ]. Often, the output of S can be described by the convolution
formula
fO = h ∗ fI
where h is some fixed function called the impulse response of the system. The corresponding Fourier
transform H = F[h] is the system’s transfer function. By the convolution identity, the output is
also given by
f O = F −1 [H F I ]
where F I = F[ f I ] .
Two such systems are
r A delayed output system, for which
h(x) = δT (x)
for some T > 0. Then
and
f O (x) = f I ∗ δT (x) =
∞
−∞
H(ω) = e −i 2π T ω
f I (x − y)δT (y)d y = f I (x − T ) .
r [ZTF98, p. 178] or [How01, p. 477] An ideal low pass filter, for which
h(x) =
sin(2π x)
πx
and
H(y) = pulse (y)
58-8
Handbook of Linear Algebra
for some > 0. If an input function f I is periodic with Fourier series
∞
H F I = pulse ·
k∈Z c k e
i 2π ωk x
, then
c k δωk
k=−∞
=
∞
c k pulse (ωk )δωk
k−∞
=
∞
ck
k−∞
1
if |ωk | ≤ δω k =
0 if |ωk | > c k δωk .
|ωk |≤
Thus,
f O (x) = F
−1
[H F I ]|x = F
−1
|ωk |≤
c k δωk =
c k e i 2π ωk x .
x
|ωk |≤
3. [ZTF98, pp. 520–521] (Deconvolution) Suppose S is a system given by
S[ f I ] = h ∗ f I ,
but with h and H = F[h] being unknown. Since
f O = F −1 [H F I ]
where F I = F[ f I ] ,
both h and H can be determined as follows:
r Find the output f = S[ f ] for some known input f for which F = F[ f ] is never zero.
O
I
I
I
I
r Compute
H=
FO
FI
where F O = F[ f O ] .
r Compute h = F −1 [H].
Similarly, an input f I can be reconstructed from an output f O by
f I = F −1
FO
H
provided the transfer function H is known and is never zero.
58.3
The Discrete Theory
Definitions:
In all the following, N ∈ N, and the indexing of any “N items” (including rows and columns of matrices)
will run from 0 to N − 1.
An Nth order sequence is an ordered list of N complex numbers,
(c 0 , c 1 , c 2 , . . . , c N−1 ).
58-9
Fourier Analysis
Such a sequence will often be written as the column vector
c = [c 0 , c 1 , c 2 , . . . , c N−1 ]T ,
and the kth component of the sequence will be denoted by either c k or [c]k . In addition, any such sequence
will be viewed as part of an infinite repeating sequence
(. . . , c −1 , c 0 , c 1 , c 2 , . . . , c N−1 , c N , . . .)
in which c N+k = c k for all k ∈ Z.
Let c be an Nth order sequence. The (Nth order) discrete Fourier transform (DFT) FN c and the (Nth
order) inverse discrete Fourier transform (inverse DFT) FN−1 c of c are the two Nth order sequences
given by
N−1
1 c k e −i 2πkn/N
[FN c]n = √
N k=0
and
N−1
1 [FN−1 c]k = √
c n e i 2πnk/N .
N n=0
Additionally, the terms “discrete Fourier transform” and “inverse discrete Fourier transform” can refer to
the processes for computing these sequences from c.
Given two Nth order sequences a and b, the corresponding product ab and convolution a ∗ b are the
Nth order sequences given by
[ab]k = ak bk
N−1
1 [a ∗ b]k = √
ak− j b j .
N j =0
and
Facts:
All the following facts can be found in [BH95, Chap. 2, 3] or [How01, Chap. 38, 39].
1. Warning: Slight variations of the above summation formulas, e.g.,
N−1
c k e −i 2π kn/N
and
N−1
1 c n e i 2π nk/N ,
N
n=0
k=0
are also often used to define, respectively, FN c and FN−1 c (or even FN−1 c and FN c).
2. FN and FN−1 are one-to-one linear transformations from C N onto C N . They preserve the standard
inner product and are inverses of each other.
3. Letting w = e −i 2π/N , the matrices (with respect to the standard basis for C N ) for the Nth order
discrete Fourier transforms (also denoted by FN and FN−1 ) are
⎡
1
⎢1
⎢
⎢
⎢1
1
⎢
FN = √ ⎢1
N⎢
⎢
⎢ ..
⎣.
1
1
1
···
w 1·1
w 1·2
w 1·3
···
2·1
2·2
2·3
···
w 2(N−1) ⎥
w 3(N−1)
..
.
w
w
⎤
1
w
w
w 3·1
..
.
w 3·2
..
.
w 3·3
..
.
···
(N−1)1
(N−1)2
(N−1)3
···
w
w
..
.
1
w 1(N−1) ⎥
⎥
w (N−1)(N−1)
⎥
⎥
⎥
⎥
⎥
⎥
⎦
58-10
Handbook of Linear Algebra
and
⎡
FN−1
⎢1
⎢
⎢
1
1 ⎢
⎢
= √ ⎢1
N⎢
⎢
⎢ ..
⎣.
1
w
⎤
1
1
1
···
w −1·1
w −1·2
w −1·3
···
w −2·1
w −2·2
w −2·3
···
w −2(N−1) ⎥
w −3·1
..
.
w −3·2
..
.
w −3·3
..
.
···
w −3(N−1)
..
.
−(N−1)1
−(N−1)2
−(N−1)3
···
1
w
w
..
.
1
w −1(N−1) ⎥
⎥
⎥
⎥
⎥.
⎥
⎥
⎥
⎦
w −(N−1)(N−1)
∗
4. FN and FN−1 are symmetric matrices, and are complex conjugates of each other (i.e., FN = FN−1 ).
5. (Convolution identities) Let v, V, w, and W be Nth order sequences with V = FN v and W = FN w.
Then
FN [vw] = V ∗ W
FN [v ∗ w] = VW.
and
Examples:
1. The matrices for the four lowest order DFTs are F1 = [1],
1
1
F2 = √
2 1
1
⎡
−1
,
2
√
1 ⎢
F3 = √ ⎣2 −1 − 3
√
2 3
2 −1 + 3
and
⎡
F4 =
2
1
1
1⎢
⎢1 −i
⎢
2 ⎣1 −1
1
i
1
1
2
⎤
i ⎥
⎥
⎥.
−1⎦
1
1
−1
−i
2. The DFT of v = [1, 2, 3, 4]T is V = [V0 , V1 , V2 , V3 ]T where
⎡ ⎤
⎡
⎤⎡ ⎤ ⎡
⎤
V0
1 1
1
1
1
5
⎢ V ⎥ 1 ⎢1 −i 1
⎢ ⎥ ⎢
⎥
i ⎥
⎢ 1⎥
⎢
⎥ ⎢2⎥ ⎢ 2 + i ⎥
⎢ ⎥= ⎢
⎥⎢ ⎥ = ⎢
⎥.
⎣ V2 ⎦ 2 ⎣1 −1 1 −1⎦ ⎣3⎦ ⎣ −1 ⎦
1 i
−1 −i
4
−1 − i
V3
Applications:
1. [CG00, p. 153] Assume
f (z) =
2M
fn zn
n=0
is the product of two known Mth order polynomials
g (z) =
M
k=0
g k zk
and
h(z) =
M
k=0
⎤
√ ⎥
−1 + 3⎦,
√
−1 − 3
hk zk .
58-11
Fourier Analysis
The coefficients of f form a sequence f = [ f 0 , f 1 , f 2 , . . . , f 2M ]T of order N = 2M + 1.
These coefficients can be computed from the coefficients of g and h using either of the following
approaches:
r Let g = [g , g , g , . . . , g ]T and h = [h , h , h , . . . , h ]T be the Nth order sequences in
0
1
2
2M
0
1
2
2M
which g k and h k are the corresponding coefficients of the polynomials g and h when k ≤ M and
are zero when M < k ≤ N. The coefficients of f can then be computed by
fk =
2M
g k− j h j =
√
N[g ∗ h]k
for k = 0, 1, 2, . . . , 2M .
j =0
r Let F = [F , F , F , . . . , F ]T be the Nth order DFT of f, and note that, for n = 0, 1, 2, . . .,
0
1
2
2M
N − 1,
N−1
1 1
1
f k e −i 2πkn/N = √ f (w n ) = √ h(w n )g (w n )
Fn = √
N k=0
N
N
where w = e −i 2π/N . Thus, the coefficients of f can be computed by first using the formulas for
g and h to compute
1
F n = √ h(w n )g (w n )
N
for n = 0, 1, 2, . . . , N − 1 ,
and then taking the inverse DFT of F = [F 0 , F 1 , F 2 , . . . , F N−1 ]T .
2. [BH95, p. 247] Consider finding the solution v = [v 0 , v 1 , v 2 , . . . , v N−1 ]T to the difference equation
αv k−1 + βv k + αv k+1 = f k
for n = 0, 1, 2, . . . , N − 1
where α and β are constants and f = [ f 0 , f 1 , f 2 , . . . , f N−1 ]T is a known Nth order sequence. If
the boundary conditions are periodic (i.e., v 0 = v N and v −1 = v N−1 ), then setting
v = FN−1 V
f = FN−1 F
and
yields the sequence v that satisfies the periodic boundary conditions, and whose DFT, V, satisfies
N−1
C n e i 2πnk/N = 0
for k = 0, 1, 2, . . . , N − 1
n=0
where
2πn
C n = 2α cos
+ β Vn − F n .
N
From this, it follows that the solution v to the original difference equation is the DFT of the sequence
[V0 , V1 , V2 , . . . , VN−1 ]T given by
Vn =
2α cos
F
2πn n N
+β
,
provided the denominator never vanishes.
3. A difference equation of the form just considered with periodic boundary conditions arises when
considering the steady-state temperature distribution on a uniform ring containing heat sources
and sinks. In this case, the temperature v, as a function of angular position θ, is modeled by the
58-12
Handbook of Linear Algebra
one-dimensional Poisson’s equation
d 2v
= f (θ) ,
dθ 2
where f (θ) describes the heat source/sink density at angular position θ. The discrete analog of this
equation is
v k−1 − 2v k + v k+1 = f k
for n = 0, 1, 2, . . . , N − 1
where
v = [v 0 , v 1 , v 2 , . . . , v N−1 ]T
and
f = [ f 0 , f 1 , f 2 , . . . , f N−1 ]T
describe, respectively, the temperatures at N evenly space positions around the ring, and the net
sources and sinks of thermal energy about these positions. Setting
v = FN−1 V
and
f = FN−1 F
and applying the formulas given above yields
2π n
2 cos
− 2 Vn = F n
N
for n = 0, 1, 2, . . . , N − 1.
The coefficient on the left is nonzero when n = 0. Thus,
Vn =
2 cos
F
2πnn N
−2
for n = 1, 2, . . . , N − 1.
However, for n = 0,
N−1
N−1
1 1 f k e −i 2πk·0/N = √
fk ,
0 · V0 = F 0 = √
N k=0
N k=0
pointing out that, for an equilibrium temperature distribution to exist, the net applied heat energy,
N−1
k=0 f k , must be zero. Assuming this, the last equation then implies that V0 is arbitrary, which,
since
N−1
N−1
1 1 v k e −i 2πk·0/N = √
vk
V0 = √
N k=0
N k=0
means that the given conditions are not sufficient to determine the average temperature throughout
the ring.
58.4
Relating the Functional and Discrete Theories
Definitions:
Let f be a (continuous) function on a finite interval [0, L ]. For any N ∈ N and x > 0 satisfying
Nx ≤ L , the corresponding (Nth order) sampling (with spacing x) is the sequence
f = [ f 0 , f 1 , f 2 , . . . , f N−1 ]T
where f k = f (kx) .
Corresponding to this is the scaled sampling
f = [ f 0, f 1, f 2, . . . , f N−1 ]T
58-13
Fourier Analysis
and the discrete approximation
∞
f =
f k δkx
k=−∞
where, in both,
fk =
f k x
for k = 0, 1, 2, . . . , N − 1
f k+N
in general
.
Facts:
1. [How01, pp. 713–715] Let f be the discrete approximation of a continuous function f based on
an Nth order sampling with spacing x. Then f is a periodic array of delta functions with spacing
x and index period N that approximates f over the interval (0, Nx). In particular, for any
other function ψ which is continuous on [0, Nx],
Nx
Nx
ψ(x) f (x) d x ≈
0
ψ(x) f (x) d x =
N−1
0
ψ(kx) f k .
k=0
Attempting to approximate f with f outside the interval (0, Nx), however, cannot be justified
unless f is also periodic with period Nx.
2. [How01, chap. 38] Let
f = [ f 0 , f 1 , f 2 , . . . , f N−1 ]T
and
F = [F 0 , F 1 , F 2 , . . . , F N−1 ]T
and
F =
be two Nth order sequences, and let
f =
∞
f k δkx
∞
F n δnω
n=−∞
k=−∞
be two corresponding periodic arrays of delta functions with index period N, and with spacings
x and ω satisfying xω = 1/N. Then
√
F = F[ f ] ⇐⇒ F = ω N FN f .
In particular, if ω = x = N −1/2 , then
F = F[ f ]
⇐⇒
F = FN f .
3. [How01, pp. 719–723] Suppose f is a continuous, piecewise smooth function that vanishes outside
the finite interval (0, L ). Let F = F[ f ], and let f be the Nth order scaled sampling of f with
spacing x chosen so that L = (N − 1)x. Then, for each n ∈ Z,
√
F (nω) ≈ N [FNf]n
where ω = (Nx)−1 . The error in this approximation is bounded by
(max | f | + 2π|n|ω max | f |)
L2N
.
2(N − 1)2
In practice, this bound may significantly overestimate the actual error.
58-14
Handbook of Linear Algebra
4. [BH95, pp. 181–188]
If f is a continuous, piecewise smooth, periodic function with period p and
Fourier series n∈Z c n e i 2π ωn x , then
1
c n ≈ √ [FN f]n
N
for −
N
N
<n≤
2
2
where f is the Nth order sampling of f over [0, p] with spacing x = p/N. Moreover, if f is
m-times differentiable on R and f (m) is piecewise continuous and piecewise monotone for some
m ∈ N, then the error in this approximation is bounded by
β
N m+1
for some positive constant β independent of n and N.
Examples:
1. Let
f =
∞
f k δk/3
k=−∞
be the periodic array with index period N = 4 and coefficients
f 0 = 0,
f 1 = 1,
f 2 = 2,
and
f3 = 1 .
The Fourier transform of f is then the periodic array
F =
∞
F n δnω
n=−∞
with index period 4, spacing ω = 3/4, and coefficients given by
⎡ ⎤
⎡
⎤⎡ ⎤
⎡ ⎤
F0
1 1
1 1
0
6
⎢F ⎥
⎢
⎥
⎥
⎢
⎥
⎢
√
1
−i
1
i
1
3⎢
⎢ 1⎥
⎥⎢ ⎥ 1 ⎢ 3 ⎥
⎢ ⎥ = ω 4 F4 f = ⎢
⎥⎢ ⎥ = ⎢ ⎥ .
4 ⎣1 −1
⎣ F2 ⎦
1 −1⎦ ⎣2⎦ 2 ⎣ 0 ⎦
1 i
−1 −i
F3
−3
1
2. The 8th order sampling with spacing x = 1/6 of
⎧
⎪
⎨ x
f (x) = 1 − x
⎪
⎩
0
for 0 ≤ x ≤ 1/2
for 1/2 ≤ x ≤ 1
otherwise
is f = [ f 0 , f 1 , f 2 , . . . , f 7 ]T with
⎧
for k = 0, 1, 2, 3
⎪
⎨ k/6
k
.
= 1 − k/6 for k = 4, 5, 6
fk = f
⎪
6
⎩
0
for k = 7
That is,
f = 0,
T
1 2 3 2 1
, , , , , 0, 0
6 6 6 6 6
.
58-15
Fourier Analysis
Multiplying this by the spacing x = 1/6 yields the scaled sampling
T
f = 0, 1 , 2 , 3 , 2 , 1 , 0, 0 .
36 36 36 36 36
Letting F = F[ f ], ω = (Nx)−1 = 3/4, and w = e −i 2π/8 ,
F (nω) ≈
√
8 [F8f]n
=
1 0n
0w + 1w 1n + 2w 2n + 3w 3n + 2w 4n + 1w 5n + 0w 6n + 0w 7n
36
=
1 3n −2n
w w
+ 2w −n + 3 + 2w n + w 2n ,
36
which simplifies to
F (nω) ≈
√
nπ !
nπ !"
1
8 [F8f]n = e −i 3nπ/4 3 + 2 cos
+ 4 cos
.
36
2
4
The error bound on this approximation is
E n = (max | f | + 2π|n|ω max | f |)
L2N
4 + 3π|n|
.
=
2(N − 1)2
36
In fact, the Fourier transform of f is easily found to be
F (ω) = e −i π ω
sin(π ω/2)
πω
2
.
So
F (nω) = F
3n
4
= e −i 3nπ/4
4 sin(3nπ/8)
3nπ
2
,
and the actual error in the approximation is
εn = F (nω) −
√
8 [F8f]n sin(nπ/8)
1 =
64
36 nπ
2
nπ !
nπ !
− 3 − 2 cos
− 4 cos
.
2
4 For n = 0 through 4, the actual and approximate values of F (nω) (to five decimal places) are
F (0ω) = 0.25000
F (1ω) = − 0.10872(1 + i )
F (2ω) = 0.02252i
F (3ω) = 0.00207(1 − i )
F (4ω) = − 0.01126
√
8 [F8f]0 = 0.25000
√
8 [F8f]1 = − 0.11448(1 + i )
√
8 [F8f]2 = 0.02778i
√
8 [F8f]3 = 0.00337(1 − i )
√
8 [F8f]3 = − 0.02778.
58-16
Handbook of Linear Algebra
The corresponding error bounds and actual errors (to five decimal places) are
E 0 = 0.11111
ε0 = 0.00000
E 1 = 0.37291
ε1 = 0.00815
E 2 = 0.63471
ε2 = 0.00526
E 3 = 0.89651
ε3 = 0.00183
E 4 = 1.15831
ε4 = 0.01652.
Applications:
1. Let f be a continuous, piecewise smooth function to be sampled that is nonzero only in the interval
(0, 1) and which satisfies
| f (x)| ≤
1
4
| f (x)| ≤ 1
and
∀x ∈ R.
Consider the problem of approximating F (nω) where F = F[ f ], ω = 1/2, and n is any
integer from 0 to 10. For each N ∈ N, set
x =
2
1
=
Nω
N
L = (N − 1)x =
and
2(N − 1)
,
N
and let f = [ f 0, f 1, f 2, . . . , f N−1 ]T be the corresponding scaled sampling of f ,
f k = f (kx)x = f
2k
N
·
2
.
N
Then, for n = 0, 1, 2, . . . , 10,
F (nω) ≈
√
N−1
2 2k
N [F N f]n =
f
e −i 2πkn/N
N
N
k=0
with an error bound of
(max | f | + 2π |n|ω max | f |)
L2N
4 + 10π
18
=
<
.
2(N − 1)2
2N
N
2. [CG00, p. 157] (Deconvolution) Assume g = f ∗h where g and h are known by Nth order samplings
g and h, each with spacing x. The deconvolution formula for finding f ,
f = F −1 [V ]
where V =
F[g ]
G
=
,
H
F[h]
is approximated by
f (kx) ≈
√
N FN−1 V
"
k
is the sequence given by
where, letting ω = (Nx)−1 , V
√
Vn = ω √
"
"n
N FN [
g]
N FN [
h]
n
for n = 0, 1, 2, . . . , N − 1 .
58-17
Fourier Analysis
This reduces to
f (kx) ≈
1
√
x N
"
FN−1 V
k
where, letting G = FN g and H = FN h,
Vn =
Gn
Hn
for n = 0, 1, 2, . . . , N − 1 .
This requires, of course, that each Hn be nonzero.
58.5
The Fast Fourier Transform
A fast Fourier transform (FFT) is any of a number of algorithms for computing DFTs in which the number
of arithmetic computations is greatly reduced through clever use of symmetries and cyclic structures in
the DFT matrices. The FFT described here is the standard radix 2 FFT for computing the Nth order
“alternative” discrete Fourier transform D N v given by
[D N v]n =
√
N[FN v]n =
N−1
v k e −i 2πnk/N .
k=0
For additional details on implementing a radix 2 FFT and descriptions of other FFTs, see [BH95,
chap. 10], [CG00], or [Wal96].
Algorithms:
Two algorithms for computing the alternative DFT
V = [V0 , V1 , V2 , . . . , VN−1 ]T
of an Nth order sequence
v = [v 0 , v 1 , v 2 , . . . , v N−1 ]T
are given. The first is the “first level” radix 2 FFT illustrating the basic concepts. The second is the more
complete radix 2 FFT.
1. (First level radix 2 FFT) This requires N = 2M for some M ∈ N:
r First, split v into two Mth order sequences v O and v E composed, respectively, of the even-indexed
and odd-indexed elements of v,
E
v E = [v 0E , v 1E , v 2E , . . . , v M−1
]T = [v 0 , v 2 , v 4 , . . . , v 2M−2 ]T
and
O
v O = [v 0O , v 1O , v 2O , . . . , v M−1
]T = [v 1 , v 3 , v 5 , . . . , v 2M−1 ]T .
r Then compute the Mth order alternative DFTs
E
]T = V E = D N v E
[V0E , V1E , V2E , . . . , VM−1
and
O
[V0O , V1O , V2O , . . . , VM−1
]T = V O = D N v O .
58-18
Handbook of Linear Algebra
r Then construct the Nth order sequence V = [V , V , V , . . . , V
T
E
O
0
1
2
N−1 ] from V and V via the
butterfly relations
Vn =
VnE + w n VnO
E
Vn−M
−w
n−M
O
Vn−M
for n = 0, 1, 2, . . . , M − 1
for n = M, M + 1, M + 2, . . . , N − 1
where w = e −i 2π/N .
2. (Full radix 2 FFT) This requires N = 2 P for some P ∈ N. To simplify the description, let
σk ∈
{0}
{E , O}
if k = 0
if k = 1, 2, 3, . . . , P
.
r Recursively, split v into the sequences in the set
{vσ0 σ1 σ2 ···σ K : K = 0, 1, 2, . . . , P }
where vσ0 = v and, letting = σ0 σ1 σ2 · · · σ K and M = 2 P −K −1 ,
T
T
E
v E = v 0 E , v 1 E , v 2 E , . . . , v M−1
T
T
= v 0 , v 2 , v 4 , . . . , v 2M−2
and
O
v O = v 0 O , v 1 O , v 2 O , . . . , v M−1
= v 1 , v 3 , v 5 , . . . , v 2M−1
.
r For each of the 2 P first order sequences in {vσ0 σ1 σ2 ···σ P = [v σ1 σ2 ···σ P ]} , set
Vσ1 σ2 ···σ P = [V σ1 σ2 ···σ P ] = [v σ1 σ2 ···σ P ] = vσ1 σ2 ···σ P
r Set V = Vσ0 where
{Vσ0 σ1 σ2 ···σ K : K = 0, 1, 2, . . . , P }
is the set of sequences recursively constructed via the butterfly relations
Vn =
Vn E + w n Vn O
E
Vn−M
−w
n−M
O
Vn−M
for n = 0, 1, 2, . . . , M − 1
for n = M, M + 1, M + 2, . . . , 2M − 1
where = σ0 σ1 σ2 · · · σ K , M = 2 P −K −1 , and w = e (−i 2π/N)2 .
K
Facts:
All the following facts except those with specific reference can be found in [BH95], [CG00], or [Wal96].
1. Let V be the sequence computed from an Nth order sequence v by either of the two algorithms
above. Then V = D N v. Multiplying V by N −1/2 yields FN v.
2. For corresponding “fast” algorithms for computing the inverse DFT of v, replace w = e −i 2π/N and
K
K
w = e (−i 2π/N)2 , respectively, with w = e i 2π/N and w = e (i 2π/N)2 in the above algorithms.
P
3. [BH95, p. 393] Assume N = 2 1 for some P ∈ N. Using just the defining formulas, the
computation of an Nth order DFT requires approximately 2N 2 arithmetic operations ((N − 1)2
multiplications and N(N − 1) additions). Using the above described radix 2 FFT allows this same
DFT to be computed with approximately (3N/2) log2 N arithmetic operations (N log2 N additions
and (N/2) log2 N multiplications, plus another N multiplications by N −1/2 if it is desired to convert
58-19
Fourier Analysis
the alternative DFT to the one previously described). Thus, the ratio of the number of arithmetic
operations used to compute an Nth order DFT by the two methods is
3 log2 N
operations if using the FFT
≈
.
operations if not using the FFT
4N
4. Assume N = 2 P 1 for some P ∈ N. Simply using the basic formulas, the computation of the
discrete convolution of two sequences of order N,
N−1
1 v k− j w j ,
[v ∗ w]k = √
N j =0
requires approximately 2N 2 arithmetic operations ((N + 1)N multiplies and (N − 1)N additions).
On the other hand, this convolution can be computed using the convolution identity
v ∗ w = FN−1 [VW]
where
V = FN v
and W = FN w .
Employing this identity along with the radix 2 FFT to compute the three DFTs involved requires
approximately
r 3 × (3N/2) log N arithmetic operations for the DFTs,
2
r Plus N multiplications for the product VW,
for a total of approximately (9N/2) log2 N arithmetic operations. The ratio of operations required
to compute this convolution by the two methods is
9 log2 N
operations if using identities and the FFT
≈
.
operations if not using identities and the FFT
4N
Examples:
1. Consider finding the DFT of
v = [v 0 , v 1 , v 2 , v 3 ]T = [1, 2, 3, 4]T
using the full radix 2 FFT. The order here is N = 4 (so P = 2). Recursively splitting v into the even
and odd index subsequences yields
v E = [1, 3]T
and
v O = [2, 4]T
and then
v E E = [1],
v E O = [3],
v O E = [2],
and v O O = [4].
(The σ0 = 0 superscript is being suppressed.) Setting Vσ1 σ2 = vσ1 σ2 yields
V E E = [1],
V E O = [3],
V O E = [2],
and V O O = [4].
58-20
Handbook of Linear Algebra
Applying the first round of butterfly relations (in which K = 1, M = 2 P −K −1 = 1 and w =
1
e (−i 2π/4)2 = −1),
V0E = V0E E + w 0 V0E O = 1 + (−1)0 3 = 4
V1E = V0E E − w 1−1 V0E O = 1 − (−1)0 3 = −2
V0O = V0O E + w 0 V0O O = 2 + (−1)0 4 = 6
V1O = V0O E − w 1−1 V0O O = 2 − (−1)0 4 = −2 ,
yields the second order sequences
V E = [4, −2]T
and
V O = [6, −2]T .
Applying the second (and last) round of butterfly relations (in which K = 0, M = 2 P −K −1 = 2
0
and w = e (−i 2π/4)2 = −i ),
V0 = V0E + w 0 V0O = 4 + (−i )0 6 = 10
V1 = V1E + w 1 V1O = −2 + (−i )1 (−2) = −2 + 2i
E
O
V2 = V2−2
− w 2−2 V2−2
= 4 − (−i )0 6 = −2
E
O
− w 3−2 V3−2
= −2 − (−i )1 (−2) = −2 − 2i ,
V3 = V3−2
yields the alternative DFT
V = Dv = [ 10 , −2 + 2i , −2 , −2 − 2i ]T .
2. Consider the computation of an Nth order DFT. When N = 1024 = 210 , the ratio of operations
required to compute this DFT with and without using the FFT is
3 · 10
operations if using the FFT
≈
≈ 7.3 × 10−3 .
operations if not using the FFT
4 · 210
When N = 11,048,576 = 220 , this ratio is
3 · 20
operations if using the FFT
≈
≈ 1.4 × 10−5 .
operations if not using the FFT
4 · 220
3. Consider computing the convolution of two Nth order sequences. When N = 1024 = 210 , the ratio
of operations required to compute this convolution with and without using the FFT and identities
is
9 · 10
operations if using identities and the FFT
≈
≈ 2.2 × 10−2 .
operations if not using identities and the FFT
4 · 210
When N = 11,048,576 = 220 , this ratio is
9 · 20
operations if using identities and the FFT
≈
≈ 4.3 × 10−5 .
operations if not using identities and the FFT
4 · 220
Fourier Analysis
58-21
References
[BH95] W. Briggs and V. Henson. The DFT: An Owner’s Manual for the Discrete Fourier Transform. SIAM,
Philadelphia, 1995.
[BC01] J. Brown and R. Churchill. Fourier Series and Boundary Value Problems, 6th ed., McGraw-Hill, New
York, 2001.
[CG00] E. Chu and A. George. Inside the FFT Black Box. CRC Press, Boca Raton, FL, 2000.
[Col88] D. Colton. Partial Differential Equations: An Introduction. Random House, New York, 1988.
[How01] K. Howell. Principles of Fourier Analysis. Chapman & Hall/CRC, Boca Raton, FL, 2001.
[Str94] R. Strichartz. A Guide to Distribution Theory and Fourier Transforms. CRC Press, Boca Raton, FL,
1994.
[Wal96] J. Walker. Fast Fourier Transforms, 2nd ed., CRC Press, Boca Raton, FL, 1996.
[ZTF98] R. Ziemer, W. Trantor, and D. Fannin. Signals and Systems: Continuous and Discrete, 4th ed.,
Prentice Hall, Upper Saddle River, NJ, 1998.