3

DISCRETE FOURIER ANALYSIS

The Fourier transform and Fourier series techniques are useful for analyzing continuous signals such as the graph in Figure 3.1. However, for many applications, the signal is a discrete data set (see Figure 3.2), such as the signal coming from a compact disc player. A discrete version of the Fourier transform is needed to analyze discrete signals.

3.1 THE DISCRETE FOURIER TRANSFORM

To motivate the idea behind the discrete Fourier transform, we numerically approximate the coefficients of a Fourier series for a continuous function f(t). The trapezoidal rule for approximating the integral (2π)^–1 F(t) dt with step size h = 2π/n is

where Y_j := F(hj) = F(2πj / n), j = 0,..., n. If F(t) is 2π -periodic, then Y₀ = Y_n and the preceding formula becomes

Applying this formula to the computation of the kth complex Fourier coefficient gives

Therefore

(3.1)

where

The sum on the right side of Eq. (3.1) involves yj, which is the value of f at t = 2πj / n. The values of f at the other t values are not needed. This sum will be used in the next section as the definition of the Fourier transform of a discrete signal, whose values may only be known at a discrete set of time nodes t = 2πj / n for j = 0,..., n – 1.

The right side of Eq. (3.1) is unchanged if k is replaced by k + n since = e-2πi = 1 Thus this expression does not approximate a_k for k ≤ n, because the at are not n-periodic. In fact, this expression only approximates a_k for k that are relatively small compared to n because the trapezoidal rule algorithm only provides accurate numerical approximations if the step size, h = 2π / n, is small relative to the frequency k.

3.1.1 Definition of Discrete Fourier Transform

Let S_n be the set of n-periodic sequences of complex numbers. Each element, in Sn, can be thought of as a periodic discrete signal where y_j is the value of the signal at a time node t = t_j. The sequence yj is n -periodic if y_k+n = y_k for any integer k. The set S_n forms a complex vector space under the operations of entry-by-entry addition and entry-by-entry multiplication by a scalar. If and then the 7th component of {x + y} is X_j + y_j and the jth component of c{x} is cx_j . Here, n should be thought of as the number of time nodes corresponding to the discrete signal of interest.

Definition 3.1 Let . The discrete Fourier transform of y is the sequence (F_n{y})k = "%, where

In detail, the discrete Fourier transform is the sequence

The formula for the discrete Fourier transform is analogous to the formula for the kth Fourier coefficient with the sum over j taking the place of the integral over t (see Theorem 1.18). As stated in the previous section, if the yj arise as values from a continuous signal, f, defined on the interval [0, 2π], then the kth Fourier coefficient of f, namely a_k, is approximated by

for k small relative to n.

The computation of the discrete Fourier transform is equivalent to the following matrix computation:

where y = (y₀,..., y_n–1)^T and and where

(3.2)

Note that F_n is a symmetric matrix.

3.1.2 Properties of the Discrete Fourier Transform

The discrete Fourier transform of a signal in S_n produces a discrete signal that also belongs to S_n as the following lemma shows.

Lemma 3.2 F_n is a linear operator from S_n to S_n.

Proof. To show that F_n{y} belongs to S_n we must show that F_n{y}_k is n-periodic. We have

Therefore k+n = k and so this sequence is n-periodic. Finally, since matrix

multiplication is a linear process, F_n is a linear operator.

Our next tak is to compute the inverse of the discrete Fourier transform. We have already computed the inverse Fourier transform since Theorem 2.1 tells us how to recover the function f from its Fourier transform. The inverse discrete Fourier transform is analogous; it allows us to recover the original discrete signal, y, from its discrete Fourier transform, = (y). Computing the inverse discrete Fourier transform is equivalent to computing the inverse of the matrix The following theorem gives an easy formula for the inverse of this matrix.

Theorem 3.3 Suppose y = {y_k} is an element of S_n. Let F_n(y) = . That is,

where w = exp(27u/n). Then y = F -_n 1(y) is given by

Remark. In matrix terminology, Theorem 3.3 asserts that the inverse discrete Fourier transform is given by matrix multiplication with the matrix F_n [F is defined in Eq. (3.2)], that is,

Since = ()(y), the equation y = F_n() is equivalent to

(3.3)

Since F_n is symmetric, this equation means that the matrix F_n / is unitary (its inverse equals the conjugate of its transpose).

Proof. To establish Eq. (3.3), we must show

Since w = e^{2πi / n}, clearly w = W^–1 and so

To sum this expression over k = 0,... ,n — 1, we use the following identity [see Eq. (1.26)]:

Set z = w^l–j ; note that z_n = 1 because wⁿ = 1. Also note that w^{l–j J} = 1 unless j = l for 0 ≤ j, l ≤ n – 1. The previous equation becomes

