Learn to Use the Discrete Fourier Transform

• Matlab
• Basics
• Tutorials

Discrete-time sequences arise in many ways:  a sequence could be a signal captured by an analog-to-digital converter; a series of measurements; a signal generated by a digital modulator; or simply the coefficients of a digital filter.  We may wish to know the frequency spectrum of any of these sequences.  The most-used tool to accomplish this is the Discrete Fourier Transform (DFT), which computes the discrete frequency spectrum of a discrete-time sequence [1 – 3].  The DFT is easily calculated using software, but applying it successfully can be challenging.  This article provides Matlab examples of some techniques you can use to obtain useful DFT’s.

Review:  The DFT Formula

For a discrete-time sequence x(n), the DFT is defined as:

$$X(k)=\sum_{n=0}^{N-1} x(n)e^{-j2\pi kn/N} \qquad (1)$$

where

X(k) = discrete frequency spectrum of time sequence x(n)

N = number of samples of x(n) and X(k)

n = 0: N-1 = time index

k = 0: N-1 = frequency index

Equation 1 calculates a single spectral component or frequency sample X(k).  To find the whole spectrum over k = 0 to N-1, Equation 1 must be evaluated N times.

We see that, by definition, the DFT applies to a finite-length sequence of N samples.  Equation 1 does not contain variables for time and frequency, but uses time and frequency indices n and k instead. The frequency index is sometimes referred to as “frequency bins.” For sample time of T_s, the discrete time variable is given by:

t = nT_s                                 (2)

For sample frequency f_s = 1/T_s, the discrete frequency variable is given by:

f = k*f_s/N                      (3)

While x(n) is normally a real sequence, X(k) is complex.  For real x(n), the real part of X(k) is even with respect to f = f_s/2, and the imaginary part is odd.

Example 1.  Simple Sequence

This very simple example will let us get our feet wet.  Let x(n) be a time sequence containing 8 samples, as follows:

x = [1 2 3 4 3 2 1 0]

Here is the Matlab code to find the DFT.  The Matlab function fft(x,N) finds the N-point DFT using a Fast-Fourier Transform algorithm [4].  The algorithm provides an efficient calculation of the DFT, provided that N is a power of 2.

    N= 8;                    % number of time samples
    n= 0:N-1;                % time index
    x= [1 2 3 4 3 2 1 0];    % time sequence
    k= 0:N-1;                % frequency index
    X= fft(x,N)              % DFT of x

    X = [16  -4.83–j4.83  0  .828–j.828  0  .828+j.828  0  -4.83+j4.83]

We should note, that although indices n and k are 0, 1, … 7, Matlab indexes the vectors x and X starting with 1:  e.g. x(1), x(2), … x(8).  In this article, we’ll plot vectors as if they are indexed starting with 0.  The time sequence is plotted in Figure 1a, and the complex DFT is plotted vs k in Figure 1b.  Figure 1a was plotted using stem(n,x).  As shown, the DFT is N = 8 points long, with k = 0: N-1 = 0: 7.

We are usually interested in the magnitude of X(k), that is:

$$ |X(k)|= \sqrt{Re\{X(k)\}^2 + Im\{X(k)\}^2}  \qquad (4) $$

This is computed in Matlab using the command abs(X).  Suppose our x(n) has a sample frequency of 100 Hz.  We can calculate the DFT magnitude and the frequency vector (Equation 3) using the following Matlab code:

    Xmag= abs(X);       % magnitude of X
    fs= 100;            % Hz sample frequency
    f= k*fs/N;          % Hz frequency vector (Eq 3)

Before plotting the magnitude spectrum, recall that the spectrum over f = 0 to f_s of any sampled signal includes an image above f_s/2 that is symmetric with respect to the spectrum below f_s/2.  Xmag is plotted vs f in Figure 2.  The spectrum is indeed symmetric, although we do not compute the value for f = f_s (k = N).  If we were to evaluate Equation 1 for k = N, we would find that X(N) = X(0).  Given the spectrum’s symmetry, we often choose to plot X(f) over 0 to f_s/2 instead of 0 to f_s.

The spectrum plot of Figure 2 contains only 8 frequency samples.  It is quite easy to generate more frequency samples for a more detailed plot, as we’ll see in Example 2.

Figure 1.  a) Discrete-time sequence x(n).  b) DFT X(k)

Figure 2.  Magnitude of the DFT, f_s = 100 Hz.

Spectrum Magnitude at k = 0 (f = 0)

From Equation 1, the value of X(k) for k = 0 is:

$$X(0)=\sum_{n=0}^{N-1} x(n) \qquad (5)$$

