[8.7 MATLAB: THE](#page-13-0) DISCRETE FOURIER TRANSFORM

← Back to LINEAR SYSTEMS AND SIGNALS Overview

8.7 MATLAB: THE DISCRETE FOURIER TRANSFORM

As an idea, the discrete Fourier transform (DFT) has been known for hundreds of years. Practical computing devices, however, are responsible for bringing the DFT into common use. MATLAB is capable of DFT computations that would have been impractical just a few decades ago.

8.7-1 Computing the Discrete Fourier Transform

The MATLAB command fft(x) computes the DFT of a vector x that is defined over (0 ≤ n ≤ N⁰ −1) (Problem 8.7-1 considers how to scale the DFT to accommodate signals that do not begin at n = 0.) As its name suggests, the function fft uses the computationally more efficient fast Fourier transform algorithm when it is appropriate to do so. The inverse DFT is easily computed by using the ifft function.

^† Actually, N0/2 is a conservative figure because some multiplications corresponding to the cases of W^{r N0} = 1,j, and so on, are eliminated.

To illustrate MATLAB’s DFT capabilities, consider 50 points of a 10 Hz sinusoid sampled at fs = 50 Hz and scaled by T = 1/fs.

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

In this case, the vector x contains exactly 10 cycles of the sinusoid. The fft command computes the DFT.

X = fft(x);

Since the DFT is both discrete and periodic, fft needs to return only the N⁰ discrete values contained in the single period (0 ≤ f < fs).

While Xr can be plotted as a function of r, it is more convenient to plot the DFT as a function of frequency f . A frequency vector, in hertz, is created by using N⁰ and T.

f = (0:N_0-1)/(T*N_0); stem(f,abs(X),‘k.’); >> axis([0 50 -0.05 0.55]); xlabel(‘f [Hz]’); ylabel(‘|X(f)|’);

As expected, Fig. 8.26 shows content at a frequency of 10 Hz. Since the time-domain signal is real, X(f) is conjugate symmetric. Thus, content at 10 Hz implies equal content at −10 Hz. The content visible at 40 Hz is an alias of the −10 Hz content.

Often, it is preferred to plot a DFT over the principal frequency range (−fs/2 ≤ f < fs/2). The MATLAB function fftshift properly rearranges the output of fft to accomplish this task.

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

When we use fftshift, the conjugate symmetry that accompanies the DFT of a real signal becomes apparent, as shown in Fig. 8.27.

Since DFTs are generally complex-valued, the magnitude plots of Figs. 8.26 and 8.27 offer only half the picture; the signal’s phase spectrum, shown in Fig. 8.28, completes it.

>> stem(f-1/(T*2),fftshift(angle(X)),'k.');
>> axis([-25 25 -1.1*pi 1.1*pi]); xlabel('f [Hz]'); ylabel('\angle X(f)');

Figure 8.26 |X(f)| computed over (0 ≤ f < 50) by using fft.

Figure 8.27 |X(f)| displayed over (−25 ≤ f < 25) by using fftshift.

Figure 8.28 X(f) displayed over (−25≤f <25).

Since the signal is real, the phase spectrum necessarily has odd symmetry. Additionally, the phase at ±10 Hz is zero, as expected for a zero-phase cosine function. More interesting, however, are the phase values found at the remaining frequencies. Does a simple cosine really have such complicated phase characteristics? The answer, of course, is no. The magnitude plot of Fig. 8.27 helps identify the problem: there is zero content at frequencies other than ±10 Hz. Phase computations are not reliable at points where the magnitude response is zero. One way to remedy this problem is to assign a phase of zero when the magnitude response is near or at zero.

8.7-2 Improving the Picture with Zero Padding

DFT magnitude and phase plots paint a picture of a signal’s spectrum. At times, however, the picture can be somewhat misleading. Given a sampling frequency fs = 50 Hz and a sampling interval T = 1/fs, consider the signal

$$y[n] = Te^{j2\pi \left(10\frac{1}{3}\right)n}$$

This complex-valued, periodic signal contains a single positive frequency at 10^{1 3} Hz. Let us compute the signal’s DFT using 50 samples.

>> y = T*exp(j*2*pi*(10+1/3)*n*T); Y = fft(y);
>> stem(f-25,fftshift(abs(Y)),'k.');
>> axis([-25 25 -0.05 1.05]); xlabel('f [Hz]'); ylabel('|Y(f)|');

Figure 8.29 |Y(f)| using 50 data points.

Figure 8.30 |Yzp(f)| over 5 ≤ f ≤ 15 using 50 data points padded with 550 zeros.

In this case, the vector y contains a noninteger number of cycles. Figure 8.29 shows the significant frequency leakage that results. Also notice that since y[n] is not real, the DFT is not conjugate symmetric.

In this example, the discrete DFT frequencies do not include the actual 10^{1 3} Hz frequency of the signal. Thus, it is difficult to determine the signal’s frequency from Fig. 8.29. To improve the picture, the signal is zero-padded to 12 times its original length.

>> y_zp = [y,zeros(1,11*length(y))]; Y_zp = fft(y_zp);
>> f_zp = (0:12*N_0-1)/(T*12*N_0);
>> stem(f_zp-25,fftshift(abs(Y_zp)),'k.');
>> axis([-25 25 -0.05 1.05]); xlabel('f [Hz]'); ylabel('|Y_{zp}(f)|');

Figure 8.30, zoomed in to 5 ^{≤ f ≤} 15, correctly shows the peak frequency at 10^{1 3} Hz and better represents the signal’s spectrum.

It is important to keep in mind that zero padding does not increase the resolution or accuracy of the DFT. To return to the picket fence analogy, zero padding increases the number of pickets in our fence but cannot change what is behind the fence. More formally, the characteristics of the sinc function, such as main beam width and sidelobe levels, depend on the fixed width of the pulse, not on the number of zeros that follow. Adding zeros cannot change the characteristics of the sinc function and thus cannot change the resolution or accuracy of the DFT. Adding zeros simply allows the sinc function to be sampled more finely.

8.7-3 Quantization

A B-bit analog-to-digital converter (ADC) samples an analog signal and quantizes amplitudes by using 2*^B* discrete levels. This quantization results in signal distortion that is particularly noticeable for small B. Typically, quantization is classified as symmetric or asymmetric and as either rounding or truncating. Let us investigate rounding-type quantizers.

The quantized output x^q of an asymmetric rounding converter is given as†

$$x_{\rm q} = \frac{x_{\rm max}}{2^{B-1}}\lfloor \frac{x}{x_{\rm max}}2^{B-1} + \frac{1}{2}\rfloor$$

The quantized output x^q of a symmetric rounding converter is given as

$$x_{\rm q} = \frac{x_{\rm max}}{2^{B-1}}\left(\lfloor \frac{x}{x_{\rm max}}2^{B-1}\rfloor + \frac{1}{2}\right)$$

Program CH8MP1 quantizes a signal using one of these two rounding quantizer rules and also ensures no more than 2*^B* output levels.

function [xq] = CH8MP1(x,xmax,B,method)
% CH8MP1.m : Chapter 8, MATLAB Program 1
% Function M-file quantizes x over (-xmax,xmax) using 2^b levels.
% Uses rounding rule, supports symmetric and asymmetric quantization
% INPUTS: x = input signal
% xmax = maximum magnitude of signal to be quantized
% B = number of quantization bits
% method = default 'sym' for symmetrical, 'asym' for asymmetrical
% OUTPUTS: xq = quantized signal
if (nargin<3),
   disp('Insufficient number of inputs.'); return
elseif (nargin==3),
   method = 'sym';
elseif (nargin>4),
   disp('Too many inputs.'); return
end
x(abs(x)>xmax)=xmax*sign(x(abs(x)>xmax)); % Limit amplitude to xmax
switch lower(method)
   case 'asym'
       xq = xmax/(2^(B-1))*floor(x*2^(B-1)/xmax+1/2);
       xq(xq>=xmax)=xmax*(1-2^(1-B)); % Ensure only 2^B levels
   case 'sym'
       xq = xmax/(2^(B-1))*(floor(x*2^(B-1)/xmax)+1/2);
       xq(xq>=xmax)=xmax*(1-2^(1-B)/2); % Ensure only 2^B levels

^† Large values of x may return quantized values x^q outside the 2*^B* allowable levels. In such cases, x^q should be clamped to the nearest permitted level.

otherwise
        disp('Unrecognized quantization method.'); return
end

Several MATLAB commands require discussion. First, the nargin function returns the number of input arguments. In this program, nargin is used to ensure that a correct number of inputs is supplied. If the number of inputs supplied is incorrect, an error message is displayed and the function terminates. If only three input arguments are detected, the quantization type is not explicitly specified and the program assigns the default symmetric method.

As with many high-level languages such as C, MATLAB supports general switch/case structures† :

switch switch_expr,
case case_expr,
      statements;
...
otherwise,
      statements;

end

CH8MP1 switches among cases of the string method. In this way, method-specific parameters are easily set. The command lower is used to convert a string to all lowercase characters. In this way, strings such as SYM, Sym, and sym are all indistinguishable. Similar to lower, the MATLAB command upper converts a string to all uppercase.

The floor command rounds input values to the nearest integer toward minus infinity. Mathematically, it computes ·. To accommodate different types of rounding, MATLAB supplies three other rounding commands: ceil, round, and fix. The ceil command rounds input values to the nearest integers toward infinity, ( ·”); the round command rounds input values toward the nearest integer; the fix command rounds input values to the nearest integer toward zero. For example, if x = [-0.5 0.5];, floor(x) yields [-1 0], ceil(x) yields [0 1], round(x) yields [-1 1], and fix(x) yields [0 0]. Finally, CH8MP1 checks and, if necessary, corrects large values of x^q that may be outside the allowable 2*^B* levels.

To verify operation, CH8MP1 is used to determine the transfer characteristics of a symmetric 3-bit quantizer operating over (−10,10).

>> x = (-10:.0001:10); xsq = CH8MP1(x,10,3,'sym');

>> plot(x,xsq,'k'); axis([-10 10 -10.5 10.5]); grid on;

>> xlabel('Quantizer input'); ylabel('Quantizer output');

Figure 8.31 shows the results. Clearly, the quantized output is limited to 2*^B* = 8 levels. Zero is not a quantization level for symmetric quantizers, so half of the levels occur above zero and half of the levels occur below zero. In fact, symmetric quantizers get their name from the symmetry in quantization levels above and below zero.

By changing the method in CH8MP1 from ‘sym’ to ‘asym’, we obtain the transfer characteristics of an asymmetric 3-bit quantizer, as shown in Fig. 8.32. Again, the quantized output is limited to 2*^B* = 8 levels, and zero is now one of the included levels. With zero as a quantization

^† A functionally equivalent structure can be written by using if, elseif, and else statements.

Figure 8.31 Transfer characteristics of a symmetric 3-bit quantizer.

Figure 8.32 Transfer characteristics of an asymmetric 3-bit quantizer.

level, we need one fewer quantization level above zero than there are levels below. Not surprisingly, asymmetric quantizers get their name from the asymmetry in quantization levels above and below zero.

There is no doubt that quantization can change a signal. It follows that the spectrum of a quantized signal can also change. While these changes are difficult to characterize mathematically, they are easy to investigate by using MATLAB. Consider a 1 Hz cosine sampled at fs = 50 Hz over 1 second.

x = cos(2*pi*n*T); X = fft(x); T = 1/50; N_0 = 50; n = (0:N_0-1);

Upon quantizing by means of a 2-bit asymmetric rounding quantizer, both the signal and spectrum are substantially changed.

>> xaq = CH8MP1(x,1,2,'asym'); Xaq = fft(xaq);
>> subplot(2,2,1); stem(n,x,'k'); axis([0 49 -1.1 1.1]);
>> xlabel('n');ylabel('x[n]');
>> subplot(2,2,2); stem(f-25,fftshift(abs(X)),'k'); axis([-25,25 -1 26])
>> xlabel('f');ylabel('|X(f)|');
>> subplot(2,2,3); stem(n,xaq,'k');axis([0 49 -1.1 1.1]);

Figure 8.33 Signal and spectrum effects of quantization.

>> xlabel('n');ylabel('x_{aq}[n]');

>> subplot(2,2,4); stem(f-25,fftshift(abs(fft(xaq))),'k'); axis([-25,25 -1 26]);

>> xlabel('f');ylabel('|X_{aq}(f)|');

The results are shown in Fig. 8.33. The original signal x[n] appears sinusoidal and has pure spectral content at ±1 Hz. The asymmetrically quantized signal xaq[n] is significantly distorted. The corresponding magnitude spectrum |Xaq(f)| is spread over a broad range of frequencies.