INTRODUCTION

Consider a function f that is defined and continuous on the interval [0, 2p]. If x₀, x₁, x₂, … , x_n, … are equally spaced points in the interval, then the corresponding function values f₀, f₁, f₂, … , f_n, … shown in FIGURE 15.6.1 are said to represent a discrete sampling of the function f. The notion of discrete samplings of a function is important in the analysis of continuous signals.

In this section, the complex or exponential form of a Fourier series plays an important role in the discussion. A review of Section 12.4 is recommended.

Discrete Fourier Transform

Consider a function f defined on the interval [0, 2p]. From (11) of Section 12.4 we saw that f can be written in a complex Fourier series,

(1)

where the ω = 2π/2p = π/p is the fundamental angular frequency and 2p is the fundamental period. In the discrete case, however, the input is f₀, f₁, f₂, … , which are the values of the function f at equally spaced points x = nT, n = 0, 1, 2, … . The number T is called the sampling rate or the length of the sampling interval.^* If f is continuous at T, then the sample of f at T is defined to be the product f(x)δ(x − T), where δ(x − T) is the Dirac delta function (see Section 4.5). We can then represent this discrete version of f, or discrete signal, as the sum of unit impulses acting on the function at x = nT:

(2)

If we apply the Fourier transform to the discrete signal (2), we have

(3)

By the sifting property of the Dirac delta function (see Section 4.5), (3) is the same as

(4)

The expression F(α) in (4) is called the discrete Fourier transform (DFT) of the function f. We often write the coefficients f(nT) in (4) as f(n) or f_n. It is also worth noting that since e^iαx is periodic in α and e^iαT = e^i(αT+2π) = e^i(α+2π/T)T, we only need to consider the function for α in [0, 2π/T]. Let N = 2π/T. This places x in the interval [0, 2π]. So, because we sample over one period, the sum in (4) is actually finite.

Now consider the function values f(x) at N equally spaced points, x = nT, n = 0, 1, 2, … , N − 1, in the interval [0, 2π]; that is, f₀, f₁, f₂, … , f_N−1. The (finite) discrete Fourier series f(x) = using these N terms gives us

If we let and use the usual laws of exponents, this system of equations is the same as

(5)

If we use matrix notation (see Sections 8.1 and 8.2), then (5) is

(6)

Let the N × N matrix in (6) be denoted by the symbol F_N. Given the inputs f₀, f₁, f₂, … , f_N−1, is there an easy way to find the Fourier coefficients c₀, c₁, c₂, … , c_N−1? If is the matrix consisting of the complex conjugates of the entries of F_N and if I denotes the N × N identity matrix, then we have

It follows from (6) and the last equation that

Discrete Transform Pair

Recall from Section 15.4 that in the Fourier transform pair we use a function f(x) as input and compute the coefficients that give the amplitude for each frequency k (c_k in the case of periodic functions of period 2π) or we compute the coefficients that give the amplitude for each frequency α (F(α) in the case of nonperiodic functions).

Also, given these frequencies and coefficients, we could reconstruct the original function f(x). In the discrete case, we use a sample of N values of the function f(x) as input and compute the coefficients that give the amplitude for each sampled frequency. Given these frequencies and coefficients, it is possible to reconstruct the n sampled values of f(x). The transform pair, the discrete Fourier transform pair, is given by

(7)

where

EXAMPLE 1 Discrete Fourier Transform

Let N = 4 so that the input is f₀, f₁, f₂, f₃ at the four points x = 0, π/2, π, 3π/2. Since ω₄ = e^iπ/2 = = i, the matrix F₄ is

Hence from (7), the Fourier coefficients are given by c = f :

If we let f₀, f₁, f₂, f₃ be 0, 2, 4, 6, respectively, we find from the preceding matrix product that

Note that we obtain the same result using (4); that is, F(α) = (nT)e^iαnT, with T = π/2, α = 0, 1, 2, 3. The graphs of |c_n|, n = 0, 1, 2, 3, or equivalently |F(α)| for α = 0, 1, 2, 3, are given in FIGURE 15.6.2. ≡

Finding the coefficients involves multiplying by matrices F_n and _n. Because of the nature of these matrices, these multiplications can be done in a very computationally efficient manner, using the Fast Fourier Transform (FFT), which is discussed later in this section.

Heat Equation and Discrete Fourier Series

If the function f in the initial-value problem

(8)