X(0) is real and is the dc component of the spectrum.  For the sequence x = [1 2 3 4 3 2 1 0] of Example 1, Equation 5 yields X(0) = 16.  This value appears in the spectral magnitude plot of Figure 2.  For a sequence having X(0) as the largest component, it is often convenient to scale the amplitude such that X(0) = 1.  This is done by dividing x by its sum:  in this case, the scaled sequence would be:

x = [1 2 3 4 3 2 1 0]/16               (6)

Example 2.  Zero Padding

Equation 3 states that f = k*f_s/N.  Thus, the frequency step between spectrum samples is:

$$ \Delta f= f_s/N \qquad (7) $$

We can reduce the frequency step by increasing the number of time samples N in x(n).  This can be done by appending zeros to x(n), an operation called zero padding.  Instead of using x = [1 2 3 4 5 6 7 0], we could use x = [1 2 3 4 5 6 7 0 0 0 0 …..].   The following Matlab code repeats Example 1 with two changes:  we append 56 more zeros (for a total length N = 64), and we scale the amplitude of x(n) by 1/16 as was done in Equation 6.  Sample frequency is still 100 Hz.

    N= 64;                  % number of time samples
    n= 0:N-1;               % time index
    x= [1 2 3 4 3 2 1 0 zeros(1,56)]/16;    % sequence with zero pad
    k=0:N-1;                % frequency index
    X= fft(x,N);            % DFT of x
    Xmag= abs(X);           % magnitude of X
    fs= 100;                % Hz sample frequency
    f= k*fs/N;              % Hz frequency vector

(Note that the Matlab fft(x,N) function automatically zero-pads x to length N.  Thus, we could simply have used x = [1 2 3 4 3 2 1]/16 in the above code).  The zero padded version of x(n) is shown in Figure 3 (left), and the DFT is plotted vs f in Figure 3 as dots (right).  The frequency step is now 100 Hz/64 = 1.5625 Hz.  Compare this to the DFT without zero padding of Figure 2, which has a frequency step of 100 Hz/8 = 12.5 Hz.  The solid line shown in Figure 3 is the Continuous Fourier Transform (CFT) of x(n).  We see that zero padding has filled-in more detail of the CFT.  The amplitude of the spectrum is 1.0 at f = 0, due to the scaling factor of 1/16.

Figure 3.  Left:  x(n) with zero pad and amplitude scaling of 1/16.  Right:  |X(k)|.

We can also view our example sequence x(n) as an FIR filter’s coefficients/impulse response                h(n)= [1 2 3 4 3 2 1]/16.  Then Figure 3 (right) is the filter’s magnitude response |H(k)|.  We see that the zero-padding technique allows us to create a detailed plot of |H(k)|.

In most cases, it’s useful to plot the spectrum’s magnitude using a decibel (dB) amplitude scale.  We define the dB-magnitude spectrum as follows:

$$ XdB(k)= 10*log_{10}(|X(k)|^2)  \qquad(8) $$

or

$$ XdB(k)= 20*log_{10}|X(k)|  \qquad(9) $$

The Matlab code for Equation 9 is:

    XdB= 20*log10(Xmag);        % dB magnitude of X

XdB is shown in Figure 4, where we have plotted the spectrum just over f = 0 to f_s/2.

Figure 4.  dB-magnitude response XdB, f = 0 to f_s/2.

Example 3.  Sinewave

From Equation 1, the DFT, by definition, applies to a finite-length sequence of N samples.  It is possible to apply the DFT to a longer sequence truncated to N samples, such as a periodic sequence.  But truncating without any other conditioning usually gives unacceptable results.

Let’s consider a sine wave sampled at 128 Hz, with N = 128 samples.  Let the sine frequency be 15 Hz.  Then the sequence’s length in cycles is exactly 15.  Here is the Matlab code to generate the signal and find the DFT:

    fs= 128;               % Hz sample frequency
    Ts= 1/fs;              % s sample time
    N= 128;                % number of samples
    f0= 15;                % Hz frequency of sine
    n= 0:N-1;              % time index
    x= 2/N*sin(2*pi*f0*n*Ts);    % discrete-time sinusoid
    X = fft(x,N);           % DFT
    Xmag= abs(X);           % magnitude of X
    XdB= 20*log10(Xmag);    % dB magnitude
    k= 0:N-1;               % frequency index
    f= k*fs/N;              % Hzfrequency vector

