The Discrete Fourier Transform: Frequency axes

Contents

Background

The DFT defined earlier is agnostic to sample rate — i.e., the concept of `frequency’ is in a relative sense, normalised by the total number of samples, N. However, when dealing with real-world audio and acoustics signals it can be helpful to interpret the DFT bins in terms of either cyclic frequency (cycles per second, aka Hz), normalised cyclic frequency, or normalised radian frequency.

A couple of foundational points:

  • Define a sample rate F, which tells us how many times per second we sample the signal x[n]. For example, F = 48000 implies a sample rate of 48kHz, or 48000 samples per second.
  • By extension, F tells us the sample period T (aka the time delay between consecutive samples) via T = 1/F. T is measured in seconds (per sample).

The DFT in terms of cyclic frequency (Hz)

Let us define

(1)   \begin{align*} f_k = k\frac{F}{N}, \hspace{0.5cm} t_n = nT \end{align*}

Here, f_k is in units of cycles per second (i.e., Hz), and t_n is in units of seconds. Each value of f_k is referred to as a DFT analysis frequency, since it sets the frequency of the given complex exponential phasor onto which we project x[n], as per equation 3.

The units of {f}_k are understood as cycles per second (Hz) via a units analysis, whereby

(2)   \begin{align*} {f}_k = k\frac{F}{N} \hspace{0.5cm} \frac{\mbox{[samples/second]}}{\mbox{[samples]}} \rightarrow \frac{\mbox{[cycles (implied)]}}{\mbox{[second]}} \end{align*}

and where we note that k is a dimensionless index. The `cycles’ term is implied by the context in which {f}_k is used, i.e., to determine the repetition rate of a complex exponential phasor, since the top half of the ratio in equation 2 is dimensionless.

This allows us to write the DFT in the form

(3)   \begin{align*} X[f_k] = \sum_{n=0}^{N-1} x[n] e^{-2\pi i f_{k} t_{n}}, \hspace{0.5cm} k = 0, 1, 2, ... , N-1 \end{align*}

