Discrete Fourier transform
|
|
In mathematics, the discrete Fourier transform (DFT), sometimes called the finite Fourier transform, is a Fourier transform widely employed in signal processing and related fields to analyze the frequencies contained in a sampled signal, solve partial differential equations, and to perform other operations such as convolutions. The DFT can be computed efficiently in practice using a fast Fourier transform (FFT) algorithm.
The sequence of n complex numbers x0, ..., xn−1 are transformed into the sequence of n complex numbers f0, ..., fn−1 by the DFT according to the formula:
- <math>f_j = \sum_{k=0}^{n-1} x_k e^{-\frac{2 \pi i}{n} j k} \quad \quad j = 0, \dots, n-1<math>
where e is the base of the natural logarithm, i is the imaginary unit (<math>i^2=-1<math>), and π is Pi. The transform is sometimes denoted by the symbol <math>\mathcal{F}<math>, as in <math>\mathbf{f} = \mathcal{F}(\mathbf{x})<math> or <math>\mathcal{F} \mathbf{x}<math>.
The inverse discrete Fourier transform (IDFT) is given by
- <math>x_k = \frac{1}{n} \sum_{j=0}^{n-1} f_j e^{\frac{2\pi i}{n} j k} \quad \quad k = 0,\dots,n-1.<math>
Note that the normalization factor multiplying the DFT and IDFT (here 1 and 1/n) and the signs of the exponents are merely conventions, and differ in some treatments. The only important thing is that the DFT and IDFT have opposite-sign exponents, and that the product of their normalization factors be 1/n. A normalization of <math>1/\sqrt{n}<math> for both the DFT and IDFT makes the transforms unitary, which has some theoretical advantages, but it is often more practical in numerical computation to perform the scaling all at once as above (and a unit scaling can be convenient in other ways).
(The convention of a negative sign in the exponent is often convenient because it means that <math>f_j<math> is the amplitude of a "positive frequency" <math>2\pi j/n<math>. Equivalently, the DFT is often thought of as a matched filter: when looking for a frequency of +1, one correlates the incoming signal with a frequency of −1.)
In the following discussion the terms "sequence" and "vector" will be considered interchangeable.
Properties
Completeness
The discrete Fourier transform is an invertible, linear transformation
- <math>\mathcal{F}:\mathbf{C}^n \to \mathbf{C}^n<math>
with C denoting the set of complex numbers. In other words, for any n ≥ 0, any n-dimensional complex vector has a DFT and an IDFT which are in turn n-dimensional complex vectors.
Orthogonality
The vectors exp(2πi jk/n) form an orthogonal basis over the set of n-dimensional complex vectors:
- <math>\sum_{k=0}^{n-1}
\left(e^{ 2\pi ijk/n }\right) \left(e^{-2\pi ij'k/n}\right) =n~\delta_{jj'} <math>
where δjk is the Kronecker delta.
The Plancherel theorem and Parseval's theorem
If Xj and Yj are the DFTs of xk and yk respectively then we have the Plancherel theorem:
- <math>\sum_{k=0}^{n-1} x_k y^*_k = \frac{1}{n} \sum_{j=0}^{n-1} X_j Y^*_j<math>
where the star denotes complex conjugation. Parseval's theorem is a special case of the Plancherel theorem and states:
- <math>\sum_{k=0}^{n-1} |x_k|^2 = \frac{1}{n} \sum_{j=0}^{n-1} |X_j|^2.<math>
Convolution theorem and cross-correlation theorem
The cyclic convolution x*y of the two vectors x = xj and y = yk is the vector x*y with components
- <math>(\mathbf{x*y})_k = \sum_{j=0}^{n-1} x_j y_{k-j} \quad \quad k = 0,\dots,n-1<math>
where we continue y cyclically so that
- <math>y_{-j} = y_{n-j}\quad\quad~~~~~~~~~~ j = 0, ..., n-1<math>
The discrete Fourier transform turns cyclic convolutions into component-wise multiplication. That is, if <math>z_k = (\mathbf{x*y})_k<math> then
- <math>Z_j=X_j Y_j \quad \quad~~~~~~~~~~ j = 0,\dots,n-1<math>
where capital letters (X, Y, Z) represent the DFTs of sequences represented by small letters (x, y, z). Note that if a different normalization convention is adopted for the DFT (e.g., the unitary normalization), then there will in general be a constant factor multiplying the above relation.
The direct evaluation of the convolution summation, above, would require <math>O(n^2)<math> operations, but the DFT (via an FFT) thus provides an <math>O(n\log n)<math> method to compute the same thing. Conversely, convolutions can be used to efficiently compute DFTs via Rader's FFT algorithm and Bluestein's FFT algorithm.
See also: Convolution theorem
In an analogous manner, it can be shown that if <math>z_k<math> is the cross-correlation of <math>x_j<math> and <math>y_j<math>:
- <math>z_k=(\mathbf{x\star y})_k = \sum_{j=0}^{n-1}x_j^*\,y_{j+k}<math>
where the sum is again cyclic in j, then the discrete Fourier transform of <math>z_k<math> is:
- <math>Z_k = X_k^*\,Y_k<math>
where capital letters are again used to signify the discrete Fourier transform.
Relationship to the trigonometric interpolation polynomial
The function
- <math>p(t) = \frac{f_0}{n} + \frac{f_1}{n} e^{it} + \frac{f_2}{n} e^{2it} + \cdots + \frac{f_{n-1}}{n} e^{(n-1)it}<math>
whose coefficients fj /n are given by the DFT of xk, above, is called the trigonometric interpolation polynomial of degree n − 1. It is the unique function of this form that satisfies the property: p(2πk/n) = xk for k = 0, ..., n − 1.
The DFT as a unitary transformation
With unitary normalization constants, the DFT becomes a unitary transformation. In a real vector space, a unitary transformation can be thought of as simply a rigid rotation of the coordinate system, and all of the properties of a rigid rotation can be found in the unitary DFT. Using unitary normalization, we can express the DFT as:
- <math> x_{k'} = \sum_{j=0}^{n-1}\mathcal{F}_{k'j}x_j<math>
where primed coordinates indicate the components of the vector x in the transformed space and
- <math>\mathcal{F}_{k'j}=\frac{e^{2\pi i jk'/n}}{\sqrt{n}} <math>
The orthogonality of the DFT is now expressed as an orthonormality condition:
- <math>\sum_{m=0}^{n-1}\mathcal{F}_{j'm}\mathcal{F}_{k'm}^*=\delta_{j'k'}<math>
The Plancherel theorem is expressed as:
- <math>\sum_{j=0}^{n-1}x_j y_j^* = \sum_{j'=0}^{n-1}x_{j'} y_{j'}^*<math>
which is the statement that the dot product of two vectors is preserved under a unitary DFT transformation. This means that the angle between two vectors is preserved as well, as is the length of a vector. The fact that the length of a vector is preserved is just Parseval's theorem:
- <math>\sum_{j=0}^{n-1}|x_j|^2 = \sum_{j'=0}^{n-1}|x_{j'}|^2<math>
Another way of looking at the unitary DFT is to note that in the above discussion, the DFT has been expressed as a Vandermonde matrix:
- <math>\mathcal{F} = \frac{1}{\sqrt{n}}
\begin{bmatrix}
\omega_n^{0 \cdot 0} & \omega_n^{0 \cdot 1} & \ldots & \omega_n^{0 \cdot (n-1)} \\
\omega_n^{1 \cdot 0} & \omega_n^{1 \cdot 1} & \ldots & \omega_n^{1 \cdot (n-1)} \\
\vdots & \vdots & \ddots & \vdots \\
\omega_n^{(n-1) \cdot 0} & \omega_n^{(n-1) \cdot 1} & \ldots & \omega_n^{(n-1) \cdot (n-1)} \\
\end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ \vdots \\ x_{n-1} \\ \end{bmatrix} <math>
where
- <math>\omega_n = e^{-2 \pi i/n}<math>
is a primitive nth root of unity. This matrix is a unitary matrix with
- <math>\det(\mathcal{F})=1<math>
where det() is the determinant function and
- <math>\mathcal{F}^*=\mathcal{F}^{-1}<math>
The real DFT
If x0, ..., xn−1 are real numbers, as they often are in the applications, then fj = fn−j*, where the star denotes complex conjugation. Hence, the full information in this case is already contained in the first half (roughly) of the sequence f0, ..., fn−1. In this case, the "DC" element f0 is purely real, and for even n the "Nyquist" element fn/2 is also real, so there are exactly n non-redundant real numbers in the first half + Nyquist element of the complex output f. Using Euler's formula, the interpolating trigonometric polynomial can then be interpreted as a sum of sine and cosine functions.
Generalized DFT
It is possible to shift the transform sampling in time and/or frequency domain by some real shifts a and b, respectively. This is sometimes known as a generalized DFT (or GDFT) and has analogous properties to the ordinary DFT:
- <math>f_j = \sum_{k=0}^{n-1} x_k e^{-\frac{2 \pi i}{n} (j+b) (k+a)} \quad \quad j = 0, \dots, n-1<math>
Most often, shifts of <math>1 \over 2<math> (half a sample) are used. While the ordinary DFT corresponds to a periodic signal in both time and frequency domains, <math>a=\frac{1}{2}<math> produces a signal that is anti-periodic in frequency domain (<math>f_{j+n} = - f_j<math>) and vice-versa for <math>b=\frac{1}{2}<math>. Thus, the specific case of <math>a = b = \frac{1}{2}<math> is known as an odd-time odd-frequency discrete Fourier transform (or O2 DFT). Such shifted transforms are most often used for symmetric data, to represent different boundary symmetries, and for real-symmetric data they correspond to different forms of the discrete cosine and sine transforms.
The discrete Fourier transform can be viewed as a special case of the z-transform, evaluated on the unit circle in the complex plane.
2-D transform
For digital image processing, the 2-D transform is used to find the frequency content of an image. The transform is defined as:
- <math>F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) e^{- 2 \pi i \left(\frac{u x}{M} +\frac{v y}{N} \right) } \quad \quad u = 0, \dots, M-1; v = 0, \dots, N-1<math>
and the inverse transform is defined as:
- <math>f(x, y) = \frac{1}{M N} \sum_{u=0}^{M-1} \sum_{v=0}^{N-1} F(u, v) e^{2 \pi i \left(\frac{u x}{M} +\frac{v y}{N} \right) } \quad \quad x = 0, \dots, M-1; y = 0, \dots, N-1<math>
where
- <math>f(x, y)<math> is a 2-D signal (e.g., an image) with <math>x<math> as the <math>x^{th}<math> column of <math>f<math>; and <math>y<math> as the <math>y^{th}<math> row of <math>f<math>
- <math>F(u, v)<math> is the 2-D frequency spectrum of <math>f(x, y)<math>
The form of the 2-D DFT can be simplified by using matrices
- <math>F = W f^T W<math>
where
- <math>F = \begin{bmatrix}
F(0, 0) & F(1, 0) & \ldots & F(M-1, 0) \\ F(0, 1) & F(1, 1) & \ldots & F(M-1, 1) \\ \vdots & \vdots & \ddots & \vdots \\ F(0, N-1) & F(1, N-1) & \ldots & F(M-1, N-1)
\end{bmatrix}<math>
- <math>W = \begin{bmatrix}
\omega_n^{0 \cdot 0} & \omega_n^{0 \cdot 1} & \ldots & \omega_n^{0 \cdot (n-1)} \\
\omega_n^{1 \cdot 0} & \omega_n^{1 \cdot 1} & \ldots & \omega_n^{1 \cdot (n-1)} \\
\vdots & \vdots & \ddots & \vdots \\
\omega_n^{(n-1) \cdot 0} & \omega_n^{(n-1) \cdot 1} & \ldots & \omega_n^{(n-1) \cdot (n-1)} \\
\end{bmatrix}<math> (this is the unitary matrix shown above)
- <math>f = \begin{bmatrix}
f(0, 0) & f(1, 0) & \ldots & f(M-1, 0) \\ f(0, 1) & f(1, 1) & \ldots & f(M-1, 1) \\ \vdots & \vdots & \ddots & \vdots \\ f(0, N-1) & f(1, N-1) & \ldots & F(M-1, N-1)
\end{bmatrix}<math>
Matrix form derivation
The forward transform can be reformulated into matrix notation
- <math>F(u, v) = \sum_{x=0}^{M-1} \sum_{y=0}^{N-1} f(x,y) e^{- 2 \pi i \left(\frac{u x}{M} \frac{v y}{N} \right) }<math>
- <math>F(u, v) =
\sum_{x=0}^{M-1}
\left[
\sum_{y=0}^{N-1} f(x,y) e^{- 2 \pi i \frac{v y}{N}}
\right]
e^{- 2 \pi i \frac{u x}{M} }
<math>
- <math>
\begin{matrix}
F(u, v) = \sum_{x=0}^{M-1} &
\underbrace{
\left[
\sum_{y=0}^{N-1} f(x,y) e^{- 2 \pi i \frac{v y}{N} }
\right]
} & e^{- 2 \pi i \frac{u x}{M} }
\\ & \mathcal{F}_y \{f(x, y)\} &
\end{matrix}
<math>
- <math>F_v = W f<math>
- <math>F(u, v) = \sum_{x=0}^{M-1} F_v(x, v) e^{- 2 \pi i \frac{u x}{M} }<math>
- <math>F = W F_v^{T}<math>
- <math>F = W f^T W<math>
Equation 1 is the 2-D DFT definition from above.
Equation 2 is equation 1 with <math>e^{- 2 \pi i \frac{u x}{M} }<math> pulled out of the inner summation.
Equation 3 shows that the inner summation is a 1-D DFT taken with respect to the <math>y<math>-dimension (the columns and denoted <math>\mathcal{F}_y\{f(x, y)\}<math>) of the signal <math>f(x, y)<math>.
Equation 4 is the matrix form of a 1-D DFT of the underbraced portion in equation 3.
Equation 5 is with <math>F_v(x, v)<math> substituted for the underbraced portion of equation 3. <math>F_v(x, v)<math> is the <math>v^{th}<math> frequency of the <math>x^{th}<math> column of <math>f(x, y)<math>. Collecting the column 1-D DFTs of <math>f(x, y)<math> forms the matrix <math>F_v<math>.
Equation 6 is derived after realizing that equation 4 is a 1-D DFT of <math>F_v^T<math>.
Equation 7 is simply equation 4 substituted into equation 6 and noting that <math>W = W^T<math> because <math>W<math> is symmetric.
In summation:
- A column of <math>F_v<math> represents the 1-D DFT of the <math>x^{th}<math> column of the signal <math>f(x, y)<math>. So, <math>F_v<math> can be created by combining the 1-D DFT column vectors taken with respect to each column of <math>f(x, y)<math>. Consequently, the transposition of <math>F_v<math> in <math>F = W F_v^T<math> and pre-multiplication of W is the same as taking the 1-D DFT of the rows of <math>F_v<math>.
Applications
The DFT has seen wide usage across a large number of fields; we only sketch a few examples below (see also the references at the end). All applications of the DFT depend crucially on the availability of a fast algorithm to compute discrete Fourier transforms and their inverses, a fast Fourier transform.
Signal analysis
Suppose a signal x(t) is to be analyzed. Here, t stands for time, which varies over the interval [0,T], and, in the case of a sound signal, x(t) is the air pressure at time t.
The signal is sampled at n equidistant points to get the n real numbers x0 = x(0), x1 = x(h), x2 = x(2h), ..., xn−1 = x((n−1)h), where h = T/n and n is even.
Then the discrete Fourier transform f0, ..., fn−1 is computed and the numbers fn/2 + 1, ..., fn−1 are discarded (they are redundant for real signals).
Then f0/n approximates the average value of the signal over the interval, and for every j = 1, ..., n/2, the argument (see complex number) arg(fj) represents the phase and the modulus |fj|/n represents one half of the amplitude of the component of the signal having frequency j/T (see wave).
The reason behind this interpretation is that the fj are approximations to the coefficients occurring in the Fourier series expansion of x(t). In general, the problem of using the DFT of discrete samples to approximate the Fourier transform of an infinite, continuous signal is called spectral estimation, and involves many more details than are described here. (For example, one often wants to window the data in order to reduce the distortion caused by the periodic boundary conditions implicit in the DFT.)
For an example see frequency spectrum.
Data compression
The field of digital signal processing relies heavily on operations in the frequency domain (i.e. on the Fourier transform). For example, several lossy image and sound compression methods employ the discrete Fourier transform: the signal is cut into short segments, each is transformed, and then the Fourier coefficients of high frequencies, which are assumed to be unnoticeable, are discarded. The decompressor computes the inverse transform based on this reduced number of Fourier coefficients. (Compression applications often use a specialized form of the DFT, the discrete cosine transform.)
Partial differential equations
Discrete Fourier transforms, especially in more than one dimension, are often used to solve partial differential equations. The reason is that it expands the signal in complex exponentials eikx, which are eigenfunctions of differentiation: d/dx eikx = ik eikx. Thus, in the Fourier representation, a linear differential equation with constant coefficients is transformed into an easily solvable algebraic equation. One then uses the inverse DFT to transform the result back into the ordinary spatial representation. Such an approach is called a spectral method.
Multiplication of large integers
The fastest known algorithms for the multiplication of large integers or polynomials are based on the discrete Fourier transform: the sequences of digits or coefficients are interpreted as vectors whose convolution needs to be computed; in order to do this, they are first Fourier-transformed, then multiplied component-wise, then transformed back.
Some discrete Fourier transform pairs
In the following table <math>\omega_n<math> stands for <math>\exp(-2\pi i/n)<math>
| <math>x_j=\frac{1}{n}\sum_{k=0}^{n-1}f_k \omega_n^{-jk}<math> | <math>f_k=\sum_{j=0}^{n-1}x_j \omega_n^{jk}<math> |
|---|---|
| <math>a^j\,<math> | <math>\frac{1-a^n\omega_n^{kn}}{1-a\omega_n^k}<math> |
| <math>{n-1 \choose j}\,<math> | <math>(1+\omega_n^{k})^{n-1}\,<math> |
References
- E. O. Brigham, The Fast Fourier Transform and Its Applications (Prentice-Hall, Englewood Cliffs, NJ, 1988).
- A. V. Oppenheim, R. W. Schafer, and J. R. Buck, Discrete-Time Signal Processing (Prentice-Hall, 1999).
- S. W. Smith, The Scientist and Engineer's Guide to Digital Signal Processing (http://www.dspguide.com/pdfbook.htm), (California Technical Publishing, San Diego, 2nd edition, 1999)de:Diskrete Fourier-Transformation
es:Transformada de Fourier Discreta pl:Dyskretna transformata Fouriera zh:离散傅里叶变换