The sinewave and its one-sided spectrum magnitude plotted over 0 to f_s/2 are shown in Figure 5.  The factor of 2/N multiplying the sine causes the spectral peak to occur at 1.0 = 0 dB.  As desired, the spectrum over f = 0 to fs/2 contains a single component at 15 Hz.  We got this fortunate result because the sinewave contained an integer number of cycles.

Now consider a 15.7 Hz sine wave sampled at 128 Hz, with N = 128 samples.  The sequence’s length in cycles is 15.7.  The Matlab code to compute the DFT is the same as above, except:

    f0= 15.7;        % Hz frequency of sine

The resulting sinewave and spectrum magnitude are shown in Figure 6.  Unfortunately, the spectrum has significant components at other than 15 or 16 Hz – this phenomenon is called spectral leakage.  We could say the leakage occurs because the DFT bins occur at … 14, 15, 16, … Hz, and there is no bin at 15.7 Hz.  Or, we can realize that we have truncated an infinite sequence, keeping N samples.  This is equivalent to multiplying the samples by a rectangular function of length N and amplitude 1.  Basically, the DFT is modulating the sinewave with this undesired rectangular function.  We call this function a rectangular window function.  In the next example, we’ll reduce spectral leakage by introducing a non-rectangular window function.

Figure 5.  15 cycles of a sinewave at 15 Hz and its one-sided spectrum magnitude, f_s = 128 Hz, N = 128.

Figure 6.  15.7 cycles of a sinewave at 15.7 Hz and its one-sided spectrum magnitude, f_s = 128 Hz, N = 128.

Example 4.  Windowing

Given that a rectangular window causes the spectral leakage shown in Figure 6, it is worth asking if we can substitute another function that reduces the leakage.  Let’s try the Blackman window function [5], defined as:

$$win(n)=0.42-0.5cos(\frac{2\pi n}{N}) + .08cos(\frac{4\pi n}{N}),\; n=0:N-1 \qquad (10) $$

The following Matlab code is the same as that of the previous example, except for two lines added after the calculation of x.  Here, we compute a Blackman window function of length N and multiply each sample of x by the corresponding sample of the window function, using the Matlab vector multiply operator (.*).  The function blackman(N,’periodic’)’ computes window values according to Equation 10.  (Note:  the transpose operator ‘ converts the blackman column vector to a row vector).

    fs= 128;                % Hz sample frequency
    Ts= 1/fs;               % s sample time
    N= 128;                 % total number of samples
    f0= 15.7;               % Hz frequency of sine
    n= 0:N-1;               % time index
    x= 2/N*sin(2*pi*f0*n*Ts);        % discrete-time sine (row vector)
    win= blackman(N,’periodic’)';    % window function (row vector)
    x_win= x.*win;           % windowed sinewave
    X = fft(x_win,N);        % DFT
    Xmag= abs(X);            % magnitude of X
    XdB= 20*log10(Xmag);
    k= 0:N-1;                % frequency index
    f= k*fs/N;               % Hz frequency vector

Figure 7 shows the sinewave x (top), the window function (middle), and the product of x and the window function (bottom).  We can say that x is modulated by the window function.  The Blackman window is one of many possible window functions.  I discussed window functions in earlier posts here [6] and here [7].

Figure 8 shows the DFT of the windowed sinewave (f₀ = 15.7 Hz).  Recalling that multiplication (modulation) in the time domain is equivalent to convolution in the frequency domain, we see that the spectrum is a convolution of a single spectrum line at f₀ with the window’s spectrum.  The window’s spectrum is shown in Figure 9, along with the spectrum of the rectangular window.  Given the more confined spectrum of the Blackman window, we would expect a better result for the Blackman-windowed DFT.  And in fact, compared to Figure 6, the spectrum of Figure 8 is closer to the ideal sine spectrum.  The peak amplitude is lower than 0 dB due to windowing.

The width of the spectral component depends on the number of samples N in the window and the sequence x; we can reduce the width (increase resolution) by increasing N.  (Try changing N to 512 in the above code).  Note that we cannot gain resolution by zero padding.  Zero padding merely fills-in detail of the window’s translated spectrum.

The window magnitude responses of Figure 9 were computed by taking the DFT of the zero-padded window functions.  We learned the techniques involved in Example 2.  (Note that the rectangular window function is just ones(1,128) and the spectrum shown is the zero-padded DFT of win = ones(1,128)/128).

Finally, keep in mind that windowing can be used in finding the spectrum of any truncated signal, not just sinewaves.

Figure 7.  Top:  Sinewave at 15.7 Hz, N = 128, f_s = 128 Hz.  

Middle:  Blackman window function.

Bottom:  Sinewave multiplied by window function.