Thus

as desired.

Both F_n and are linear transformations from S_n to itself. Here are some additional properties that are analogous to the corresponding properties for the Fourier transform given in Theorems 2.6 and 2.10. The derivations of these properties will be left as exercises.

Theorem 3.4 The following properties hold for the discrete Fourier transform.

• Shifts or translations. If y ∈ S_n and zk = yK+1, then F[z]_j = w^j F[y]j, where w = e^{2πi / n}.
• Convolutions. If y ∈ S_n and z ∈ S_n, then the sequence defined by is also in S_n. The sequence y * z is called the convolution of the sequences y and z.

• The Convolution Theorem: F [y * z] = F[y]kF[z]k.
• If y ε S_n is a sequence of real numbers, then

or

Since n-k = k, only the ₀ … n/2-1 need to be calculated (assuming n is even). Further savings in calculations are possible by using the Fast Fourier Transform algorithm, which we discuss in the next section.

At the end of Section 3.1, we indicated that if y_k = f(k/2πn) for a continuous function, f, then _k is an accurate approximation of the kth Fourier coefficient of f, a_k, when k is small relative to n. The equation _n-k = k further emphasizes this point. For example, if n is large and k = 1, then _n–1 = ₁, which approximates a₁—a low-frequency coefficient. Therefore _n–1 has no resemblance to a_n–1—a high-frequency coefficient.

3.1.3 The Fast Fourier Transform

Since the discrete Fourier transform is such a valuable tool in signal analysis, an efficient algorithm for computing the discrete Fourier and inverse Fourier transforms is very desirable. The computation of the discrete Fourier transform is equivalent to multiplying the sequence y (written as an n × 1 column vector) by the n × n matrix _n, an operation that requires n² multiplications. However, the fast Fourier transform is an algorithm that takes advantage of the special form of the matrix F_n to reduce the number of multiplications to roughly 5n log₂ n. If n = 1000, then this simplification reduces the number of multiplications from a million to about 50,000. The relative savings increase as n gets larger.

The fast Fourier transform algorithm will require an even number of nodes n = 2N. We consider a sequence y = (y₀, y_2N–1) (extended periodically, with period n = 2N). The y>k are calculated via

Splitting the preceding sum into even and odd indices yields

Recall that w = e^{2πi / n} and n = 2N. Let W := exp(2πi / N) = w². We obtain

The preceding equation expresses yît in terms of two discrete Fourier transforms with n replaced by ./V:

for 0 ≤ k ≤ 2N – 1. Further savings are possible. In the last equation, replace k by k + N and use the following facts:

The result is that for 0 ≤ k ≤ N – 1 we have

(3.4)

(3.5)

Thus we have described (F_2Ny)k = _k for 0 ≤ k ≤ 2N – 1 in terms of FN[{y₀ ,y₂ , . . . , y2_N–2}]_k and F _N [{y1 , y3 , . . . , y2_N–₁}]_k for 0 ≤ k ≤ N – 1. Similar formulas can be derived for the inverse discrete Fourier transform; they are

for 0 ≤ k ≤ N – 1 (the factor of 1/2 appears because the inversion formula contains a factor of 1/n).

In matrix terms the computation of the discrete Fourier transform is given by

where FN is the matrix defined in Eq. (3.2), I_N is the N ~x N identity matrix and Dn is the diagonal matrix with the entries 1, w, w²,..., w^{N -1} on the diagonal.

The number of multiplications required to compute the discrete Fourier transform from its definition is about n² = 4N² [since the discrete Fourier transform requires multiplication by the n × n matrix F_n given in Eq. (3.2)]. The number of operations required to compute _k and _k+N in Eqs. (3.4) and (3.5) is 2N² + N + 1 (since FN[{y₀ , y₂, . . ., y_2N–2}]k and FN[{y₁, y₃, . . ., y_2N–1}]k both require N² multiplications). The number of multiplications has almost been cut in half (provided N is large). If N is a multiple of two, then the process can be repeated by writing F _N as a product involving Fn/2. This would cut the number of multiplications by almost another factor of 1/2 or almost one-quarter of the original number.

To keep going, we need to choose the original value of n to be a power of 2, say n = 2^L. Then we can iterate the preceding procedure L times. The Lth iteration involves multiplication by the 1 × 1 matrix which is just the scalar 1. How many multiplications are involved with the full implementation of this algorithm? The answer is found as follows. Let Ki be the number of real multiplications required to compute F[y] by the preceding method for n = 2^L (and so N = 2^L-1). From the formulas derived previously, we see that to compute F[y], we need to compute F [jeven] and F [jodd]. The number of multiplications required is proportional to 2Ki-1 + N = 2Ki-1 + 2^L-1. Thus, Kl is related to A-L-1 via

