8.4 The discrete Fourier transform

8.4 The discrete Fourier transform#

We now have everything we need. We substitute the discrete analysis frequencies \(\omega_k\) into our work-in-progress transform and apply the simplification \(e^{-j\omega_k n \Delta t} = e^{-2\pi j k n / N}\) from the previous section:

\[\hat{X}(\omega_k) \propto \sum_{n=0}^{N-1} x[n]\, e^{-j \omega_k n \Delta t} = \sum_{n=0}^{N-1} x[n]\, e^{-2\pi j k n / N}.\]

The result is a clean expression that depends only on the samples and the indices. This is the discrete Fourier transform.

Definition 23 (Discrete Fourier transform)

The discrete Fourier transform of a length-\(N\) signal \(x[n]\) is the length-\(N\) sequence

\[\texttt{DFT}(x)[k] \triangleq \sum_{n=0}^{N-1} x[n]\, e^{-2\pi j k n / N}, \qquad k \in \{0, 1, \ldots, N-1\}.\]

Intuitively, the DFT does exactly what the Fourier transform did, just over a finite set of frequencies. For each of the \(N\) bins \(k\) (the name for these discrete analysis frequencies), it synthesizes a phasor at \(\omega_k\), multiplies it by the signal to measure their similarity, and sums the result. We are effectively searching a finite set of bins for frequencies that resemble the signal.

Another observation is that the DFT is no longer a function of the sampling rate \(f_s\) at all. Instead, it depends only on \(N\), the number of input samples. This is perhaps counterintuitive, since we originally defined the bin frequencies \(\omega_k\) by dividing the sampling rate by \(N\). The sampling rate has not vanished, though: it is still what tells us the physical frequency, in Hz, that each bin corresponds to. To see this, it helps to work out how far apart adjacent bins are.

Definition 24 (DFT bin spacing)

The DFT bins for \(N\) input samples are evenly spaced in frequency, covering \(f_s\) over \(N\) bins. Recall that \(N = T \cdot f_s\), where \(T\) is duration in seconds. Accordingly:

\[\Delta f = \frac{f_s}{N} = \frac{f_s}{T f_s} = \frac{1}{T}.\]

This yields two equivalent forms:

\[\boxed{\; \Delta f = \frac{f_s}{N} \;} \qquad \text{and} \qquad \boxed{\; \Delta f = \frac{1}{T} \;} \qquad \text{(both in Hz).}\]

These two forms highlight some important properties of the DFT. The bin spacing \(\Delta f = 1/T\) depends only on the duration \(T\) of the analyzed segment, not on the sampling rate. Analyzing a longer stretch of audio always gives finer frequency resolution, no matter what \(f_s\) is.

The number of bins, on the other hand, is \(N = T f_s\), which grows with the sampling rate. So for a fixed duration, raising the sampling rate gives you more bins (extending the analysis up to a higher Nyquist frequency), but it does not pack the bins any closer together.

Real and imaginary parts#

Applying Euler’s formula to the definition splits the DFT into a real and an imaginary part, exactly as with the continuous transform:

\[\texttt{DFT}(x)[k] = R[k] + j\, I[k], \qquad R[k] = \sum_{n=0}^{N-1} x[n] \cos\!\left(\tfrac{2\pi k n}{N}\right), \qquad I[k] = -\sum_{n=0}^{N-1} x[n] \sin\!\left(\tfrac{2\pi k n}{N}\right).\]

As before, we usually care about the amplitude spectrum and phase spectrum, obtained by converting each complex bin to polar form:

\[A[k] = \sqrt{R^2[k] + I^2[k]}, \qquad \phi[k] = \tan^{-1}\!\frac{I[k]}{R[k]}.\]

Intuition: the “winding” view#

The following interactive example makes the “multiply by a phasor and sum” intuition concrete, in the spirit of the winding visualization from Chapter 5. Adjust the frequency of a real input sinusoid and the frequency of the probing phasor, and watch the wound-up signal and its center of mass in the complex plane. When the probe frequency matches a bin containing signal energy, the center of mass swings far from the origin:

# hide
# no-output
from IPython.utils.io import capture_output
with capture_output():
    %pip install -q plotly anywidget

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import ipywidgets as widgets
import icm_plotly

# hide
# autorun
def figure():
    fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.2)

    # left: the signal wound around the complex plane, colored by time so you
    # can watch time elapse along the wound path (dark = start, bright = end)
    fig.add_scatter(name="Wound signal", x=[0], y=[0], mode="lines+markers",
                    line=dict(color="rgba(150,150,150,0.35)", width=1),
                    marker=dict(size=5, color=[0], colorscale="Viridis",
                                cmin=0, cmax=1, showscale=False),
                    row=1, col=1)
    # left: the center of mass (proportional to the DFT bin), drawn last so it
    # sits on top of the colored path; a white outline keeps it legible
    fig.add_scatter(name="Center of mass", x=[0], y=[0], mode="markers",
                    marker=dict(size=16, color="#C41230",
                                line=dict(color="white", width=2)),
                    row=1, col=1)
    # right: the real input sinusoid
    fig.add_scatter(name="Input", x=[0], y=[0], mode="lines", row=1, col=2)

    fig.update_xaxes(range=[-1.2, 1.2], scaleanchor="y", scaleratio=1, row=1, col=1,
                     title_text="Real", fixedrange=True)
    fig.update_yaxes(range=[-1.2, 1.2], scaleanchor="x", scaleratio=1, row=1, col=1,
                     title_text="Imaginary", fixedrange=True)
    fig.update_xaxes(range=[0, 4], row=1, col=2, title_text="Time (s)", fixedrange=True)
    fig.update_yaxes(range=[-1.2, 1.2], row=1, col=2, title_text="Amplitude", fixedrange=True)
    fig.update_layout(showlegend=False)
    return fig