Notice the following:

  • Inserting the expressions for f_k and t_n into equation 3 returns the original DFT expression of equation ??.
  • For k=0 we have f_0 = 0, which corresponds to the 0 Hz (aka `DC’, aka `zero frequency’) complex exponential.
  • For k=1 we have f_1 = \frac{F}{N}, which is the frequency of the first non-zero frequency bin, and by extension, the frequency resolution of the N-point DFT as a whole, i.e., the spacing between each pair of bins, expressed in Hz.
  • f_{k} spans the range from 0:F (strictly speaking, from 0:F-df, since it needs to stop 1 bin below F). Note that some writers denote this as f_k \in [0:F), where `[‘ specifies inclusion of 0, and `)’ the exclusion of F (i.e., stopping 1 bin below), and the suffix _k reminds us that the range is discrete.
  • Since the DFT operates on a finite length, discrete time signal, the available analysis frequencies f_k are said to be band limited (i.e., they span the range from [0:F), and discrete valued (i.e., there are exactly N of them). The band limiting property is interesting, because strictly speaking we can define f_k in terms of any range of values that span [0:F). For example, the range f_k \in [-F/2,F/2) is often used, in which the `upper half’ of the analysis frequency range is mapped to negative frequency.

When writing your own DSP code, such as in Matlab, all of the above makes it straightforward to produce a frequency axis in Hz for a given use of the DFT. Here is a simple example, written in Matlab:

N     = 6400;                % Length of discrete time signal
x     = 2*(rand(N, 1)-0.5);  % Dummy discrete time signal
F     = 48000;               % Sample rate
T     = 1/F;                 % Sample period
df    = F/N;                 % Frequency resolution (in Hz)
fVec  = [0:df:F-df];         % Frequency vector (spans full bandwidth, in Hz)
x_dft = fft(x);              % Transform the signal into the discrete frequency domain

plot(fVec,abs(x_dft))        % Plot the DFT magnitude against frequency (in Hz)

The DFT in terms of normalised cyclic frequency (cycles per sample)

Another useful way to write the DFT equation is in terms of a normalised cyclic frequency, \hat{f}, with units of cycles per sample. In this case, \hat{f} defines the number of cycles of phase swept out by a given complex exponential phasor per tick of the sample clock.

We define \hat{f_k} as

(4)   \begin{align*} \hat{f_k} = \frac{k}{N} \end{align*}

The units of \hat{f}_k are understood as cycles per sample via a units analysis, whereby

(5)   \begin{align*} \hat{f}_k = \frac{k}{N} \hspace{0.5cm} \frac{\mbox{[cycles (implied)]}}{\mbox{[samples]}} \end{align*}

and where we note that k is a dimensionless index. The `cycles’ term is implied by the context in which \hat{f}_k is used, i.e., to determine the repetition rate of a complex exponential phasor, since the top half of the ratio is dimensionless.

This permits a DFT expressed as

(6)   \begin{align*} X[\hat{f}_k] = \sum_{n=0}^{N-1} x[n] e^{-2\pi i \hat{f_k} n}, \hspace{0.5cm} k = 0, 1, 2, ... , N-1 \end{align*}

The normalised frequency, \hat{f_k} spans the range from 0:1 (strictly speaking, from 0:\frac{N-1}{N} since it needs to stop 1 bin below the normalised frequency limit). As earlier, we can denote this as \hat{f_k} \in [0:1).

Notice the following:

  • As expected \hat{f}_0 = 0 corresponds to the zero-frequency complex exponential.
  • \hat{f}_{N-1} defines a complex exponential that almost, but not quite, completes one complete cycle on the complex plane for each tick of the sample clock. If we imagine further increasing k to specify \hat{f}_N = 1, we see that the complex exponential in equation 6 reduces to e^{-2\pi i n}, which always evaluates to 1 for integer n. The corresponding X[\hat{f}_N] then simply returns the sum of x[n] over its N values, which is equivalent to projection onto the zero-frequency complex exponential. In other words, \hat{f}_N is the alias of \hat{f}_0, and provides no new information about x[n].

Here is a simple example, written in Matlab, that illustrates the creation of a normalised cyclic frequency axis:

N     = 6400;                  % Length of discrete time signal
x     = 2*(rand(N, 1)-0.5);    % Dummy discrete time signal
fVec  = (1/N)*[0:N-1]';        % Frequency vector (spans 0:1, in cycles/sample)
xdft  = fft(x);                % Transform the signal into the discrete frequency domain

plot(fVec,abs(xdft))           % Plot the DFT magnitude against normalised frequency

The DFT in terms of normalised radian frequency (radians per sample)

Yet another useful way to write the DFT equation is in terms of a normalised radian frequency, \hat{\omega}, with units of radians per sample. In this case, \hat{\omega}_k defines the number of radians of phase swept out by a given complex exponential phasor per tick of the sample clock.

We define \hat{\omega} as

(7)   \begin{align*} \hat{\omega}_k = k\frac{2\pi}{N} = 2\pi \hat{f}_k \end{align*}

where we first used \hat{f}_k in equation 4. The units of \hat{\omega}_k are understood as radians per sample via a units analysis, whereby

(8)   \begin{align*} \hat{\omega}_k = k\frac{2\pi}{N} \hspace{0.5cm} \frac{\mbox{[radians]}}{\mbox{[samples]}} \end{align*}

and where we note that k is a dimensionless index.

This allows us to express the DFT as

(9)   \begin{align*} X[\hat{\omega}_k] = \sum_{n=0}^{N-1} x[n] e^{-i \hat{\omega}_k n}, \hspace{0.5cm} k = 0, 1, 2, ... , N-1 \end{align*}

The normalised radian frequency, \hat{\omega}_k spans the range from 0:2\pi (strictly speaking, from 0:2\pi\frac{N-1}{N} since it needs to stop 1 bin below the normalised frequency limit). As above, we can denote this by \hat{\omega}_k \in [0:2\pi).

Here is a simple example, written in Matlab, that illustrates the creation of a normalised radian frequency axis:

N     = 6400;                  % Length of discrete time signal
x     = 2*(rand(N, 1)-0.5);    % Dummy discrete time signal
fVec  = (2*pi/N)*[0:N-1]';     % Frequency vector (spans 0:2*pi, in radians/sample)
xdft  = fft(x);                % Transform the signal into the discrete frequency domain

plot(fVec,abs(xdft))           % Plot the DFT magnitude against normalised radian frequency