When L = 0, n = 20 = 1 and no multiplications are required; thus, K0 = 0. Inserting L = 1 in the last equation, we find that K1 ~ 1. Similarly, setting L = 2 then yields K2 ~2x21. Similarly, K₃ ≈ 3 × 2², K4 ~ 4 × 2³, and so on. The general formula is Ki s ≈ L × 2^L-1. Setting n =2^L and noting that L = log₂w, we see that the number of multiplications required is proportional to nlog₂w. The actual proportionality constant is about 5.

Similar algorithms can be obtained when n = N₁N² ... N_m, although the fastest one is obtained in the case discussed previously. For a discussion of this and related topics, consult Briggs and Henson (1995).

Example 3.5 A plot of the absolute value of the discrete Fourier coefficients for the function f(t) = t +1² is given in Figure 3.3 (this graph was produced by MATLAB’s fft routine). The horizontal axis is k = 0,1,2,... ,n, with n = 64, and the dots represent (k, |k|). Since the Fourier coefficients decay with frequency (by the Riemann–Lebesgue Lemma), the |t| decrease with k until k = n/2. The graph of the |_k| is symmetric about the line k = n/2 = 32 since [_n–k| = |_k| (so the values of |_n/2+1|,..., |_n/2+|,...,|_n–1| are the same as _n/2–1,...,₁|).

Example 3.6 We consider the signal in Figure 3.4 which is generated by the function

over the interval 0 ≤ t ≤ 2π. We wish to filter the high-frequency noise from this signal. The signal is first discretized by evaluating it at 2⁸ = 256 equally spaced points on the interval 0 ≤ t ≤ 2π. Then the fast Fourier transform is used to generate the discrete Fourier coefficients yL k = 0,..., 255. Judging from Figure 3.4, the noise has a frequency that is larger than 5 (cycles per 2π interval). Thus we keep only k for 0 ≤ k ≤ 5 and set _k = 0 for 6 ≤ k ≤ 128. By Theorem 3.4, _k = 0 for 128 ≤ k ≤ 250. Applying the inverse fast Fourier transform to the filtered yL we assemble the filtered signal whose graph is given in Figure 3.5.

Example 3.7 We consider the signal given in Figure 3.6 which is generated by the function

over the interval 0 ≤ t ≤ 2π. We wish to compress this signal by taking its fast Fourier transform and then discard the small Fourier coefficients. As in the previous example, we sample the signal at 2⁸ = 256 equally spaced nodes. We apply the fast Fourier transform to generate _k and set 80% of the _k equal to zero (the smallest 80%). Taking the inverse fast Fourier transform of the new _k, we assemble the compressed signal which is graphed in Figure 3.7. Notice the

Gibbs effect since the partial Fourier series is periodic and so its value at 0 and 2π are the same. ¦

The relative error between the original signal y and the compressed signal yc can be computed to be

(see the Appendix C to see how to compute the relative error using Matlab).

3.1.4 The FFT Approximation to the Fourier Transform

A very important application of the FFT is the computation of a numerical approximation to the Fourier transform of a function f defined on an interval a ≤ t ≤ b from samples taken at n nodes, spaced a constant interval T units apart (i.e., at a, a + T, a + 2T,...,a + (n — 1)T = b – T, where T = (b – a)/n). We assume that the function f is zero outside of the interval [a, b], continuous on a ≤ t ≤ b, and that it satisfies f (b) = f (a). Since most signals are not periodic, this last assumption will usually result in a discontinuity in the signal at t = b, thus causing a Gibbs effect when the function is reconstructed from the Fourier transform of the data.

The Fourier transform of f is given by

We want to change variables so that the time interval runs from 0 to 2π rather than a to b. Letting θ = 2π(or t = a + (b — a)9/2n), f is then given by

Note that if we define the function g and the frequency l via

then f(ω_k) has the form

where gt is the kth Fourier series coefficient for g.

The relationship between the discrete Fourier transform (DFT) and Fourier series coefficients is one that we have dealt with before. First, let

Thus the y_j ’s are the samples of f taken at times a + T_j. The DFT approximation to is _k/n. This implies that the DFT for the y_j’s is related to f via

(3.6)

Formula (3.6) accomplishes our goal of approximating the Fourier transform of f at the frequencies u>t by the discrete Fourier transform of the values of f at n equally spaced nodes on a ≤ t ≤ b.

3.1.5 Application: Parameter Identification

As an application of the discrete Fourier transform, we consider the problem of predicting the sway of a building resulting from wind or other external forces. Buildings, especially very tall ones, sway and vibrate in the wind. This is common and is hard to avoid. Too much sway and vibration can cause structural damage, though. How do we find out how a structure behaves in a strong wind? A rough model for the horizontal displacement u due to a wind exerting a force f is

(3.7)

where a, b, and c are positive constants, with b² — 4ac < 0. This is a standard problem in ordinary differential equations. The twist here is that we don’t know what a, b, and c are. Instead, we have to identify them experimentally.

One way to do this is to hit a wall with a jackhammer and record the response. The impulse force exerted by the jackhammer is modeled by f(t) = f₀δ(t), where S(t) is Dirac’s δ-function, and f₀ is the strength of the impulse. Think of the δ-function as a signal with total strength equal to one and lasts for a very short time. Technically speaking, the δ-function is the limit as h → 0 of f_h whose graph is given in Figure 3.8. The corresponding impulse response u has the form [cf. Boyce and DiPrima (1977, Section 6.4)]

(3.8)

where and .

The idea now is to take measurements of the displacements (or sway) at various intervals after the building is hit by the jackhammer. This will give a (discretized) graph of the function û In theory, the frequency ω can be determined from this graph along with the decay constant μ and the amplitude f₀/(ωa). Since f₀ is known (the strength of the jackhammer), the parameter a can be determined from the amplitude f₀/(ωa) and ω. The parameters b and c can then be determined by V4ac — b² b the equations and (two equations and two unknowns).

The trouble with this approach is that the actual data (the measurements of the displacements) will look somewhat different from the theoretical response u in Eq. (3.8). The measured data will contain noise that adds high frequencies to the response. The linear differential equation and the δ-function used in the preceding model are only rough approximations to the corresponding physical system and actual applied force of the jackhammer. Small nonlinearities in the system and an impulse lasting a finite amount of time will add oscillations to the response. Thus, the frequency ω may be difficult to infer directly from the graph of the data.

In spite of all of the additional factors mentioned previously, an approximate value of ω can still be found. This is where the discrete Fourier transform plays a key role. The frequency ω will be the largest frequency component present in the data, and this largest frequency can be found easily by computing the discrete Fourier transform of the data. An approximate value for the decay constant µ can also be determined by analyzing the amplitude of the data. Then approximate values for the parameters a, b, and c can be determined as indicated previously. These ideas are further investigated in exercises 3-6.

3.1.6 Application: Discretizations of Ordinary Differential Equations

As another application of the DFT, we consider the differential equation in (3.7) and take f to be a continuous, 2π -periodic function of t. There is a well-known analical method for finding the unique periodic solution to this equation [cf. Boyce and DiPrima (1977, Section 3.7.2)], provided that f is known for all t. On the other hand, if we only know f at the points t_j = 2πj/n, for some integer n ≤ 1, this method is no longer applicable.

Instead of trying to work directly with the differential equation itself, we will work with a discretized version of it. There a many ways of discretizing. With h = 2π / n, the one that we use here involves the following approximations:

We let t_k = 2πk / n and u_k = u(t_k) for 0 ≤ k ≤ n. Since u is periodic, u_n = u₀. At t = t_k, the preceding difference formulas become

(3.9)

(3.10)

for 1 ≤ k ≤ n – 1. Substituting these values into the differential equation (3.7) and collecting terms yields the following difference equation (see exercise 10):

(3.11)

where f_k = f(2πk / n), β = ch² + bh – 2a, and γ = a – bh.

Suppose that u ε S_n is solution to the difference equation (3.11). Let õ = F[u] and f = F[f], and w = exp(2πi / n). By taking the DFT of both sides and using the first part of Theorem 3.4, we can show

(3.12)

provided that aw^j + γw^j is never 0 (exercise 12). This gives us û We obtain u by applying the inverse DFT to û

Of course, there are also direct methods for finding the Uk [by solving the linear equations in (3.11) for the uk]. However, often information on the frequency components of the solution is desired. In this case, Eq. (3.12) provides the discrete Fourier transform of the solution directly.

3.2 DISCRETE SIGNALS

Sampling invariably produces digitized signals, not continuous ones. The signal becomes a sequence, x = (... x_–2 , x_–1 , x₀ , x₁,...), and the index k in x_k corresponds to the discrete time. The graph of a discrete signal is called a time series. In this section we discuss discrete versions of many of the Fourier analysis topics discussed in Chapter 2.

3.2.1 Time-Invariant, Discrete Linear Filters

This section should be considered as the discrete analog of Section 2.3.1, where continuous time-invariant filters were discussed. We first define the time translation operator T_p on a sequence x by

In words, T_p takes a sequence x and shifts it p units to the right. For example, if x_k = (1 + k²)^–1, then y = T_p(x) has entries y_k = (1 + (k – p)²)^–1. Note that T_p is a linear operator.

To avoid questions concerning convergence of infinite series, we will temporarily assume that our sequences are finite; that is, we consider a sequence x such that x_k = 0 for all |k| > K, where K is some positive number.

Definition 3.8 A linear operator F : x → y that takes a sequence x into another sequence y is time-invariant if F (T_p(x)) = T_p (F(x)).

This condition is simply the discretized version of the time-invariance condition associated with continuous-time filters discussed (Definition 2.13).

Any linear operator F that takes sequences to sequences is completely determined by what it does to the unit sequences, eⁿ, where

Why? A sequence x can be written as the sum x = ∑_n∈Z ^x_n^en, because its kth component is the sum ∑_n∈z ^x_n^en_k = x_k. The linearity of F then implies

(3.13)

So knowing the values of F on eⁿ completely determines F.

Let fⁿ = F (eⁿ) be the sequence that eⁿ is transformed into by F. If F is time-invariant, then

Therefore,

(3.14)

We also have (T_p(fⁿ))_k = from the definition of T_p. Therefore (3.14) becomes . If we set n = 0, then . Replacing p by n then gives

On the other hand, from (3.13) we have

Since f_kⁿ = f_k⁰_-n, we have

Compare this formula to that in Theorem 2.17 for time-invariant operators in the continuous case. There, the operator had the form of a convolution of two continuous functions (see Definition 2.9). If we replace the integral in the continuous convolutions with a sum, we obtain the preceding formula. Let us therefore define the discrete version of the convolution of two sequences.

Definition 3.9 (Discrete Convolution). Given the sequences x and y, the convolution x * y is defined to be

(3.15)

provided that the series involved are absolutely convergent.

Our analysis in this section can be summarized as follows.

Theorem 3.10 If F is a time-invariant linear operator acting on sequences, it has the form of a convolution; namely, there is a sequence f such that

provided that the series involved are absolutely convergent. Conversely, if F(x) = f *x, then F is a discrete, time-invariant linear operator.

We again call such convolution operators discrete filters. The sequence f, which satisfies F(e0) = f and is thus a response to an “impulse” at discrete time 0, is called the impulse response (IR). (Refer to the discussion following Theorem 2.17.) If f has an infinite number of nonvanishing terms, it is called an infinite impulse response (IIR); if it has only a finite number of nonzero terms, it is a finite impulse response (FIR).

3.2.2 Z-Transform and Transfer Functions

In this section we generalize the discrete Fourier transform to infinite sequences in I². The resulting transform is called the Z-transform. There is a close relationship between the Z-transform and the complex form of Fourier series. We will encounter the Z-transform in Chapter 7.

Recall that I² is the space of all (complex-valued) sequences having finite energy; that is, all sequences x = (.. .x_–1,x₀,x₁, ...) with . The inner product of two sequences x and y in I² is given by 〈x,y〉 =

Definition 3.11 The Z-transform of a sequence x = (... x_–1 ,x₀,x₁,...) ∈l² is the function : [–π, π] → C:

Often, one uses z = e^iϕ in the equation defining (·), so becomes

This is where the “Z” in the transform’s name comes from.

The Z-transform generalizes the definition of the discrete Fourier transform for finite sequences. To see this, suppose x = x₀, x₁ ,..., x_n_1 is a finite sequence (so Xj = 0 for j 0 and j > n). Inserting the angle ϕ = 2itk/n into the Z-transform

gives

where . The right side is the definition of x_k - the Kth coefficient of the discrete Fourier transform of x.

Connection with Fourier Series. There is an important relation between the Z-transform and Fourier series. For a given f ∈ L²[–it, re], its Fourier series is

where x_n is the nth Fourier coefficient of f :

From the definition of x~ we have

(3.16)

Thus, a Fourier series expansion is a process that converts a function f ∈ L²[—it, 7t] into a sequence, {x_n} ∈ I²; and the Z-transform takes the sequence {x_n} back to the function f. The following diagram illustrates this relationship.

Isometry Between L² and I ². A further connection to Fourier series is given by Parseval’s Theorem [see Eq. (1.42)]. If x_n and y_n are the Fourier coefficients of

f and g ∈ L²[—tc, tc], then

With the change of variables ϕ → — ϕ, we obtain

or in other words:

where x = (... x_–1, x₀, x1 ,... ) and y = (... y_–1, y₀, y₁,...) We summarize this discussion in the following theorem.

Theorem 3.12 The Z-transform is an isometry between I² and L²[—tc, tc] (preserves the inner products up to a scalar factor of 1/2tc). In particular, if x = (...x_–1 ,x₀,x₁,...) and y = (...y)_–1,y₀,y₁, ...) are sequences in I², then

