6.6 NUMERICAL [COMPUTATION OF](#page-12-0) Dn

← Back to LINEAR SYSTEMS AND SIGNALS Overview

6.6 NUMERICAL COMPUTATION OF Dn

We can compute Dn numerically by using the DFT (the discrete Fourier transform discussed in Sec. 8.5), which uses the samples of a periodic signal x(t) over one period. The sampling interval is T seconds. Hence, there are N⁰ =T0/T number of samples in one period T0. To find the relationship between Dn and the samples of x(t), consider Eq. (6.19) and write

$$D_n = \frac{1}{T_0} \int_{T_0} x(t)e^{-jn\omega_0 t} dt$$

=

$$\lim_{T \to 0} \frac{1}{N_0 T} \sum_{k=0}^{N_0 - 1} x(kT)e^{-jn\omega_0 kT} T$$

=

$$\lim_{T \to 0} \frac{1}{N_0} \sum_{k=0}^{N_0 - 1} x(kT)e^{-jn\Omega_0 k}$$

(6.55)

where x(kT) is the kth sample of x(t) and

$$N_0 = \frac{T_0}{T} \quad \text{and} \quad \Omega_0 = \omega_0 T = \frac{2\pi}{N_0}$$

In practice, it is impossible to make T → 0 in computing the right-hand side of Eq. (6.55). We can make T small, but not zero, which will cause the data to increase without limit. Thus, we shall ignore the limit on T in Eq. (6.55) with the implicit understanding that T is reasonably small. Nonzero T will result in some computational error, which is inevitable in any numerical evaluation of an integral. The error resulting from nonzero T is called the aliasing error, which is discussed in more detail in Ch. 8. Thus, we can express Eq. (6.55) as

$$D_n \approx \frac{1}{N_0} \sum_{k=0}^{N_0 - 1} x(kT) e^{-jn\Omega_0 k}$$

(6.56)

Since 0N⁰ = 2π, we know that ejn0(k+N0) = ejn0*^k*, and it follows that

$$D_{n+N_0}=D_n$$

The periodicity property Dn+N⁰ = Dn means that beyond n = N0/2, the coefficients represent the values for negative n. For instance, when N⁰ = 32, D¹⁷ = D−15, D¹⁸ = D−14,…,D³¹ = D−1. The cycle repeats again from n = 32 on.

We can use the efficient FFT (the fast Fourier transform discussed in Sec. 8.6) to compute the right-hand side of Eq. (6.56). We shall use MATLAB to implement the FFT algorithm. For

660 CHAPTER 6 CONTINUOUS-TIME SIGNAL ANALYSIS: THE FOURIER SERIES

this purpose, we need samples of x(t) over one period starting at t = 0. In this algorithm, it is also preferable (although not necessary) that N⁰ be a power of 2, (i.e., N⁰ = 2*^m* where m is an integer).

EXAMPLE 6.16 Numerical Computation of Fourier Spectra

Numerically compute and then plot the exponential Fourier spectra for the periodic signal in Fig. 6.2a (Ex. 6.1).

The samples of x(t) start at t = 0 and the last (N0th) sample is at t = T⁰ − T. At the points of discontinuity, the sample value is taken as the average of the values of the function on two sides of the discontinuity. Thus, the sample at t = 0 is not 1 but (e−π/² +1)/2 = 0.604. To determine N0, we require that Dn for n ≥ N0/2 be negligible. Because x(t) has a jump discontinuity, Dn decays rather slowly as 1/n. Hence, a choice of N⁰ = 200 is acceptable because the (N0/2)nd (100th) harmonic is about 1% of the fundamental. However, we also require N⁰ to be a power of 2. Hence, we shall take N⁰ = 256 = 28.

First, the basic parameters are established.

T_0 = pi; N_0 = 256; T = T_0/N_0; t = (0:T:T*(N_0-1))’; >> x = exp(-t/2); x(1) = (exp(-pi/2)+1)/2;

Next, the DFT, computed by means of the fft function, is used to approximate the exponential Fourier spectra up to n = N0/2. To facilitate comparison with previous plots of Dn, we only plot the results over −5 ≤ n ≤ 5.

>> D_n = fft(x)/N_0; n = [-N_0/2:N_0/2-1]';
>> clf; subplot(1,2,1); stem(n,abs(fftshift(D_n)),'.k');
>> axis([-5 5 0 .6]); xlabel('n'); ylabel('|D_n|');
>> subplot(1,2,2); stem(n,angle(fftshift(D_n)),'.k');
>> axis([-5 5 -2 2]); xlabel('n'); ylabel('\angle D_n [rad]');

As shown in Fig. 6.28, the resulting approximation is visually indistinguishable from the true Fourier series spectra shown in Fig. 6.12 or Fig. 6.13.