Figure 8.  Spectrum of 15.7 Hz sinewave using Blackman window.  f_s = 128 Hz, N = 128.

Figure 9.  Magnitude responses of Blackman window function and Rectangular

window function, f_s = 128 Hz, window length = 128.

Summary

Here’s a summary of the DFT techniques we’ve learned.  In Example 1, we found the DFT of a short sequence of N = 8 time samples.  The magnitude spectrum contained N = 8 frequency samples (Figure 2).  In Example 2, we used zero padding to compute a more detailed spectrum with 64 frequency samples (Figure 3).

By definition (Equation 1), the DFT applies to a finite-length time sequence of N samples.  In Example 3, we found the DFT of an infinite sequence (sinewave) truncated to N = 128 time samples.  When the sine frequency did not fall on a frequency bin, the spectrum was degraded by spectral leakage (Figure 6).  To reduce leakage, in Example 4 we used windowing of the time sequence.  The improved spectrum is shown in Figure 8.  Windowing replaced the rectangular window (inherent in truncation of the time sequence) with a window having a more compact spectrum (Figure 9).

Additional Topics

For a DFT X(k), The power spectrum or power spectral density (PSD) is given by:

$$ P_{bin}= \frac{2}{N^2}|X(k)|^2 \quad watts/bin \qquad (11) $$

Equation 11 applies to a one-sided (0 to f_s/2) spectrum, and we assume a load resistance of 1 ohm.  I discussed the power spectrum in an earlier post [8].  PSD averaging is used for signals with a random component, such as digitally modulated signals or any signal containing Gaussian noise.  In an earlier post, I presented a Matlab function (based on the pwelch function) that computes the PSD with averaging [9].  I also discussed PSD averaging here [10].  Welch’s classic paper [11] discusses the most popular method of PSD averaging.

I have discussed windowing briefly in this article; I gave more detailed presentations in two earlier posts here [6] and here [7].  Some windows of merit, in addition to Blackman, are Hanning, Kaiser, Chebyshev, and Flattop [12].

Given that x(n) is a sampled signal, its DFT is periodic [13].  Our computation of the DFT normally covers only one period, but it’s easy to extend the summation range of Equation 1 to more than one period, as shown in an earlier post [14].

References

1.  Lyons, Richard, Understanding Digital Signal Processing, 3^rd Ed., Prentice Hall, 2011, Ch 3.

2.  Jones, Douglas and Selesnick, Ivan, The DFT, FFT, and Practical Spectral Analysis, Rice University, 2007, https://www.scribd.com/document/485235793/the-dft-fft-and-practical-spectral-analysis-2-1

3.  Mitra, Sanjit, Digital Signal Processing, 2^nd Ed., McGraw Hill, 2001, Section 3.2

4.  Matlab website, “fft”,https://www.mathworks.com/help/matlab/ref/fft.html#buuutyt-5

5.  Oppenheim, Alan and Shafer, Ronald, Discrete-Time Signal Processing, Prentice Hall, 1989, p. 448.

6.  Robertson, Neil, “The Discrete Fourier Transform and the Need for Window Functions,”, DSPrelated.com, November, 2021, https://www.dsprelated.com/showarticle/1433.php

7.  Robertson, Neil, “Evaluate Window Functions for the Discrete Fourier Transform,” DSPrelated.com, December, 2018, https://www.dsprelated.com/showarticle/1211.php

8.  Robertson, Neil, “The Power Spectrum,” DSPrelated.com, October, 2016, https://www.dsprelated.com/showarticle/1004.php

9.  Robertson, Neil, “A Simplified Matlab Function for Power Spectral Density,” DSPrelated.com, March, 2020, https://www.dsprelated.com/showarticle/1333.php

10.  Robertson, Neil, “Use Matlab Function pwelch to Find Power Spectral Density – or Do It Yourself,” DSPrelated.com, January 2019, https://www.dsprelated.com/showarticle/1221.php

11.  Welch, Peter, “The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms,” IEEE Trans. Audio and Electroacoustics, vol. AU-15, p 70 – 73, June 1967.

12.  Matlab Website, “Windows,” https://www.mathworks.com/help/signal/ug/windows.html

13.  Lyons, Richard, Understanding Digital Signal Processing, 3^rd Ed., Prentice Hall, 2011, section 3.14.

14.  Robertson, Neil, “Demonstrating the Periodic Spectrum of a Sampled Signal Using the DFT,” DSPrelated.com, March, 2019, https://www.dsprelated.com/showarticle/1235.php

                      September, 2024                       Neil Robertson

Sign in