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 Yj := F(hj) = F(2πj / n), j = 0,..., n. If F(t) is 2π -periodic, then Y0 = Yn and the preceding formula becomes
Applying this formula to the computation of the kth complex Fourier coefficient gives
Therefore
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 ak for k ≤ n, because the at are not n-periodic. In fact, this expression only approximates ak 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 Sn be the set of n-periodic sequences of complex numbers. Each element, in Sn, can be thought of as a periodic discrete signal where yj is the value of the signal at a time node t = tj. The sequence yj is n -periodic if yk+n = yk for any integer k. The set Sn 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 Xj + yj and the jth component of c{x} is cxj . 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 (Fn{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 ak, 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 = (y0,..., yn–1)T and and where
Note that Fn is a symmetric matrix.
3.1.2 Properties of the Discrete Fourier Transform
The discrete Fourier transform of a signal in Sn produces a discrete signal that also belongs to Sn as the following lemma shows.
Lemma 3.2 Fn is a linear operator from Sn to Sn.
Proof. To show that Fn{y} belongs to Sn we must show that Fn{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, Fn 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 = {yk} is an element of Sn. Let Fn(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 Fn [F is defined in Eq. (3.2)], that is,
Since = ()(y), the equation y = Fn() is equivalent to
Since Fn is symmetric, this equation means that the matrix Fn / is unitary (its inverse equals the conjugate of its transpose).
Proof. To establish Eq. (3.3), we must show
Since w = e2π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 = wl–j ; note that zn = 1 because wn = 1. Also note that wl–j J = 1 unless j = l for 0 ≤ j, l ≤ n – 1. The previous equation becomes
Thus
as desired.
Both Fn and are linear transformations from Sn 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 ∈ Sn and zk = yK+1, then F[z]j = wj F[y]j, where w = e2πi / n.
- Convolutions. If y ∈ Sn and z ∈ Sn, then the sequence defined by is also in Sn. 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 ε Sn is a sequence of real numbers, then
or
Since n-k = k, only the 0 … 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 yk = f(k/2πn) for a continuous function, f, then k is an accurate approximation of the kth Fourier coefficient of f, ak, 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 = 1, which approximates a1—a low-frequency coefficient. Therefore n–1 has no resemblance to an–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 n2 multiplications. However, the fast Fourier transform is an algorithm that takes advantage of the special form of the matrix Fn to reduce the number of multiplications to roughly 5n log2 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 = (y0, y2N–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 = e2πi / n and n = 2N. Let W := exp(2πi / N) = w2. 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
Thus we have described (F2Ny)k = k for 0 ≤ k ≤ 2N – 1 in terms of FN[{y0 ,y2 , . . . , y2N–2}]k and F N [{y1 , y3 , . . . , y2N–1}]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), IN is the N ~x N identity matrix and Dn is the diagonal matrix with the entries 1, w, w2,..., wN -1 on the diagonal.
The number of multiplications required to compute the discrete Fourier transform from its definition is about n2 = 4N2 [since the discrete Fourier transform requires multiplication by the n × n matrix Fn given in Eq. (3.2)]. The number of operations required to compute k and k+N in Eqs. (3.4) and (3.5) is 2N2 + N + 1 (since FN[{y0 , y2, . . ., y2N–2}]k and FN[{y1, y3, . . ., y2N–1}]k both require N2 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 = 2L. 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 = 2L (and so N = 2L-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 + 2L-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, K3 ≈ 3 × 22, K4 ~ 4 × 23, and so on. The general formula is Ki s ≈ L × 2L-1. Setting n =2L and noting that L = log2w, we see that the number of multiplications required is proportional to nlog2w. The actual proportionality constant is about 5.
Similar algorithms can be obtained when n = N1N2 ... Nm, 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 +12 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,...,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 28 = 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 28 = 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 yj ’s are the samples of f taken at times a + Tj. The DFT approximation to is k/n. This implies that the DFT for the yj’s is related to f via
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
where a, b, and c are positive constants, with b2 — 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) = f0δ(t), where S(t) is Dirac’s δ-function, and f0 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 fh whose graph is given in Figure 3.8. The corresponding impulse response u has the form [cf. Boyce and DiPrima (1977, Section 6.4)]
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 f0/(ωa). Since f0 is known (the strength of the jackhammer), the parameter a can be determined from the amplitude f0/(ωa) and ω. The parameters b and c can then be determined by V4ac — b2 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 tj = 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 tk = 2πk / n and uk = u(tk) for 0 ≤ k ≤ n. Since u is periodic, un = u0. At t = tk, the preceding difference formulas become
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):
where fk = f(2πk / n), β = ch2 + bh – 2a, and γ = a – bh.
Suppose that u ε Sn 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
provided that awj + γwj 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.
Sampling invariably produces digitized signals, not continuous ones. The signal becomes a sequence, x = (... x–2 , x–1 , x0 , x1,...), and the index k in xk 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 Tp on a sequence x by
In words, Tp takes a sequence x and shifts it p units to the right. For example, if xk = (1 + k2)–1, then y = Tp(x) has entries yk = (1 + (k – p)2)–1. Note that Tp 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 xk = 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 (Tp(x)) = Tp (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, en, where
Why? A sequence x can be written as the sum x = ∑n∈Z xnen, because its kth component is the sum ∑n∈z xnenk = xk. The linearity of F then implies
So knowing the values of F on en completely determines F.
Let fn = F (en) be the sequence that en is transformed into by F. If F is time-invariant, then
Therefore,
We also have (Tp(fn))k = from the definition of Tp. 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 fkn = fk0-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 I2. 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 I2 is the space of all (complex-valued) sequences having finite energy; that is, all sequences x = (.. .x–1,x0,x1, ...) with . The inner product of two sequences x and y in I2 is given by 〈x,y〉 =
Definition 3.11 The Z-transform of a sequence x = (... x–1 ,x0,x1,...) ∈l2 is the function : [–π, π] → C:
Often, one uses z = eiϕ 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 = x0, x1 ,..., xn_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 xk - 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 ∈ L2[–it, re], its Fourier series is
where xn 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 ∈ L2[—it, 7t] into a sequence, {xn} ∈ I2; and the Z-transform takes the sequence {xn} back to the function f. The following diagram illustrates this relationship.
Isometry Between L2 and I 2. A further connection to Fourier series is given by Parseval’s Theorem [see Eq. (1.42)]. If xn and yn are the Fourier coefficients of
f and g ∈ L2[—tc, tc], then
With the change of variables ϕ → — ϕ, we obtain
or in other words:
where x = (... x–1, x0, x1 ,... ) and y = (... y–1, y0, y1,...) We summarize this discussion in the following theorem.
Theorem 3.12 The Z-transform is an isometry between I2 and L2[—tc, tc] (preserves the inner products up to a scalar factor of 1/2tc). In particular, if x = (...x–1 ,x0,x1,...) and y = (...y)–1,y0,y1, ...) are sequences in I2, 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 = {fn} and x = {xn} are sequences in I2. 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 : l2 → l2 is the operator F* : l2 → l2 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 I2, 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.
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 xk.
The graph of xk 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
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, xk = 0. To plot xk 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
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 hk’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 sx be the starting value of k considered. We thus assume that xk = 0 if k ≤ sx. In addition, we let lx be the last value of k we consider; again, we assume that xk = 0 if k > lx. Also, we will let nx be the length of the stretch between the starting and last xk’s, so nx = lx – sx + 1. For example, if we only wish to work with nonzero xk’s in the sequence x defined previously in Eq. (3.17), then sx = 0, lx = 5, and nx = 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 sx*y, lx*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 ≤ ly – sy, the index n – k – sy satisfies
If n – ly > lx, then xn–k–sy = 0 and (x * y)n = 0. Similarly, if n – sy ≤ sx, we have (x * y)n = 0. In addition, by direct substitution we see that (x * y)sy+sx = xsxys and (x * y)ly+lx =xlxly . It follows that the starting index for x * y is sx*y = sx + sy, that the last index is lx*y = lx + ly, and that the number of terms we need to compute is nx*y = lx + ly – sx – sy + 1 = nx — 1 + ny – 1 + 1= nx + ny – 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.
Some of the following problems require the fast Fourier transform (e.g., Maple or Matlab’s FFT routine).
1. Show that Tp{en) = en+p, where Tp and en 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 /0 = 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 yk = 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 y0 • • • ym and y256-m ... Y256 for some low value of m. Filter out the high-frequency terms by setting yk = 0 for m ≤ k ≤ 255 – m with m = 6; then apply the inverse fast Fourier transform to this new set of yk to compute the yk (now filtered); plot the new values of yk 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 I2 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 = ch2 + bh – 2a, and y = a – bh. Show that a wj+ β + γ is never 0. (Hint: Show that the quadratic equation az2 + ftz + y = 0 has no solutions for which |z| = 1.)
12. Derive Eq. (3.12), assuming that a wi + β + γ 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 D2[u] = u″. A discretized version of this operator acts on the periodic sequences in Sn. If uk is n-periodic, then define L[u] via
(a) Show that L maps Sn into itself.
(b) For n = 4, show the matrix M4 that represents L is
(c) Observe that M4 is self-adjoint and thus diagonalizable. Find its eigenvalues and corresponding eigenvectors. How are these related to the matrix columns for the matrix F4 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 al = Al+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 xl = xl+1,1 and similarly for yl = Yl+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, xk = 0 for k < 3 and k > –4, and Xk = 1 for –4 ≤ k ≤ 3. For yk, assume that yk = 0 for k < –2 and for yk < –8. When –8 ≤ k ≤ –2, yk = 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 = eiw, 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 Sh = –0n, l = n, and hk = 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.