Convolutions. Now we discuss the relationship between convolution operators and the Z-transform. We have already mentioned that convolution operators resemble their continuous counterparts. The following theorem carries this analogy one step further (compare with Theorem 2.10).

Theorem 3.13 Suppose f = {f_n} and x = {x_n} are sequences in I². Then

The function f{ϕ) is the Z-transform of the sequence f and is also called the transfer function associated to the operator F, where F(x) = f *x. Often, we identify the convolution operator F with its associated sequence f, and therefore we also identify their transfer functions (i.e., F with f).

Proof. Using the definition of the Z-transform, we have

Writing e^–inϕ = e–i(n–k)ϕ we obtain

where we have switched the order of summation. By replacing n – k by the index m in the second sum, this equation becomes

as desired.

Adjoint of Convolution Operators. Recall that the adjoint of an operator F : l² → _l² is the operator F* : l² → l² which is defined by (see Definition 0.29):

The next theorem calculates a formula for the sequence associated to the adjoint of a convolution operator. This theorem will be needed in Chapter 7.

Theorem 3.14 Suppose F is the convolution operator associated with the sequence f. Then F*, is the convolution operator associated to the sequence . The transfer function for F* is F(ϕ).

Proof. From the definition of convolution and I², we have

The last sum on the right (under the conjugate sign) is (f* _* y)(k), where f*(m) = f_–m Therefore, we obtain

