[9.7-1 Computing the Discrete-Time Fourier Series](#page-14-0)

← Back to LINEAR SYSTEMS AND SIGNALS Overview

9.7-1 Computing the Discrete-Time Fourier Series

Within a scale factor, the DTFS is identical to the DFT. Thus, methods to compute the DFT can be readily used to compute the DTFS. Specifically, the DTFS is the DFT scaled by 1/N0. As an example, consider a 50 Hz sinusoid sampled at 1000 Hz over one-tenth of a second.

>> T = 1/1000; N_0 = 100; n = (0:N_0-1)';
>> x = cos(2*pi*50*n*T);

The DTFS is obtained by scaling the DFT.

>> X = fft(x)/N_0; f = (0:N_0-1)/(T*N_0);
>> stem(f-1/(2*T),fftshift(abs(X)),'k.');
>> axis([-500 500 -0.05 0.55]); xlabel('f [Hz]'); ylabel('|X(f)|');

Figure 9.19 shows a peak magnitude of 0.5 at ±50 Hz. This result is consistent with Euler’s representation

$$\cos(2\pi 50nT) = \frac{1}{2}e^{j2\pi 50nT} + \frac{1}{2}e^{-j2\pi 50nT}$$

Lacking the 1/N⁰ scale factor, the DFT would have a peak amplitude 100 times larger.

The inverse DTFS is obtained by scaling the inverse DFT by N0.

>> x = real(ifft(X)*N_0); stem(n,x,'k.');
>> axis([0 99 -1.1 1.1]); xlabel('n'); ylabel('x[n]');

Figure 9.20 confirms that the sinusoid x[n] is properly recovered. Although the result is theoretically real, computer round-off errors produce a small imaginary component, which the real command removes.

Figure 9.19 DTFS computed by scaling the DFT.

Figure 9.20 Inverse DTFS computed by scaling the inverse DFT.

Although MATLAB’s fft command provides an efficient method to compute the DTFS, other important computational methods exist. A matrix-based approach is one popular way to implement Eq. (9.4). Although not as efficient as an FFT-based algorithm, matrix-based approaches provide insight into the DTFS and serve as an excellent model for solving similarly structured problems.

To begin, define WN⁰ = ej⁰ , which is a constant for a given N0. Substituting WN⁰ into Eq. (9.4) yields

$$\mathcal{D}_r = \frac{1}{N_0} \sum_{n=0}^{N_0 - 1} x[n] W_{N_0}^{-nr}$$

An inner product of two vectors computes Dr.

$$\mathcal{D}_r = \frac{1}{N_0} \begin{bmatrix} 1 & W_{N_0}^{-r} & W_{N_0}^{-2r} & \dots & W_{N_0}^{-(N_0-1)r} \end{bmatrix} \begin{bmatrix} x[0] \\ x[1] \\ x[2] \\ \vdots \\ x[N_0-1] \end{bmatrix}$$

Stacking the results for all r yields

\begin{bmatrix}\n\begin{bmatrix}\n\mathcal{D}_0 \\ \mathcal{D}_1 \\ \mathcal{D}_2 \\ \vdots \\ \mathcal{D}_{N_0-1}\n\end{bmatrix} = \frac{1}{N_0} \begin{bmatrix}\n1 & 1 & 1 & \cdots & 1 \\ 1 & W_{N_0}^{-1} & W_{N_0}^{-2} & \cdots & W_{N_0}^{-(N_0-1)} \\ 1 & W_{N_0}^{-2} & W_{N_0}^{-4} & \cdots & W_{N_0}^{-2(N_0-1)} \\ \vdots & \vdots & \vdots & \cdots & \vdots \\ 1 & W_{N_0}^{-(N_0-1)} & W_{N_0}^{-2(N_0-1)} & \cdots & W_{N_0}^{-(N_0-1)^2}\n\end{bmatrix} \begin{bmatrix}\nx[0] \\ x[1] \\ x[2] \\ \vdots \\ x[N_0-1]\n\end{bmatrix}

In matrix notation, this equation is compactly written as

$$\mathbf{D} = \frac{1}{N_0} \mathbf{W}_{N_0} \mathbf{x}$$

Since it is also used to compute the DFT, matrix W*^N*⁰ is often called a DFT matrix.

Let us create an anonymous function to compute the N0-by-N⁰ DFT matrix W*^N*⁰ . Although not used here, the signal-processing toolbox function dftmtx computes the same DFT matrix, although in a less obvious but more efficient fashion.

W = @(N_0) (exp(-j*2*pi/N_0)).^((0:N_0-1)’*(0:N_0-1));

While less efficient than FFT-based methods, the matrix approach correctly computes the DTFS.

>> X = W(N_0)*x/N_0; stem(f-1/(2*T),fftshift(abs(X)),'k.');
>> axis([-500 500 -0.05 0.55]); xlabel('f [Hz]'); ylabel('|X(f)|');

The resulting plot is indistinguishable from Fig. 9.19. Problem 9.7-1 investigates a matrix-based approach to compute Eq. (9.3), the inverse DTFS.

9.7-2 Measuring Code Performance

Writing efficient code is important, particularly if the code is frequently used, requires complicated operations, involves large data sets, or operates in real time. MATLAB provides several tools for assessing code performance. When properly used, the profile function provides detailed statistics that help assess code performance. MATLAB help thoroughly describes the use of the sophisticated profile command.

A simpler method of assessing code efficiency is to measure execution time and compare it with a reference. The MATLAB command tic starts a stopwatch timer. The toc command reads the timer. Sandwiching instructions between tic and toc returns the elapsed time. For example, the execution time of the 100-point matrix-based DTFS computation is

>> tic; W(N_0)*x/N_0; toc
    Elapsed time is 0.004417 seconds.

Different machines operate at different speeds with different operating systems and with different background tasks. Therefore, elapsed-time measurements can vary considerably from machine to machine and from execution to execution. For relatively simple and short events like the present case, execution times can be so brief that MATLAB may report unreliable times or fail to register an elapsed time at all.

To increase the elapsed time and therefore the accuracy of the time measurement, a loop is used to repeat the calculation.

>> tic; for i=1:100, W(N_0)*x/N_0; end; toc
    Elapsed time is 0.173388 seconds.

This elapsed time suggests that each 100-point DTFS calculation takes a little under 2 milliseconds. What exactly does this mean, however? Elapsed time is only meaningful relative to some reference. Let us see what difference occurs by precomputing the DFT matrix, rather than repeatedly using our anonymous function.

>> W100 = W(100); tic; for i=1:100, W100*x/N_0; end; toc
   Elapsed time is 0.001199 seconds.

892 CHAPTER 9 FOURIER ANALYSIS OF DISCRETE-TIME SIGNALS

Amazingly, this small change makes a hundredfold change in our computational efficiency! Clearly, it is much better to precompute the DFT matrix.

To provide another example, consider the time it takes to compute the same DTFS using the FFT-based approach.

tic; for i=1:100, fft(x)/N_0; end; toc Elapsed time is 0.000399 seconds.

With this as a reference, our fastest matrix-based computations appear to be several times slower than the FFT-based computations. This difference becomes more dramatic as N⁰ is increased. Since the two methods provide identical results, there is little incentive to use the slower matrix-based approach, and the FFT-based algorithm is generally preferred. Even so, the FFT can exhibit curious behavior: adding a few data points, even the artificial samples introduced by zero padding, can dramatically increase or decrease execution times. The tic and toc commands illustrate this strange result. Consider computing the DTFS of 1015 random data points 100 times.

x1 = rand(1015,1); tic; for i=1:100; fft(x1)/1015; end; T1 = toc T1 = 0.0067

Next, pad the sequence with four zeros.

$$x2 = [x1; zeros(4,1)]$$

; tic; for i=1:100; fft( $x2$ )/1019; end; T2 = toc
T2 = 0.0134

The ratio of the two elapsed times indicates that adding four points to an already long sequence increases the computation time by a factor of 2. Next, the sequence is zero-padded to a length of N⁰ = 1024.

x3 = [x2;zeros(5,1)]; tic; for i=1:100; fft(x3)/1024; end; T3 = toc T3 = 0.0017

In this case, the added data decrease the original execution time by a factor of 4 and the second execution time by a factor of 8! These results are particularly surprising when it is realized that the lengths of y1, y2, and y3 differ by less than 1%.

As it turns out, the efficiency of the fft command depends on the factorability of N0. With the factor command, 1015 = (5)(7)(29), 1019 is prime, and 1024 = (2)10. The most factorable length, 1024, results in the fastest execution, while the least factorable length, 1019, results in the slowest execution. To ensure the greatest factorability and fastest operation, vector lengths are ideally a power of 2.

9.7-3 FIR Filter Design by Frequency Sampling

Finite impulse response (FIR) digital filters are flexible, always stable, and relatively easy to implement. These qualities make FIR filters a popular choice among digital filter designers. The difference equation of a length-N causal FIR filter is conveniently expressed as

$$y[n] = h_0 x[n] + h_1 x[n-1] + \dots + h_{N-1} x[n-(N-1)] = \sum_{k=0}^{N-1} h_k x[n-k]$$

The filter coefficients, or tap weights as they are sometimes called, are expressed by using the variable h to emphasize that the coefficients themselves represent the impulse response of the filter.

The filter’s frequency response is

$$H(\Omega) = \frac{Y(\Omega)}{X(\Omega)} = \sum_{k=0}^{N-1} h_k e^{-j\Omega k}$$

Since H() is a 2π-periodic function of the continuous variable , it is sufficient to specify H() over a single period (0 ≤ < 2π ).

In many filtering applications, the desired magnitude response |Hd()| is known but not the filter coefficients h[n]. The question, then, is one of determining the filter coefficients from the desired magnitude response.

Consider the design of a lowpass filter with cutoff frequency ^c = π/4. An anonymous function represents the desired ideal frequency response.

H_d = @(Omega) (mod(Omega,2*pi)<pi/4)+(mod(Omega,2*pi)>2*pi-pi/4);

Since the inverse DTFT of Hd() is a sampled sinc function, it is impossible to perfectly achieve the desired response with a causal, finite-length FIR filter. A realizable FIR filter is necessarily an approximation, and an infinite number of possible solutions exist. Thought of another way, Hd() specifies an infinite number of points, but the FIR filter only has N unknown tap weights. In general, we expect a length-N filter to match only N points of the desired response over (0 ≤ < 2π ). Which frequencies should be chosen?

A simple and sensible method is to select N frequencies uniformly spaced on the interval (0 ≤ < 2π ), (0, 2π/N, 4π/N, 6π/N,…,(N −1)2π/N). By choosing uniformly spaced frequency samples, the N-point inverse DFT can be used to determine the tap weights h[n]. Program CH9MP1 illustrates this procedure.

function [h] = CH9MP1(N,H_d);
% CH9MP1.m : Chapter 9, MATLAB Program 1
% Function M-file designs a length-N FIR filter by sampling the desired
% magnitude response H_d. Phase response is left as zero.
% INPUTS: N = desired FIR filter length
% H_d = anonymous function that defines the desired magnitude response
% OUTPUTS: h = impulse response (FIR filter coefficients)
% Create N equally spaced frequency samples:
Omega = linspace(0,2*pi*(1-1/N),N)';
% Sample the desired magnitude response and create h[n]:
H = 1.0*H_d(Omega); h = real(ifft(H));

To complete the design, the filter length must be specified. Small values of N reduce the filter’s complexity but also reduce the quality of the filter’s response. Large values of N improve the approximation of Hd() but also increase complexity. A balance is needed. We choose an intermediate value of N = 21 and use CH9MP1 to design the filter.

N = 21; h = CH9MP1(N,H_d);

To assess the filter quality, the frequency response is computed by means of program CH5MP1.

>> Omega = linspace(0,2*pi,1000); samples = linspace(0,2*pi*(1-1/N),N)';

>> H = CH5MP1(h,1,Omega);

>> subplot(2,1,1); stem([0:N-1],h,'k.'); xlabel('n'); ylabel('h[n]');

>> subplot(2,1,2);

plot(samples,H_d(samples),‘k.’,Omega,H_d(Omega),‘k:‘,Omega,abs(H),‘k’);

>> axis([0 2*pi -0.1 1.6]); xlabel('\Omega'); ylabel('|H(\Omega)|');

legend(‘Samples’,‘Desired’,‘Actual’,‘Location’,‘North’);

As shown in Fig. 9.21, the filter’s frequency response intersects the desired response at the sampled values of Hd(). The overall response, however, has significant ripple between sample points that renders the filter practically useless. Increasing the filter length does not alleviate the ripple problems. Figure 9.22 shows the case N = 41.

To understand the poor behavior of filters designed with CH9MP1, remember that the impulse response of an ideal lowpass filter is a sinc function with the peak centered at zero. Thought of another way, the peak of the sinc is centered at n = 0 because the phase of Hd() is zero. Constrained to be causal, the impulse response of the designed filter still has a peak at n = 0 but cannot include values for negative n. As a result, the sinc function is split in an unnatural way with sharp discontinuities on both ends of h[n]. Sharp discontinuities in the time domain appear as high-frequency oscillations in the frequency domain, which is why H() has significant ripple.

To improve the filter behavior, the peak of the sinc is moved to n = (N − 1)/2, the center of the length-N filter response. In this way, the peak is not split, no large discontinuities are present, and frequency response ripple is consequently reduced. From DFT properties, a cyclic shift of (N − 1)/2 in the time domain requires a scale factor of e−j(N−1)/² in the frequency

Figure 9.21 Length-21 FIR lowpass filter using zero phase.

Figure 9.22 Length-41 FIR lowpass filter using zero phase.

domain.† Notice that the scale factor e−j(N−1)/² affects only phase, not magnitude, and results in a linear phase filter. Program CH9MP2 implements the procedure.

function [h] = CH9MP2(N,H_d);
% CH9MP2.m : Chapter 9, MATLAB Program 2
% Function M-file designs a length-N FIR filter by sampling the desired
% magnitude response H_d. Phase is defined to shift h[n] by (N-1)/2.
% INPUTS: N = desired FIR filter length
% H_d = anonymous function that defines the desired magnitude response
% OUTPUTS: h = impulse response (FIR filter coefficients)
% Create N equally spaced frequency samples and use to sample H_d:
Omega = linspace(0,2*pi*(1-1/N),N)'; H = H_d(Omega);
% Define phase to shift h[n] by (N-1)/2:
H = H.*exp(-j*Omega*((N-1)/2));
H(fix(N/2)+2:N,1) = H(fix(N/2)+2:N,1)*((-1)^(N-1));
h = real(ifft(H));

^† Technically, the shift property requires (^{N −}1)/2 to be an integer, which occurs only for odd-length filters. The next-to-last line of program CH9MP2 implements a correction factor, of sorts, required to accommodate the fractional shifts desired for even-length filters. The mathematical derivation of this correction is nontrivial and is not included here. Those hesitant to use this correction factor have an alternative: simply round (N − 1)/2 to the nearest integer. Although the rounded shift is slightly off-center for even-length filters, there is usually little or no appreciable difference in the characteristics of the filter. Even so, true centering is desirable because the resulting impulse response is symmetric, which can reduce by half the number of multiplies required to implement the filter.

Figure 9.23 Length-21 FIR lowpass filter using linear phase.

Figure 9.23 shows the results for the N = 21 case using CH9MP2 to compute h[n]. As hoped, the impulse response looks like a sinc function with the peak centered at n = 10. Additionally, the frequency response ripple is greatly reduced. With CH9MP2, increasing N improves the quality of the filter, as shown in Fig. 9.24 for the case N = 41. While the magnitude response is needed to establish the general shape of the filter response, it is the proper selection of phase that ensures the acceptability of the filter’s behavior.

To illustrate the flexibility of the design method, consider a bandpass filter with passband (π/4 < || < π/2).

>> H_d = @(Omega) (mod(Omega,2*pi)>pi/4)&(mod(Omega,2*pi)<pi/2)+...
>> (mod(Omega,2*pi)>3*pi/2)&(mod(Omega,2*pi)<7*pi/4);

Figure 9.25 shows the results for N = 50. Notice that this even-length filter uses a fractional shift and is symmetric about n = 24.5.

Although FIR filter design by means of frequency sampling is very flexible, it is not always appropriate. Extreme care is needed for filters, such as digital differentiators and Hilbert transformers, that require special phase characteristics for proper operation. Additionally, if frequency samples occur near jump discontinuities of Hd(), rounding errors may, in rare cases, disrupt the desired symmetry of the sampled magnitude response. Such cases are corrected by slightly adjusting the location of problematic jump discontinuities or by changing the value of N.

Figure 9.24 Length-41 FIR lowpass filter using linear phase.

Figure 9.25 Length-50 FIR bandpass filter using linear phase.