is periodic with period 2π, the solution can be written in terms of a Fourier series for f(x). We can also approximate this solution with a finite sum

If we examine both sides of the one-dimensional heat equation in (8), we see that

and

since .

Equating these last two expressions, we have the first-order DE

The final task is to find the values c_j(0). However, recall that u(x, t) = (t)e^ikx and u(x, 0) = f(x), so c_j (0) are the coefficients of the discrete Fourier series of f(x). Compare this with Section 13.3.

Heat Equation and Discrete Fourier Transform

The initial-value problem (8) can be interpreted as the mathematical model for the temperature u(x, t) in an infinitely long bar. In Section 15.4 we saw that we can solve (8) using the Fourier transform and that the solution u(x, t) depends on the Fourier transform F(α) of f(x) (see pages 794–795). We can approximate F(α) by taking a different look at the discrete Fourier transform.

First we approximate values of the transform by discretizing the integral {f(x)} = F(α) = . Consider an interval [a, b]. Let f(x) be given at n equally spaced points

Now approximate:

If we now choose a convenient value for α, say with M an integer, we have

where recall that ω_n = e^i2π/n. This is a numerical approximation to the Fourier transform of f(x) evaluated at points with M an integer.

EXAMPLE 2 Example 1, Section 15.4—Revisited

Recall from Example 1 in Section 15.4 (with u₀ = 1) that the Fourier transform of a rectangular pulse defined by

is

The frequency spectrum is the graph of |F(α)| versus α given in FIGURE 15.6.3(a). Using n = 16 equally spaced points between a = −2 and b = 2, and M running from −6 to 6, we get the discrete Fourier transform of f(x), superimposed over the graph of |F(α)|, in Figure 15.6.3(b).

≡

Aliasing

A problem known as aliasing may appear whenever one is sampling data at equally spaced intervals. If you have ever seen a motion picture where rotating wheels appear to be rotating slowly (or even backwards!), you have experienced aliasing. The wheels may rotate at a high rate, but because the frames in a motion picture are “sampled” at equally spaced intervals, we see a low rate of rotation.

Graphing calculators also suffer from aliasing due to the way they sample points to create graphs. For example, plot the trigonometric function y = sin 20πx with frequency 10 on a Texas Instruments TI-92 and you get the nice graph in FIGURE 15.6.4(a). At higher frequencies, say y = sin 100πx with frequency 50, you get the correct amount of cycles, but the amplitudes of the graph in Figure 15.6.4(b) are clearly not 1.

On a calculator such as the Texas Instruments TI-83, the graphs in FIGURE 15.6.5 show aliasing much more clearly.

The problem lies in the fact that e^2nπi = cos 2nπ + i sin 2nπ = 1 for all integer values of n. The discrete Fourier series cannot distinguish e^inx from 1 as these functions are equal at sampled points x = . The higher frequency is seen as the lower one. Consider the functions cos () and cos (). If we sample at the points n = 0, 1, 2, … , these two functions appear the same, the lower frequency is assumed, and the amplitudes (Fourier coefficients) associated with the higher frequencies are added in with the amplitude of the low frequency. If these Fourier coefficients at large frequencies are small, however, we do not have a big problem. In the Sampling Theorem below, we will see what can take care of this problem.

Signal Processing

Beyond solving PDEs as we have done earlier, the ideas of this section are useful in signal processing. Consider the functions we have been dealing with as signals from a source. We would like to reconstruct a signal transmitted by sampling it at discrete points. The problem of calculating an infinite number of Fourier coefficients and summing an infinite series to reconstruct a signal (function) is not practical. A finite sum could be a decent approximation, but certain signals can be reconstructed by a finite number of samples.

If a signal f(x) is band-limited; that is, if the range of frequencies of the signal lie in a band −A < k < A, then the signal can be reconstructed by sampling two times for every cycle of the highest frequency present; in fact,

To justify Theorem 15.6.1, consider the Fourier transform F(α) of f(x) as a periodic extension so that F(α) is defined for all values of α, not just those in −A < k < A. Using the Fourier transform, we have

(9)

(10)

Treating F(α) as a periodic extension, the Fourier series for F(α) is

(11)

where (12)

Using (10), note that

which by (12) is equal to c_n. Substituting c_n = into (11) yields

Substituting this expression for F(α) back into (10), we have