Since (F(x), y) = (x, F*(y)}, by definition, we conclude that F*(y) = (f* * y), where as desired.

The second part of this theorem follows from the first because

as desired.

3.3 DISCRETE SIGNALS & MATLAB

It is important to be able to display, convolve, filter, and generally work numerically with a discrete signal x. In Section 3.2 we defined such a signal to be a bi-infinite sequence, and we called the integer variable k the discrete time. An equivalent way to think about x is that it is a function that assigns to k some real (or complex) number x_k.

The graph of x_k versus k is called a time series. Matlab provides several ways of plotting time series, or discrete data. The simplest is the stem plot. We let the discrete signal be

(3.17)

where the first nonzero entry corresponds to k = 0 and the last to k = 5. For values of k larger than 5 or less than 0, x_k = 0. To plot x_k for k = —2 to k = 7, which will include some zeros, we use these commands. (See Figure 3.9.)

X=[0 0 -1 2 3 -2 0 1 0 0];
dtx= -2:7; (discrete time for x)
stem(dtx,x)

The convolution of two discrete-time signals x and y is x * y, which is defined by

(3.18)

As is the case with the continuous-time convolution, x * y = y * x. The convolution is of interest in discrete-time signal processing because of its connection with linear, time-invariant filters. If H is such a filter, then there is a sequence such that HM = h * x; h is called the impulse response (IR) of the filter H. When the IR h has only a finite number of nonzero h_k’s, we say that H has finite impulse response (FIR). Otherwise, it has infinite impulse response (IIR).

