8.5 NUMERICAL [COMPUTATION OF THE](#page-13-0) FOURIER TRANSFORM: THE

← Back to LINEAR SYSTEMS AND SIGNALS Overview

8.5 NUMERICAL COMPUTATION OF THE FOURIER TRANSFORM: THE DISCRETE FOURIER TRANSFORM

Numerical computation of the Fourier transform of x(t) requires sample values of x(t) because a digital computer can work only with discrete data (sequence of numbers). Moreover, a computer can compute X(ω) only at some discrete values of ω [samples of X(ω)]. We therefore need to relate the samples of X(ω) to samples of x(t). This task can be accomplished by using the results of the two sampling theorems developed in Secs. 8.1 and 8.4.

We begin with a timelimited signal x(t) (Fig. 8.16a) and its spectrum X(ω) (Fig. 8.16b). Since x(t) is timelimited, X(ω) is nonbandlimited. For convenience, we shall show all spectra as functions of the frequency variable f (in hertz) rather than ω. According to the sampling theorem, the spectrum X(ω) of the sampled signal x(t) consists of X(ω) repeating every fs Hz, where fs = 1/T, as depicted in Fig. 8.16d.† In the next step, the sampled signal in Fig. 8.16c is repeated periodically every T⁰ seconds, as illustrated in Fig. 8.16e. According to the spectral sampling theorem, such an operation results in sampling the spectrum at a rate of T⁰ samples/Hz. This sampling rate means that the samples are spaced at f⁰ = 1/T⁰ Hz, as depicted in Fig. 8.16f.

The foregoing discussion shows that when a signal x(t) is sampled and then periodically repeated, the corresponding spectrum is also sampled and periodically repeated. Our goal is to relate the samples of x(t) to the samples of X(ω).

NUMBER OF SAMPLES

One interesting observation from Figs. 8.16e and 8.16f is that N0, the number of samples of the signal in Fig. 8.16e in one period T0, is identical to N ⁰, the number of samples of the spectrum in Fig. 8.16f in one period fs. To see this, we notice that

$$N_0 = \frac{T_0}{T} \quad N'_0 = \frac{f_s}{f_0} \quad f_s = \frac{1}{T} \quad \text{and} \quad f_0 = \frac{1}{T_0}$$

(8.10)

Using these relations, we see that

$$N_0 = \frac{T_0}{T} = \frac{f_s}{f_0} = N'_0$$

ALIASING AND LEAKAGE IN NUMERICAL COMPUTATION

Figure 8.16f shows the presence of aliasing in the samples of the spectrum X(ω). This aliasing error can be reduced as much as desired by increasing the sampling frequency fs (decreasing the sampling interval T = 1/fs). The aliasing can never be eliminated for timelimited x(t), however, because its spectrum X(ω) is nonbandlimited. Had we started with a signal having a bandlimited spectrum X(ω), there would be no aliasing in the spectrum in Fig. 8.16f. Unfortunately, such a signal is nontimelimited, and its repetition (in Fig. 8.16e) would result in signal overlapping (aliasing in the time domain). In this case, we shall have to contend with errors in signal

^† There is a multiplying constant 1/T for the spectrum in Fig. 8.16d [see Eq. (8.2)], but this is irrelevant to our discussion here.

samples. In other words, in computing the direct or inverse Fourier transform numerically, we can reduce the error as much as we wish, but the error can never be eliminated. This is true of numerical computation of the direct and inverse Fourier transforms, regardless of the method used. For example, if we determine the Fourier transform by direct integration numerically, by using Eq. (7.9), there will be an error because the interval of integration t can never be made zero. Similar remarks apply to numerical computation of the inverse transform. Therefore, we should always keep in mind the nature of this error in our results. In our discussion (Fig. 8.16), we assumed x(t) to be a timelimited signal. If x(t) were not timelimited, we would need to timelimit it because numerical computations can work only with finite data. Furthermore, this data truncation causes error because of spectral spreading (smearing) and leakage, as discussed in Sec. 7.8. The leakage also causes aliasing. Leakage can be reduced by using a tapered window for signal truncation. But this choice increases spectral spreading or smearing. Spectral spreading can be reduced by increasing the window width (i.e., more data), which increases T0, and reduces f⁰ (increases spectral or frequency resolution).

PICKET FENCE EFFECT

The numerical computation method yields only the uniform sample values of X(ω). The major peaks or valleys of X(ω) can lie between two samples and may remain hidden, giving a false picture of reality. Viewing samples is like viewing the signal and its spectrum through a “picket fence” with upright posts that are very wide and placed close together. What is hidden behind the pickets is much more than what we can see. Such misleading results can be avoided by using a sufficiently large N0, the number of samples, to increase resolution. We can also use zero padding (discussed later) or the spectral interpolation formula [Eq. (8.9)] to determine the values of X(ω) between samples.

POINTS OF DISCONTINUITY

If x(t) or X(ω) has a jump discontinuity at a sampling point, the sample value should be taken as the average of the values on the two sides of the discontinuity because the Fourier representation at a point of discontinuity converges to the average value.

DERIVATION OF THE DISCRETE FOURIER TRANSFORM (DFT)

If x(nT) and X(rω0) are the nth and rth samples of x(t) and X(ω), respectively, then we define new variables xn and Xr as

$$x_n = Tx(nT) = \frac{T_0}{N_0}x(nT)$$

\n(8.11)

and

$$X_r = X(r\omega_0)$$

where

$$\omega_0 = 2\pi f_0 = \frac{2\pi}{T_0}$$

We shall now show that xn and Xr are related by the following equations† :

$$X_r = \sum_{n=0}^{N_0 - 1} x_n e^{-jr\Omega_0 n}$$

(8.12)

and

$$x_n = \frac{1}{N_0} \sum_{r=0}^{N_0 - 1} X_r e^{ir\Omega_0 n}$$

\n(8.13)

where

$$\Omega_0 = \omega_0 T = \frac{2\pi}{N_0}$$

These equations define the direct and the inverse discrete Fourier transforms, with Xr the direct discrete Fourier transform (DFT) of xn, and xn the inverse discrete Fourier transform (IDFT) of Xr. The notation

$$x_n \Longleftrightarrow X_r$$

is also used to indicate that xn and Xr are a DFT pair. Remember that xn is T0/N⁰ times the nth sample of x(t) and Xr is the rth sample of X(ω). Knowing the sample values of x(t), we can use the DFT to compute the sample values of X(ω)—and vice versa. Note, however, that xn is a function of n (n = 0, 1, 2,…,N⁰ − 1) rather than of t and that Xr is a function of r (r = 0, 1, 2,…,N⁰ − 1) rather than of ω. Moreover, both xn and Xr are periodic sequences of period N⁰ (Figs. 8.16e, 8.16f). Such sequences are called N0*-periodic sequences*. The proof of the DFT relationships in Eqs. (8.12) and (8.13) follows directly from the results of the sampling theorem. The sampled signal x(t) (Fig. 8.16c) can be expressed as

$$\overline{x}(t) = \sum_{n=0}^{N_0 - 1} x(nT)\delta(t - nT)$$

Since δ (t −nT) ⇐⇒ e−jnω*^T* , applying the Fourier transform yields

$$\overline{X}(\omega) = \sum_{n=0}^{N_0 - 1} x(nT) e^{-jn\omega T}$$

But from Fig. 8.1f [or Eq. (8.2)], it is clear that over the interval |ω| ≤ ωs/2, X(ω), the Fourier transform of x(t) is X(ω)/T, assuming negligible aliasing. Hence,

$$X(\omega) = T\overline{X}(\omega) = T\sum_{n=0}^{N_0 - 1} x(nT)e^{-jn\omega T} \qquad |\omega| \le \frac{\omega_s}{2}$$

^† In Eqs. (8.12) and (8.13), the summation is performed from 0 to ^{N0 −}1. It is shown in Sec. 9.1-2 [Eqs. (9.6) and (9.7)] that the summation may be performed over any successive N⁰ values of n or r.

and

$$X_r = X(r\omega_0) = T \sum_{n=0}^{N_0 - 1} x(nT)e^{-nkr\omega_0 T}$$

\n(8.14)

If we let ω0T = 0, then from Eq. (8.10),

$$\Omega_0 = \omega_0 T = 2\pi f_0 T = \frac{2\pi}{N_0}$$

Also, from Eq. (8.11),

$$Tx(nT) = x_n$$

Therefore, Eq. (8.14) becomes

$$X_r = \sum_{n=0}^{N_0 - 1} x_n e^{-jr\Omega_0 n} \qquad \Omega_0 = \frac{2\pi}{N_0}$$

This is Eq. (8.12), which we set to prove.

The inverse transform relationship of Eq. (8.13) can be derived by using a similar procedure with the roles of t and ω reversed, but here we shall use a more direct proof. To prove Eq. (8.13), we multiply both sides of Eq. (8.12) by ejm0*^r* and sum over r as

$$\sum_{r=0}^{N_0-1} X_r e^{jm\Omega_0 r} = \sum_{r=0}^{N_0-1} \left[ \sum_{n=0}^{N_0-1} x_n e^{-jr\Omega_0 n} \right] e^{jm\Omega_0 r}$$

By interchanging the order of summation on the right-hand side, we have

$$\sum_{r=0}^{N_0-1} X_r e^{jm\Omega_0 r} = \sum_{n=0}^{N_0-1} x_n \left[ \sum_{r=0}^{N_0-1} e^{j(m-n)\Omega_0 r} \right]$$

As the footnote below readily shows, the inner sum on the right-hand side is zero for n = m and is N⁰ when n = m. † Thus, the outer sum will have only one nonzero term when n = m, and it is

† We show that ^N

$$\sum_{n=0}^{N_0-1} e^{jk\Omega_0 n} = \begin{cases} N_0 & k = 0, \pm N_0, \pm 2N_0, \dots \\ 0 & \text{otherwise} \end{cases}$$

(8.15)

Recall that 0N^{0 = 2}π. So ejk0*ⁿ* ⁼ 1 when ^{k =} 0,±N0,±2N0,… Hence, the sum on the left-hand side of Eq. (8.15) is N0. To compute the sum for other values of k, we note that the sum on the left-hand side of Eq. (8.15) is a geometric progression with common ratio ^{α = e}jk⁰ . Therefore, (see Sec. B.8-3)

$$\sum_{n=0}^{N_0-1} e^{jk\Omega_0 n} = \frac{e^{jk\Omega_0 N_0} - 1}{e^{jk\Omega_0} - 1} = 0 \qquad (e^{jk\Omega_0 N_0} = e^{j2\pi m} = 1)$$

810 CHAPTER 8 SAMPLING: THE BRIDGE FROM CONTINUOUS TO DISCRETE

N0xn = N0xm. Therefore,

$$x_m = \frac{1}{N_0} \sum_{r=0}^{N_0 - 1} X_r e^{im\Omega_0 r} \qquad \Omega_0 = \frac{2\pi}{N_0}$$

Because Xr is N⁰ periodic, we need to determine the values of Xr over any one period. It is customary to determine Xr over the range (0, N⁰ − 1), rather than over the range (−N0/2,(N0/2)−1).†

CHOICE OF T AND T⁰

In DFT computation, we first need to select suitable values for N⁰ and T or T0. For this purpose, we begin by deciding on B, the essential bandwidth (in hertz) of the signal. The sampling frequency fs must be at least 2B, that is,

$$\frac{f_s}{2} \ge B$$

Moreover, the sampling interval T = 1/fs [Eq. (8.10)], and

$$T \le \frac{1}{2B} \tag{8.16}$$

Once we pick B, we can choose T according to Eq. (8.16). Also,

$$f_0 = \frac{1}{T_0} \tag{8.17}$$

where f⁰ is the frequency resolution [separation between samples of X(ω)]. Hence, if f⁰ is given, we can pick T⁰ according to Eq. (8.17). Knowing T⁰ and T, we determine N⁰ from

$$N_0 = \frac{T_0}{T}$$

ZERO PADDING

Recall that observing Xr is like observing the spectrum X(ω) through a picket fence. If the frequency sampling interval f⁰ is not sufficiently small, we could miss out on some significant details and obtain a misleading picture. To obtain a higher number of samples, we need to reduce f0. Because f⁰ =1/T0, a higher number of samples requires us to increase the value of T0, the period of repetition for x(t). This option increases N0, the number of samples of x(t), by adding dummy samples of 0 value. This addition of dummy samples is known as zero padding. Thus, zero padding increases the number of samples and may help in getting a better idea of the spectrum X(ω) from its samples Xr. To continue with our picket fence analogy, zero padding is like using more, and narrower, pickets.

^† The DFT of Eq. (8.12) and the IDFT of Eq. (8.13) represent a transform in their own right, and they are exact. There is no approximation. However, xn and Xr, thus obtained, are only approximations to the actual samples of a signal x(t) and of its Fourier transform X(ω).

ZERO PADDING DOES NOT IMPROVE ACCURACY OR RESOLUTION

Actually, we are not observing X(ω) through a picket fence. We are observing a distorted version of X(ω) resulting from the truncation of x(t). Hence, we should keep in mind that even if the fence were transparent, we would see a reality distorted by aliasing. Seeing through the picket fence just gives us an imperfect view of the imperfectly represented reality. Zero padding only allows us to look at more samples of that imperfect reality. It can never reduce the imperfection in what is behind the fence. The imperfection, which is caused by aliasing, can be lessened only by reducing the sampling interval T. Observe that reducing T also increases N0, the number of samples, and is like increasing the number of pickets while reducing their width. But in this case, the reality behind the fence is also better dressed and we see more of it.

EXAMPLE 8.7 Number of Samples and Frequency Resolution

A signal x(t) has a duration of 2 ms and an essential bandwidth of 10 kHz. It is desirable to have a frequency resolution of 100 Hz in the DFT (f⁰ = 100). Determine N0.

To have f⁰ = 100 Hz, the effective signal duration T⁰ must be

$$T_0 = \frac{1}{f_0} = \frac{1}{100} = 10$$

ms

Since the signal duration is only 2 ms, we need zero padding over 8 ms. Also, B = 10,000. Hence, fs = 2B = 20,000 and T = 1/fs = 50 µs. Furthermore,

$$N_0 = \frac{f_s}{f_0} = \frac{20,000}{100} = 200$$

The fast Fourier transform (FFT) algorithm (discussed later; see Sec. 8.6) is used to compute DFT, where it proves convenient (although not necessary) to select N⁰ as a power of 2; that is, N⁰ = 2*ⁿ* (n, integer). Let us choose N⁰ = 256. Increasing N⁰ from 200 to 256 can be used to reduce aliasing error (by reducing T), to improve resolution (by increasing T⁰ using zero padding), or a combination of both.

Reducing Aliasing Error. We maintain the same T⁰ so that f⁰ = 100. Hence,

$$f_s = N_0 f_0 = 256 \times 100 = 25,600$$

and $T = \frac{1}{f_s} = 39 \,\mu s$

Thus, increasing N⁰ from 200 to 256 permits us to reduce the sampling interval T from 50 µs to 39 µs while maintaining the same frequency resolution (f⁰ = 100).

Improving Resolution. Here, we maintain the same T = 50 µs, which yields

$$T_0 = N_0 T = 256(50 \times 10^{-6}) = 12.8$$

ms and $f_0 = \frac{1}{T_0} = 78.125$ Hz

812 CHAPTER 8 SAMPLING: THE BRIDGE FROM CONTINUOUS TO DISCRETE

Thus, increasing N⁰ from 200 to 256 can improve the frequency resolution from 100 to 78.125 Hz while maintaining the same aliasing error (T = 50 µs).

Combination of Reducing Aliasing Error and Improving Resolution. To simultaneously reduce alias error and improve resolution, we could choose T = 45 µs and T⁰ = 11.5 ms so that f⁰ = 86.96 Hz. Many other combinations exist as well.

EXAMPLE 8.8 DFT to Compute the Fourier Transform of an Exponential

Use the DFT to compute (samples of) the Fourier transform of e−2*^t u*(t). Plot the resulting Fourier spectra.

We first determine T and T0. The Fourier transform of e−2*^t u*(t) is 1/(jω + 2). This lowpass signal is not bandlimited. In Sec. 7.6, we used the energy criterion to compute the essential bandwidth of a signal. Here, we shall present a simpler, but workable alternative to the energy criterion. The essential bandwidth of a signal will be taken as the frequency at which |X(ω)| drops to 1% of its peak value (see the footnote on page 736). In this case, the peak value occurs at ω = 0, where |X(0)| = 0.5. Observe that

$$|X(\omega)| = \frac{1}{\sqrt{\omega^2 + 4}} \approx \frac{1}{\omega} \qquad \omega \gg 2$$

Also, 1% of the peak value is 0.01 × 0.5 = 0.005. Hence, the essential bandwidth B is at ω = 2πB, where

$$|X(\omega)| \approx \frac{1}{2\pi B} = 0.005 \quad \Rightarrow \quad B = \frac{100}{\pi} \text{ Hz}$$

and from Eq. (8.16),

$$T \le \frac{1}{2B} = \frac{\pi}{200} = 0.015708$$

Had we used 1% energy criterion to determine the essential bandwidth, following the procedure in Ex. 7.20, we would have obtained B = 20.26 Hz, which is somewhat smaller than the value just obtained by using the 1% amplitude criterion.

The second issue is to determine T0. Because the signal is not timelimited, we have to truncate it at T⁰ such that x(T0) 1. A reasonable choice would be T⁰ = 4 because x(4) = e−⁸ = 0.0003351. The result is N⁰ = T0/T = 254.6, which is not a power of 2. Hence, we choose T⁰ = 4, and T = 0.015625 = 1/64, yielding N⁰ = 256, which is a power of 2.

Note that there is a great deal of flexibility in determining T and T0, depending on the accuracy desired and the computational capacity available. We could just as well have chosen T = 0.03125, yielding N⁰ = 128, although this choice would have given a slightly higher aliasing error.

Because the signal has a jump discontinuity at t = 0, the first sample (at t = 0) is 0.5, the averages of the values on the two sides of the discontinuity. We compute Xr (the DFT) from the samples of e−2*^t u*(t) according to Eq. (8.12). Note that Xr is the rth sample of X(ω), and these samples are spaced at f⁰ = 1/T⁰ = 0.25 Hz (ω⁰ = π/2 rad/s).

Because Xr is N⁰ periodic, Xr = X(r+256) so that X²⁵⁶ = X0. Hence, we need to plot Xr over the range r = 0 to 255 (not 256). Moreover, because of this periodicity, X−^r = X(−r+256), and the values of Xr over the range r = −127 to −1 are identical to those over the range r = 129 to 255. Thus, X−¹²⁷ = X129, X−¹²⁶ = X130,…,X−¹ = X255. In addition, because of the property of conjugate symmetry of the Fourier transform, X−^r = X^{∗ r} , it follows that X−¹ = X^{∗ 1} , X−² = X^{∗ 2} ,…,X−¹²⁸ = X^{∗ 128}. Thus, we need Xr only over the range r = 0 to N0/2 (128 in this case).

Figure 8.17 shows the computed plots of |Xr| and Xr. The exact spectra are depicted by continuous curves for comparison. Note the nearly perfect agreement between the two sets of spectra. We have depicted the plot of only the first 28 points rather than all 128 points, which would have made the figure very crowded, resulting in loss of clarity. The points are at the intervals of 1/T⁰ = 1/4 Hz or ω⁰ = 1.5708 rad/s. The 28 samples, therefore, exhibit the plots over the range ω = 0 to ω = 28(1.5708) ≈ 44 rad/s or 7 Hz.

Figure 8.17 Discrete Fourier transform of an exponential signal e−2*^t u*(t).

814 CHAPTER 8 SAMPLING: THE BRIDGE FROM CONTINUOUS TO DISCRETE

In this example, we knew X(ω) beforehand; hence we could make intelligent choices for B (or the sampling frequency fs). In practice, we generally do not know X(ω) beforehand. In fact, that is the very thing we are trying to determine. In such a case, we must make an intelligent guess for B or fs from circumstantial evidence. We should then continue reducing the value of T and recomputing the transform until the result stabilizes within the desired number of significant digits.

USING MATLAB TO COMPUTE AND PLOT THE RESULTS

Let us now use MATLAB to confirm the results of this example. First, parameters are defined and MATLAB’s fft command is used to compute the DFT.

>> T_0 = 4; N_0 = 256; T = T_0/N_0; t = (0:T:T*(N_0-1))';
>> x = T*exp(-2*t); x(1) = x(1)/2;
>> X_r = fft(x); r = [-N_0/2:N_0/2-1]'; omega_r = r*2*pi/T_0;

The true Fourier transform is also computed for comparison.

omega = linspace(-pi/T,pi/T,5001); X = 1./(j*omega+2);

For clarity, we display spectrum over a restricted frequency range.

>> subplot(1,2,1); stem(omega_r,fftshift(abs(X_r)),'k.');

>> line(omega,abs(X),'color',[0 0 0]); axis([-0.01 44 -0.01 0.51]);

• xlabel(‘\omega’); ylabel(‘|X(\omega)|’);

• subplot(1,2,2); stem(omega_r,fftshift(angle(X_r)),‘k.’);

• line(omega,angle(X),‘color’,[0 0 0]); axis([-0.01 44 -pi/2-0.01 0.01]);

• xlabel(‘\omega’); ylabel(‘\angle X(\omega)’);

The results, shown in Fig. 8.18, match the earlier results shown in Fig. 8.17.

EXAMPLE 8.9 DFT to Compute the Fourier Transform of a Rectangular Pulse

Use the DFT to compute the Fourier transform of 8 rect(t).

This gate function and its Fourier transform are illustrated in Figs. 8.19a and 8.19b. To determine the value of the sampling interval T, we must first decide on the essential bandwidth B. In Fig. 8.19b, we see that X(ω) decays rather slowly with ω. Hence, the essential bandwidth B is rather large. For instance, at B = 15.5 Hz (97.39 rad/s), X(ω) = −0.1643, which is about 2% of the peak at X(0). Hence, the essential bandwidth is well above 16 Hz if we use the 1% of the peak amplitude criterion for computing the essential bandwidth. However, we shall deliberately take B = 4 for two reasons: to show the effect of aliasing and because the use of B > 4 would give an enormous number of samples, which could not be conveniently displayed on the page without losing sight of the essentials. Thus, we shall intentionally accept approximation to graphically clarify the concepts of the DFT.

The choice of B = 4 results in the sampling interval T = 1/2B = 1/8. Looking again at the spectrum in Fig. 8.19b, we see that the choice of the frequency resolution f⁰ = 1/4 Hz is reasonable. Such a choice gives us four samples in each lobe of X(ω). In this case T⁰ = 1/f⁰ = 4 seconds and N⁰ = T0/T = 32. The duration of x(t) is only 1 second. We must repeat it every 4 seconds (T⁰ = 4), as depicted in Fig. 8.19c, and take samples every 1/8 second. This choice yields 32 samples (N⁰ = 32). Also,

$$x_n = Tx(nT) = \frac{1}{8}x(nT)$$

Since x(t) = 8 rect(t), the values of xn are 1, 0, or 0.5 (at the points of discontinuity), as illustrated in Fig. 8.19c, where xn is depicted as a function of t as well as n, for convenience.

In the derivation of the DFT, we assumed that x(t) begins at t = 0 (Fig. 8.16a), and then took N⁰ samples over the interval (0, T0). In the present case, however, x(t) begins at −1/2. This difficulty is easily resolved when we realize that the DFT obtained by this procedure is actually the DFT of xn repeating periodically every T⁰ seconds. Figure 8.19c clearly indicates that periodic repeating the segment of xn over the interval from −2 to 2 seconds yields the same signal as the periodic repeating the segment of xn over the interval from 0 to 4 seconds. Hence, the DFT of the samples taken from −2 to 2 seconds is the same as that of the samples taken from 0 to 4 seconds. Therefore, regardless of where x(t) starts, we can always take the samples of x(t) and its periodic extension over the interval from 0 to T0. In the present example, the 32 sample values are

$$x_n = \begin{cases} 1 & 0 \le n \le 3 \text{ and } 29 \le n \le 31 \\ 0 & 5 \le n \le 27 \\ 0.5 & n = 4,28 \end{cases}$$

816 CHAPTER 8 SAMPLING: THE BRIDGE FROM CONTINUOUS TO DISCRETE

Figure 8.19 Discrete Fourier transform of a gate pulse.

Observe that the last sample is at t = 31/8, not at 4, because the signal repetition starts at t = 4, and the sample at t = 4 is the same as the sample at t = 0. Now, N⁰ = 32 and ⁰ = 2π/32 = π/16. Therefore [see Eq. (8.12)],

$$X_r = \sum_{n=0}^{31} x_n e^{-jr(\pi/16)n}$$

Values of Xr are computed according to this equation and plotted in Fig. 8.19d.

The samples Xr are separated by f⁰ = 1/T⁰ Hz. In this case T⁰ = 4, so the frequency resolution f⁰ is 1/4 Hz, as desired. The folding frequency fs/2 = B = 4 Hz corresponds to r = N0/2 = 16. Because Xr is N⁰ periodic (N⁰ = 32), the values of Xr for r = −16 to n = −1 are the same as those for r = 16 to n = 31. For instance, X¹⁷ = X−15, X¹⁸ = X−14, and so on. The DFT gives us the samples of the spectrum X(ω).

For the sake of comparison, Fig. 8.19d also shows the shaded curve 8 sinc(ω/2), which is the Fourier transform of 8 rect(t). The values of Xr computed from the DFT equation show aliasing error, which is clearly seen by comparing the two superimposed plots. The error in X² is just about 1.3%. However, the aliasing error increases rapidly with r. For instance, the error in X⁶ is about 12%, and the error in X¹⁰ is 33%. The error in X¹⁴ is a whopping 72%. The percent error increases rapidly near the folding frequency (r = 16) because x(t) has a jump discontinuity, which makes X(ω) decay slowly as 1/ω. Hence, near the folding frequency, the inverted tail (due to aliasing) is very nearly equal to X(ω) itself. Moreover, the final values are the difference between the exact and the folded values (which are very close to the exact values). Hence, the percent error near the folding frequency (r = 16 in this case) is very high, although the absolute error is very small. Clearly, for signals with jump discontinuities, the aliasing error near the folding frequency will always be high (in percentage terms), regardless of the choice of N0. To ensure a negligible aliasing error at any value r, we must make sure that N⁰ r. This observation is valid for all signals with jump discontinuities.

USING MATLAB TO COMPUTE AND PLOT THE RESULTS

Once again, MATLAB lets us easily confirm the results of this example. First, parameters are defined and MATLAB’s fft command is used to compute the DFT.

T_0 = 4; N_0 = 32; T = T_0/N_0; >> x_n = [ones(1,4) 0.5 zeros(1,23) 0.5 ones(1,3)]’; >> X_r = fft(x_n); r = [-N_0/2:N_0/2-1]’; omega_r = r*2*pi/T_0;

The true Fourier transform is also computed for comparison.

>> omega = linspace(-pi/T,pi/T,5001); X = 8*sinc(omega/(2*pi));

Since it is real, we can display the resulting spectrum using a single plot.

clf; stem(omega_r,fftshift(real(X_r)),‘k.’); >> line(omega,X,‘color’,[0 0 0]);

818 CHAPTER 8 SAMPLING: THE BRIDGE FROM CONTINUOUS TO DISCRETE

xlabel(‘\omega’); ylabel(‘X(\omega)’); axis tight

The result, shown in Fig. 8.20, matches the earlier result shown in Fig. 8.19d. The DFT approximation does not perfectly follow the true Fourier transform, especially at high frequencies, because the parameter B is deliberately set too small.

8.5-1 Some Properties of the DFT

The discrete Fourier transform is basically the Fourier transform of a sampled signal repeated periodically. Hence, the properties derived earlier for the Fourier transform apply to the DFT as well.

LINEARITY If xn ⇐⇒ Xr and gn ⇐⇒ Gr, then

$$a_1x_n + a_2g_n \Longleftrightarrow a_1X_r + a_2G_r$$

The proof is trivial.

CONJUGATE SYMMETRY

From the conjugation property x∗(t) ⇐⇒ X∗(−ω), we have

x∗ ⁿ ←→ X^∗ −r

From this equation and the time-reversal property, we obtain

$$x_{-n}^* \longleftrightarrow X_r^*$$

When x(t) is real, then the conjugate-symmetry property states that X∗(ω) = X(−ω). Hence, for real xn,

$$X_r^* = X_{-r}$$

Moreover, Xr is N⁰ periodic. Hence,

$$X_r^* = X_{N_0-r}$$

Because of this property, we need compute only half the DFTs for real xn. The other half are the conjugates.

TIME SHIFTING

The time-shifting (circular shifting) property states†

$$x_{n-k} \Longleftrightarrow X_r e^{-jr\Omega_0 k}$$

Proof. We use Eq. (8.13) to find the inverse DFT of Xre−jr0*^k* as

$$\frac{1}{N_0} \sum_{r=0}^{N_0-1} X_r e^{-jr\Omega_0 k} e^{jr\omega_0 n} = \frac{1}{N_0} \sum_{r=0}^{N_0-1} X_r e^{jr\Omega_0 (n-k)} = x_{n-k}$$

FREQUENCY SHIFTING

A dual of the time-shifting property, the frequency-shifting property states

$$x_n e^{jn\Omega_0 m} \Longleftrightarrow X_{r-m}$$

Proof. This proof is identical to that of the time-shifting property except that we start with Eq. (8.12).

CIRCULAR CONVOLUTION

The circular (or periodic) convolution property states

$$x_n \circledast g_n \Longleftrightarrow X_r G_r \tag{8.18}$$

and

$$x_n g_n \Longleftrightarrow \frac{1}{N_0} X_r \circledast G_r \tag{8.19}$$

For two N0-periodic sequences xn and gn, circular (or periodic) convolution is defined by

$$x_n \circledast g_n = \sum_{k=0}^{N_0 - 1} x_k g_{n-k} = \sum_{k=0}^{N_0 - 1} g_k x_{n-k}$$

\n(8.20)

^† Time shifting is also known as circular shifting because such a shift can be interpreted as a circular shift of the N⁰ samples in the first cycle 0 ≤ n ≤ N⁰ −1.

Figure 8.21 Graphical depictions of circular convolution.

To prove Eq. (8.18), we find the DFT of the circular convolution xn-∗ gn as

$$\sum_{n=0}^{N_0-1} \left( \sum_{k=0}^{N_0-1} x_k g_{n-k} \right) e^{-j r \omega_0 n} = \sum_{k=0}^{N_0-1} x_k \left( \sum_{n=0}^{N_0-1} g_{n-k} e^{-j r \omega_0 n} \right)$$ $$= \sum_{k=0}^{N_0-1} x_k (G_r e^{-j r \Omega_0 k}) = X_r G_r$$

Equation (8.19) can be proved in the same way.

For periodic sequences, the convolution can be visualized in terms of two sequences, with one sequence fixed and the other inverted and moved past the fixed sequence, one digit at a time. If the two sequences are N⁰ periodic, the same configuration will repeat after N⁰ shifts of the sequence. Clearly the convolution xn-∗ gn becomes N⁰ periodic. Such convolution can be conveniently visualized in terms of N⁰ sequences, as illustrated in Fig. 8.21, for the case of N⁰ = 4. The inner N0-point sequence xn is clockwise and fixed. The outer N0-point sequence gn is inverted so that it becomes counterclockwise. This sequence is now rotated clockwise 1 unit at a time. We multiply the overlapping numbers and add. For example, the value of xn-∗ gn at n = 0 (Fig. 8.21) is

$$x_0g_0 + x_1g_3 + x_2g_2 + x_3g_1$$

and the value of xn-∗ gn at n = 1 is (Fig. 8.21)

$$x_0g_1 + x_1g_0 + x_2g_3 + x_3g_2$$

and so on.

8.5-2 Some Applications of the DFT

The DFT is useful not only in the computation of direct and inverse Fourier transforms, but also in other applications such as convolution, correlation, and filtering. Use of the efficient FFT algorithm, discussed shortly (Sec. 8.6), makes it particularly appealing.

LINEAR CONVOLUTION

Let x(t) and g(t) be the two signals to be convolved. In general, these signals may have different time durations. To convolve them by using their samples, they must be sampled at the same rate (not below the Nyquist rate of either signal). Let xn (0 ≤ n ≤ N¹ − 1) and gn (0 ≤ n ≤ N² − 1) be the corresponding discrete sequences representing these samples. Now,

$$c(t) = x(t) * g(t)$$

and if we define three sequences as xn = Tx(nT), gn = Tg(nT), and cn = Tc(nT), then†

$$c_n = x_n * g_n$$

where we define the linear convolution sum of two discrete sequences xn and gn as

$$c_n = x_n * g_n = \sum_{k=-\infty}^{\infty} x_k g_{n-k}$$

Because of the width property of the convolution, cn exists for 0≤ n≤ N1+N2−1. To be able to use the DFT circular convolution technique, we must make sure that the circular convolution will yield the same result as does linear convolution. In other words, the signal resulting from the circular convolution must have the same length (N¹ + N² − 1) as that of the signal resulting from linear convolution. This step can be accomplished by adding N² − 1 dummy samples of zero value to xn and N¹ −1 dummy samples of zero value to gn (zero padding). This procedure changes the length of both xn and gn to N1+N² −1. The circular convolution now is identical to the linear convolution except that it repeats periodically with period N¹ +N² −1. A little reflection will show that in such a case the circular convolution procedure in Fig. 8.21 over one cycle (0 ≤ n ≤ N¹ + N² − 1) is identical to the linear convolution of the two sequences xn and gn. We can use the DFT to find the convolution xn ∗ gn in three steps, as follows:

• 1. Find the DFTs Xr and Gr corresponding to suitably padded xn and gn.
• 2. Multiply Xr by Gr.
• 3. Find the IDFT of XrGr. This procedure of convolution, when implemented by the fast Fourier transform algorithm (discussed later), is known as fast convolution.

FILTERING

We generally think of filtering in terms of a hardware-oriented solution (e.g., building a circuit with RLC components and operational amplifiers). However, filtering also has a software-oriented solution [a computer algorithm that yields the filtered output y(t) for a given input x(t)]. This goal can be conveniently accomplished by using the DFT. If x(t) is the signal to be filtered, then Xr, the DFT of xn, is found. The spectrum Xr is then shaped (filtered) as desired by multiplying Xr by Hr, where Hr are the samples of H(ω) for the filter [Hr = H(rω0)]. Finally, we take the IDFT of XrHr to obtain the filtered output yn[yn =Ty(nT)]. This procedure is demonstrated in the following example.

^† We can show that cn ⁼ limT→⁰ xn ^∗ gn; [see 4]. Error is inherent in any numerical method used to compute convolution of continuous-time signals; since T = 0 in practice, there will be some error in this equation.

The signal x(t) in Fig. 8.22a is passed through an ideal lowpass filter of frequency response H(ω) depicted in Fig. 8.22b. Use the DFT to find the sampled version of the filter output.

Figure 8.22 DFT solution for filtering x(t) through H(ω).

We have already found the 32-point DFT of x(t) (see Fig. 8.19d). Next we multiply Xr by Hr. To find Hr, we recall using f⁰ = 1/4 in computing the 32-point DFT of x(t). Because Xr is 32-periodic, Hr must also be 32-periodic with samples separated by 1/4 Hz. This fact means that Hr must be repeated every 8 Hz or 16π rad/s (see Fig. 8.22c). The resulting 32 samples of Hr over (0 ≤ ω ≤ 16π ) are as follows:

$$H_r = \begin{cases} 1 & 0 \le r \le 7 \\ 0 & 9 \le r \le 23 \\ 0.5 & r = 8,24 \end{cases} \text{ and } 25 \le r \le 31$$

We multiply Xr with Hr. The desired output signal samples yn are found by taking the inverse DFT of XrHr. The resulting output signal is illustrated in Fig. 8.22d.

It is quite simple to verify the results of this filtering example using MATLAB. First, parameters are defined, and MATLAB’s fft command is used to compute the DFT of xn.

T_0 = 4; N_0 = 32; T = T_0/N_0; n = (0:N_0-1); r = n; >> x_n = [ones(1,4) 0.5 zeros(1,23) 0.5 ones(1,3)]’; X_r = fft(x_n);

The DFT of the filter’s output is just the product of the filter response Hr and the input DFT Xr. The output yn is obtained using the ifft command and then plotted.

>> H_r = [ones(1,8) 0.5 zeros(1,15) 0.5 ones(1,7)]';
>> Y_r = H_r.*X_r; y_n = ifft(Y_r);
>> clf; stem(n,real(y_n),'k.');
>> xlabel('n'); ylabel('y_n'); axis([0 31 -.1 1.1]);

The result, shown in Fig. 8.23, matches the earlier result shown in Fig. 8.22d. Recall, this DFT-based approach shows the samples yn of the filter output y(t) (sampled in this case at a rate ^{T = 1 8} ) over 0 ≤ n ≤ N⁰ − 1 = 31 when the input pulse x(t) is periodically replicated to form samples xn (see Fig. 8.19c).

Figure 8.23 Using MATLAB and the DFT to determine filter output.