Note that we used, in succession, an interchange of summation and integration (not always allowed, but is acceptable here), integration of an exponential function, sin θ = , and the fact that sin(−θ) = −sin θ.

So, from samples at intervals of π/A, all values of f can be reconstructed. Note that if we allow e^iAx (in other words, we allow k = A), then the Sampling Theorem will fail. If, for example, f(x) = sin Ax, then all samples will be 0 and f cannot be reconstructed, as aliasing appears again.

Band-Limited Signals

A signal that contains many frequencies can be filtered so that only frequencies in an interval survive, and it becomes a band-limited signal. Consider the signal f(x). Multiply the Fourier transform F(α) of f by a function G(α) that is 1 on the interval containing the frequencies α you wish to keep, and 0 elsewhere. This multiplication of two Fourier transforms in the frequency domain is a convolution of f(x) and g(x) in the time domain. Recall that Problem 20 in Exercises 15.4 states that

The integral on the right-hand side is called the convolution of f and g and is written f*g. The last statement can be written more compactly as

The analogous idea for Laplace transforms is in Section 4.4. The function g(x) = has as its Fourier transform the pulse function

This implies that the function (f*g)(x) is band-limited, with frequencies between −A and A.

Computing with the Fast Fourier Transform

Return to the discrete Fourier transform of f(x), where we have f sampled at n equally spaced points a distance of T apart, namely, 0, T, 2T, 3T, … , (n − 1)T. (We used T = π/n at the beginning of this section.) Substituting this, the discrete Fourier transform

becomes

For simplicity of notation, write this instead as

This should remind you of (6), where we had

or f = F_nc. The key to the FFT is properties of ω_n and matrix factorization. If n = 2^N, we can write F_n in the following way (which we will not prove):

(13)

where I_k is the k × k identity matrix and P is the permutation matrix that rearranges the matrix c so that the even subscripts are ordered on the top and the odd ones are ordered on the bottom. The matrix D is a diagonal matrix defined by

Note that each of the matrices can, in turn, be factored. In the end, the matrix F_n with n² nonzero entries is factored into the product of n simpler matrices at a great savings to the number of computations needed on the computer.

EXAMPLE 3 The FFT

Let n = 2² = 4 and let F₄ be the matrix in Example 1:

From (13), the desired factorization of F₄ is

(14)

We have inserted dashed lines in the matrices marked A and B so that you can identify the submatrices I₂, D₂, −D₂, and F₂ by comparing (14) directly with (13). You are also encouraged to multiply out the right side of (14) and verify that you get F₄. Now if c = , then

Without going into details, a computation of F_n requires n² computations, while using the matrix factorization (the FFT) means the number of computations is reduced to one proportional to n ln n. Try a few larger values of n and you will see substantial savings.

15.6 Exercises Answers to selected odd-numbered problems begin on page ANS-39.

1. Show that .
2. Prove the sifting property of the Dirac delta function:

   

   [Hint: Consider the function

   

   Use the mean value theorem for integrals and then let → 0.]

3. Find the Fourier transform of the Dirac delta function δ(x).
4. Show that the Dirac delta function is the identity under the convolution operation; that is, show f*δ = δ*f = f. [Hint: Use Fourier transforms and Problem 3.]
5. Show that the derivative of the Dirac delta function δ′(x − a) has the property that it sifts out the derivative of a function f at a. [Hint: Use integration by parts.]
6. Use a CAS to show that the Fourier transform of the function g(x) = is the pulse function

   

7. Write the matrix F₈ and then write it in factored form (13). Verify that the product of the factors is F₈. If instructed, use a CAS to verify the result.
8. Let ω_n = e^i2π/n = cos () + i sin (). Since e^i2πk = 1, the numbers , k = 0, 1, 2, … , n − 1, all have the property that ()ⁿ = 1. Because of this, , k = 0, 1, 2, … , n − 1, are called the nth roots of unity and are solutions of the equation zⁿ − 1 = 0. Find the eighth roots of unity and plot them in the xy-plane where a complex number is written z = x + iy. What do you notice?

Computer Lab Assignments

9. Use a CAS to verify that the function f * g, where f(x) = and g(x) = , is band-limited. If your CAS can handle it, plot the graphs of {f*g} and F(α)G(α) to verify the result.
10. If your CAS has a discrete Fourier transform command, choose any six points and compare the result obtained using this command with that obtained from c = ₆f.

* Note that the symbol T used here does not have the same meaning as in Section 12.4.