In practical situations, we will work only with finite sequences, but these sequences may be arbitrarily large. To compute convolutions of such sequences, we need to discuss indexing of sequences. For a finite sequence x we will let s_x be the starting value of k considered. We thus assume that x_k = 0 if k ≤ s_x. In addition, we let l_x be the last value of k we consider; again, we assume that x_k = 0 if k > l_x. Also, we will let n_x be the length of the stretch between the starting and last x_k’s, so n_x = l_x – s_x + 1. For example, if we only wish to work with nonzero x_k’s in the sequence x defined previously in Eq. (3.17), then s_x = 0, l_x = 5, and n_x = 6.

Our goal is to use Matlab’s function conv to compute x * y when x and y are finite. To do this, we first need to look at what the indexing of the convolution x * y is. That is, we want to know what s_x*y, l_x*_y, and nx*y are, given the corresponding quantities for finite x and y. When y is finite, Eq. (3.18) has the form

(3.19)

where the second equality follows from a change of index. Since 0 ≤ k ≤ l_y – s_y, the index n – k – s_y satisfies

If n – ly > l_x, then x^n–k–s_y = 0 and (x _* y)_n = 0. Similarly, if n – s_y ≤ s_x, we have (x * y)_n = 0. In addition, by direct substitution we see that (x * y)_sy+_sx = x_sxy_s and (x * y)l_y+l_x =xl_xl_y . It follows that the starting index for x * y is s_x*_y = s_x + s_y, that the last index is l_x*y = l_x + l_y, and that the number of terms we need to compute is n_x*y = l_x + l_y – s_x – s_y + 1 = n_x — 1 + n_y – 1 + 1= n_x + n_y – 1. We list these below.