def controls(fig):
    full = {"description_width": "initial"}
    # the input and probe frequency sliders share exactly the same range and step
    freq_real = widgets.FloatSlider(description="Input frequency (Hz)", min=0.0, max=10.0,
                                    value=4.0, step=0.05, style=full)
    freq_probe = widgets.FloatSlider(description="Probe frequency (Hz)", min=0.0, max=10.0,
                                     value=4.0, step=0.05, style=full)
    sample_rate = widgets.IntSlider(description="Sample rate (Hz)", min=20, max=200,
                                    value=100, step=10, style=full)
    num_samples = widgets.IntSlider(description="Number of samples N", min=50, max=800,
                                    value=400, step=50, style=full)

    # redraw whenever any slider changes
    def update(freq_real, freq_probe, sample_rate, num_samples):
        n = np.arange(num_samples)
        t = n / sample_rate
        x = np.cos(2 * np.pi * freq_real * t)                 # real input sinusoid
        wound = x * np.exp(-1j * 2 * np.pi * freq_probe * t)  # multiply by the probe phasor
        com = wound.mean()                                    # center of mass (~ DFT bin)
        with fig.batch_update():
            fig.data[0].x = wound.real
            fig.data[0].y = wound.imag
            fig.data[0].marker.color = t                      # color by elapsed time
            fig.data[0].marker.cmin = 0.0
            fig.data[0].marker.cmax = float(t[-1]) if len(t) else 1.0
            fig.data[1].x = [com.real]
            fig.data[1].y = [com.imag]
            fig.data[2].x = t
            fig.data[2].y = x
            fig.layout.xaxis2.range = [0, float(t[-1]) if len(t) else 1]

    widgets.interactive_output(update, {"freq_real": freq_real, "freq_probe": freq_probe,
                                        "sample_rate": sample_rate, "num_samples": num_samples})
    return widgets.VBox([freq_real, freq_probe, sample_rate, num_samples])


icm_plotly.show(figure, controls)

Removing redundancy#

What is the “type signature” of the DFT? It takes \(N\) real samples and returns \(N\) complex numbers, so \(\texttt{DFT} : \mathbb{R}^N \to \mathbb{C}^N\). Since a computer stores each complex number as two floats (its real and imaginary parts), we could also view it as \(\mathbb{R}^N \to \mathbb{R}^{2N}\). But that feels wasteful. The DFT, like the Fourier transform, is an invertible bijection, so turning \(N\) numbers into \(2N\) numbers must involve redundancy.

Indeed it does, and the redundancy comes from the symmetry of real signals we met in Chapter 6. Because cosine is even and sine is odd, the DFT of a real signal satisfies

\[R[k] = R[N-k] \quad (\text{even}), \qquad I[k] = -I[N-k] \quad (\text{odd}).\]

So the upper half of the bins is just a mirror image of the lower half. This is the same even/odd symmetry of the amplitude and phase spectra from Chapter 6. We can tabulate it for a small example, \(N = 8\) at \(f_s = 1000\) Hz:

Table 2 DFT bins for \(N = 8\), \(f_s = 1000\) Hz. Blue marks the \(N/2 + 1\) non-redundant bins we actually need to compute; red marks the redundant upper bins, which merely mirror the lower ones.#
\(k\)	\(\blue{0}\)	\(\blue{1}\)	\(\blue{2}\)	\(\blue{3}\)	\(\blue{4}\)	\(\red{5}\)	\(\red{6}\)	\(\red{7}\)
Frequency (Hz)	0	125	250	375	500	625	750	875
Aliased (Hz)	0	125	250	375	500	\(-375\)	\(-250\)	\(-125\)
\(R[k]\)	\(\blue{a}\)	\(\blue{b}\)	\(\blue{c}\)	\(\blue{d}\)	\(\blue{e}\)	\(\red{d}\)	\(\red{c}\)	\(\red{b}\)
\(I[k]\)	\(\purple{0}\)	\(\blue{g}\)	\(\blue{h}\)	\(\blue{i}\)	\(\purple{0}\)	\(\red{-i}\)	\(\red{-h}\)	\(\red{-g}\)

Two additional optimizations appear in the table. The imaginary part vanishes at both ends, \(I[0] = 0\) and \(I[N/2] = 0\), because \(\sin(0) = 0\) and \(\sin(\pi n) = 0\) for all integer \(n\). Counting what is left, we need only the bins \(k = 0, 1, \ldots, N/2\), which is \(N/2 + 1\) complex bins, but with two of them (\(k=0\) and \(k=N/2\)) purely real-valued. That works out to exactly \(N\) real numbers to store, matching the \(N\) real inputs. The bijection is tidy after all: \(N\) samples in, \(N\) non-redundant coefficients out.

Important

For a real-valued signal of length \(N\), the DFT has only \(N/2 + 1\) non-redundant bins, spanning \(0\) to \(f_s/2\). This is exactly what NumPy’s np.fft.rfft (“real FFT”) returns, and it is what you will use in practice.

This is an important from an efficiency perspective as well. Because of this symmetry, we only have to run half as many computations! However, this does not change the asymptotic complexity of the DFT, it only improves the constant factor runtime. More on that next.

8.4 The discrete Fourier transform

Contents

8.4 The discrete Fourier transform#

Real and imaginary parts#

Intuition: the “winding” view#

Removing redundancy#