Discrete Fourier Transform Frequently Asked Questions
Table of contents
DFT FAQ is intended for people who already have some understanding of DFT. It works well as a remainder or it can answer some commonly asked questions but it will not really teach you what DFT is or how to use it.
Before using this FAQ you shall know a bunch of things including: sampling, quantization, discretization, signal, discrete signal, complex numbers, Euler's formula, sinusoid, sine, cosine, complex sinusoid, sigma summing notation.
This page uses JavaScript library named MathJax to display mathematical equations. If you have JS disabled, you will only see raw latex text.
This FAQ was written by Paweł Sabat. You can contact me by mail at pawel@sabat.e90.biz
Transform is also called transformation. It is just a special kind of function.
DFT is a transform that takes a list of values and returns how much certain frequencies build up this signal. So in a sense DFT can take a piece of music and return how much certian frequencies occur in this signal.
DFT usually uses complex numbers to operate on.
DFT really consists of two transforms: forward DFT and inverse DFT.
Forward DFT takes signal (in time, space or other domain) and returns how much certain frequencies build this signal. It usually is reffered to as just DFT.
$$DFT(x) = y = (y_0, y_1, y_2, ..., y_{N-1})$$
$$y_k = \sum_{n=0}^{N-1}{x_n * e^{-j*2 \pi *k/N*n}}$$
By using Euler's formula you can also write
$$y_k = \sum_{n=0}^{N-1}{x_n * (cos(2\pi * k/N*n) - j*sin(2\pi * k/N * n))}$$
Inverse DFT takes how much certain frequencies build signal and gives back original signal.
$$IDFT(y) = x = (x_0, x_1, x_2, .., x_{N-1})$$
$$x_n = \frac{1}{N} * \sum_{k=0}^{N-1}{y_k * e^{j*2\pi * k/N * n}}$$
By Euler's formula:
$$x_n = \frac{1}{N} * \sum_{k=0}^{N-1}{y_k * (cos(2\pi * k/N * n) + j * sin(2\pi * k/N * n)}$$
Notice that there is a scaling factor $\frac{1}{N}$ for inverse DFT. According to mathematical theory it must be applied, so that forward and inverse DFT are inverses of each other, but the actual place may vary. Some implementations apply scaling factor when computing forward DFT, some during inverse DFT, and some (eg. KissFFT) do not apply scaling factor at all, so you must remember about it.
Notice that in forward DFT there is a minus sign at $e^{-j*2 \pi *k/N*n}$. Some implementations may choose to put minus sign in forward DFT, some at inverse DFT. If you are implementing your own DFT/FFT and getting strangely inversed result, check this one.
DFT is usually used with complex numbers (especially as FFT), but it can also be computed with real numbers. To use only real numbers you may set imaginary part of every sample to zero, which allows you to use existing DFT and FFT implementations. Some of them may also have functions to specifically deal with real valued signal (read the docs). You may also implement DFT using only real numbers - see Scientist and Engineer's Guide to Digital Signal Processing in links.
Frequencies used by DFT are $0/N, 1/N, 2/N, .., N-1/N$, but these are not in Hz. When you have complex sinusoid with equation $e^{j*2\pi * k/N * n}$ then this sinusoid does k full cycles over N samples. To get frequency in Hz multiply this frequency by sampling frequency:
$f_s$ - sampling frequency
$f_k$ - frequency of k-th complex sinusoid in Hz
$$f_k = \frac{k}{N} * f_s$$
DFT can be drived by discretization of Fourier Transform.
DFT can also be derived as a change of basis of N-dimensional vector space. You can see book Mathematics of the Discrete Fourier Transform listed in links to know more.
Interesting may also be term cross-correlation.
Fast Fourer Transform is an optimized implementation of DFT that is usually used to compute DFT. So FFT gives the same result as DFT, but usually does it faster. DFT computed from definition has computational complexity of $O(n^2)$ whereas FFT has $O(n\log{n})$.
Good description of FFT is in "Introduction to Algorithms" 3rd Edition by Charles E. Leiserson, Clifford Stein, Ronald Rivest, and Thomas H. Cormen, Chapter 30: Polynomials and the FFT.
Some popular implementations of FFT include:
There are many, many more, including implementations by processor manufacturers. More implementations are listed at FFTW benchmark page.
Keep in mind that some FFT implementations may only work for signals whose length is a power of 2.
Octave and Matlab have functions fft() and ifft() to compute dft.
# make some signal x x = [1,2] # compute dft for x y = fft(x) # compute inverse dft for y x2 = ifft(y)
Full description of these functions and other useful tools are in documentation under signal processing section (Signal Processing Toolbox in Matlab).
Code snippet below shows how you can play with DFT and sinusoids. Remember that vectors in Octave and Matlab are 1-indexed.
# how to get single full cycle of sinusoid with frequency f Hz # sampled with frequency fs Hz f = 32 fs = 128 x = 2*pi*f * (0:1/fs:1-1/fs); y = sin(x); # how to plot sampled sin plot(x, y); # how to compute dft of signal fourier = fft(y); # DFT N viariable is the same as length of signal N = length(y); # how to get dft frequencies for each dft number # each sample in dft corresponds to frequency (i-1)*fs/N # where i is 1-indexed sample number # fs - signal sampling frequency # N - number of dft points freqs = (0:(length(x)-1)) * fs/N; # how to plot dft with frequencies on x axis plot(freqs, abs(fourier));
Below is a simple implementation of DFT (forward and inverse) directly following DFT definition, that you can use as a reference. It gives the same values as Octave (within computation errors).
It is written by me and licensed under Creative Commons Zero license (public domain).
You can also download sourcecode file directly.
/*
Discrete Fourier Transform on complex numbers
compilation: cc -o cdft cdft.c -lm
*/
#include <stdio.h>
#include <math.h>
typedef struct cnum
{
float r; // real part
float i; // imaginary part
} cnum;
/*
Compute forward DFT for complex numbers
in - signal in time domain
out - signal in frequency domain
N - number of complex numbers in in, out
out[k] = sum(n=0,N-1): in[n] * e^(-j*2pi*k/N*n)
*/
void cdft(cnum* in, cnum* out, int N)
{
for(int n=0; n<N; ++n)
{
out[n].r = 0.0;
out[n].i = 0.0;
}
// for every frequency
for(int k=0; k<N; ++k)
{
// for every point
for(int n=0; n<N; ++n)
{
float ar = in[n].r;
float ai = in[n].i;
float br = cosf(2*M_PI*k*n/N);
float bi = -sinf(2*M_PI*k*n/N); // notice minus
// (ar + i * ai) * (br + i* bi) = (ar*br - ai*bi) + i (ar*bi + ai*br)
out[k].r += ar*br - ai*bi;
out[k].i += ar*bi + ai*br;
}
}
}
/*
Compute inverse DFT for complex numbers
in - signal in frequency domain
out - signal in time domain
N - number of points in in, out
out[n] = 1/N * sum(k=0,N-1): in[k] * e^(j*2pi*k/N*n)
*/
void icdft(cnum *in, cnum *out, int N)
{
for(int n=0; n<N; ++n)
{
out[n].r = 0.0;
out[n].i = 0.0;
}
// for every point
for(int n=0; n<N; ++n)
{
// for every frequency
for(int k=0; k<N; ++k)
{
float ar = in[k].r;
float ai = in[k].i;
float br = cosf(2*M_PI*k*n/N);
float bi = sinf(2*M_PI*k*n/N);
out[n].r += ar*br - ai*bi;
out[n].i += ar*bi + ai*br;
}
// divide out[n] by N
out[n].r /= N;
out[n].i /= N;
}
}
int main(void)
{
cnum in[5], out[5], in2[5];
// create some input signal
in[0].r = 1.0; // 1+2i
in[0].i = 2.0;
in[1].r = 3.0; // 3+4i
in[1].i = 4.0;
in[2].r = 2.0; // 2+4i
in[2].i = 4.0;
in[3].r = 3.0; // 3+7i
in[3].i = 7.0;
in[4].r = 2.0; // 2+5i
in[4].i = 5.0;
// compute forward dft
cdft(in, out, 5);
// display result of forward dft
printf("%f%+fi\n", out[0].r, out[0].i);
printf("%f%+fi\n", out[1].r, out[1].i);
printf("%f%+fi\n", out[2].r, out[2].i);
printf("%f%+fi\n", out[3].r, out[3].i);
printf("%f%+fi\n", out[4].r, out[4].i);
// compute inverse dft = restore signal from frequencies
icdft(out, in2, 5);
// display restored signal
printf("\nidft:\n");
printf("%f%+fi\n", in2[0].r, in2[0].i);
printf("%f%+fi\n", in2[1].r, in2[1].i);
printf("%f%+fi\n", in2[2].r, in2[2].i);
printf("%f%+fi\n", in2[3].r, in2[3].i);
printf("%f%+fi\n", in2[4].r, in2[4].i);
return 0;
}
Article written 20.11.2016