(3.20)

Matlab stores a finite sequence in a row or column vector, depending on the user’s choice. Such vectors are always indexed starting with 1 and going up to the number of entries in the vector. Let’s look at our earlier example, where we let x=[0 0 -1 2 3 -2 0 1 0 0]. To access the fourth entry, we type in x(4). When we hit the enter key, we would get ans = 2. Now, as a sequence, x is a function the discrete time, k. In the entries given above, k runs from k = –2 through k = 7. For plotting purposes and other reasons, we define a discrete time vector to carry this information, dtx=-2 : 7. For example, to find the discrete time that corresponds to the fourth entry in x, we type in dtx(4) and hit enter to get ans=i.

Let’s use this to finding x * y. Again, take x=[0 0 -1 2 3 -2 0 1 0 0] anddtx=-2 : 7. Also, let y=[1 -1 2 4] and dty=8 : 11. To find the convolution of these two and plot the result, we use these commands:

X=[0 0 -1 2 3 -2 0 1 0 0];
cltx= -2:7;
y=[1 -1 2 4];
clty=8 : 11;
z=conv(x,y);
dtz=6:18;
stem(dtz,z)

The result of this is that z=[0 0 -1 3 -1 -5 1 6 9 -9 2 4 0 0]. Apart from these entries, which are for discrete times k = 6 through k = 18, the other entries in z = x * y are all 0. So, for example, z-20 = 0 and z5 = 0. To include more entries on the plot, say for k = -5 to k = 27, one needs to pad the row vector z with zeros. There are many ways to do this. Here is one of them.

dtz=-2:27;
stem(dtz,[zeros(size(-2:5)) z zeros(size(19:27))])

Exercises 17 through 22 below illustrate the concepts presented here. They require Matlab.

EXERCISES

Some of the following problems require the fast Fourier transform (e.g., Maple or Matlab’s FFT routine).

1. Show that T_p{eⁿ) = e^n+p, where T_p and eⁿ are defined in Section 3.2.

2. Prove Theorem 3.4.

3. Derive Eq. (3.8).

4. Plot the solution u in Eq. (3.8), given a = 1, b = 2, c = 37, and /₀ = 1.

5. In the previous exercise, use MATLAB or some similar program to take 256 samples of u on the interval [0,4]. Plot the absolute value of the FFT of û Locate the largest natural frequency u>/2n. Compare with the result from the previous exercise.

6. The FFT locates the frequency reasonably well. Find a way to estimate the damping constant, /x, from FFT data.

7. Filtering. Let

Discretize f by setting y_k = f(2kπ/256), & = 1...256. Use the fast Fourier transform to compute yit for 0 ≤ k ≤ 256. Recall from Theorem 3.4 that _n-k = _k. Thus, the low-frequency coefficients are y₀ • • • ym and y256-m ... Y256 for some low value of m. Filter out the high-frequency terms by setting y_k = 0 for m ≤ k ≤ 255 – m with m = 6; then apply the inverse fast Fourier transform to this new set of y_k to compute the y_k (now filtered); plot the new values of y_k and compare with the original function. Experiment with other values of m.

8. Compression. Let tol = 1.0. In exercise 7, if |_k| ≤ tol then set yπ equal to zero. Apply the inverse fast Fourier transform to this new set of _k to compute the _k; plot the new values of _k and compare with the original function. Experiment with other values of tol. Keep track of the percentage of Fourier coefficients which have been filtered out. Matlab’s sort command is useful for finding a value for tol in order to filter-out a specified percentage of coefficients (see the compression routine in the section on Matlab code). Compute the relative I² error of the compressed signal as compared with the original signal (again, see the section on Matlab code).

9. Repeat the previous two exercises over the interval 0 ≤ x ≤ 1 with the function

where

This time, _k = f(k/256).

10. Derive Eq. (3.11) by inserting Eq. (3.9) and Eq. (3.10) into Eq. (3.7) at t = tk = 2nk/n.

11. Let a, b and c be strictly positive. As in Eq. (3.11), set ft = ch² + bh – 2a, and y = a – bh. Show that a w^j+ β + γ is never 0. (Hint: Show that the quadratic equation az² + ftz + y = 0 has no solutions for which |z| = 1.)

12. Derive Eq. (3.12), assuming that a wⁱ + β + γ is never 0.

13. Find the exact, π/3-periodic solution to u″ + 2u′ + 2u = 3 cos(6t). Compare it with the discrete approximation obtained for these values of n : n = 16, n = 64, and n = 256. You will need to use MATLAB or similar software to compute the FFT and inverse FFT involved.

14. Recall that the complex exponentials exp(mf ) are 2π-periodic eigenfunctions of the operator D²[u] = u″. A discretized version of this operator acts on the periodic sequences in S_n. If u_k is n-periodic, then define L[u] via

(a) Show that L maps S_n into itself.

(b) For n = 4, show the matrix M₄ that represents L is

(c) Observe that M₄ is self-adjoint and thus diagonalizable. Find its eigenvalues and corresponding eigenvectors. How are these related to the matrix columns for the matrix F₄ in Eq. (3.2)? Could you have diagonalized this matrix via an FFT? Explain.

(d) Generalize this result for all n. (Hint: Use the DFT on L[u]_k; recall that the FFT is really a fast DFT.)

15. (Circulants and the DFT.) An n × n matrix A is called a circulant if all of its diagonals (main, super, and sub) are constant and the indices are interpreted “modulo” n. For example, this 4 × 4 matrix is a circulant:

(a) Look at the n-periodic sequence a, where a_l = A_l+1,1, l = 0,..., n – 1. Write the entries of A in terms of the sequence a

(b) Let X be an n × 1 column vector. Show that Y = AX is equivalent to y = a _* x if x, y are the n-periodic sequences for which x_l = x_l+1,1 and similarly for y_l = Y_l+1,1 l = 0, ..., n – 1.

(c) Prove that the DFT diagonalizes all circulant matrices; that is, n^–1 A = D, where D is diagonal. What are the diagonal entries of D (i.e., what are the eigenvalues of A)?

16. Find the Z-transform for the sequence

The remaining exercises require Matlab.

17. Do each of the sets of commands listed in Section 3.3. Print the resulting plots.

18. Find the convolution of the x and y. Here, x_k = 0 for k < 3 and k > –4, and X_k = 1 for –4 ≤ k ≤ 3. For y_k, assume that y_k = 0 for k < –2 and for y_k < –8. When –8 ≤ k ≤ –2, y_k = k + 2. Find x _* y and determine the discrete time index. Plot the result with stem, again using the correct time index. Put a title on your plot by using the command title(’Convolution of x _* y’). Print the result.

19. Take t=linspace (0,2*pi,20) and x=sin(t). Do the plots stem(t,x), stem(t,x,’:r’ , ’fill’), and stem(t,x,’ - .sb’ , ’fill’). Put them all in one plot with the commands below and print the result.

subplot(1,3,1), stem(t,x)
title(’Default Stem Plot’)
subplot(1,3,2), stem(t,x,’:r’,’fill’)
title(’Filled Stem Plot’)
subplot(1,3,3), stem(t,x,’ sk’)
title(’Square Marker Stem Plot’)

20. This exercise illustrates the use of another plotting tool, stairs. Start with the following commands.

t=linspace(-pi, pi, 20); x=sin(t);
stairs(t,x)

Next, change the plot by using a “dash–dot” line instead of a solid one. We will also change the color to red: stairs(t,x,’-.r’). We will now combine stairs, stem and plot. Title the plot and print the result.

t=linspace(-pi,pi,20); x=sin(t);
tt=linspace(-pi,pi,600); xx=sin(tt);
stairs(t,x,’-.r’), hold on
stem(t,x,’:sb’,’fill’) (Dotted stems & filled circles)
plot(tt,xx,’k’), hold off

21. The stairs plots are well suited to doing plots involving the Haar scaling function and wavelet. The Haar scaling function is defined by

Use stairs to plot f(x) = 3φ (x + 1) – 2φ(x) + φ(x – 1) + φ (x – 2) on the interval –3 ≤ x ≤ 4. On the same interval, plot g(x) = φ(2x + 3) – φ(2x + 2) + 2φ(2x) – φ(2x – 3). For g, use a dash-dot pattern (-.) and make the color red. (Hint: To plot two functions, you will need to use hold on.)

22. This problem pertains to the Z-transform. Let h be a finite discrete-time signal. If we let z = e^iw, then the Z-transform is

We want to illustrate how to numerically compute and plot h. For simplicity, we will work with an h for which S_h = –0n, l = n, and h_k = 1/n for – n ≤ k ≤ n. We will do the following.

w=linspace(-pi,pi,600); z=exp(i*w); (Initialize Z-)
n=2; h=ones(1,2*n+1)/(2*n+1);
H=z."n.*polyval(h,1./z); (Compute the Z-transform.)
norm(imag(H)) (This should be small; if not, H is wrong.)
plot(w,real(H))

In addition to n = 2, do this for n = 4, and n = 5. Title, and then print, each plot.