# Introduction to the FFT with examples in MATLAB

## Introduction

Fast Fourier transform (FFT) refers to any one of a family of algorithms that can compute the discrete Fourier transform (DFT) of a signal with *n* elements in *O*(*n* lg *n*) FLOPs. The DFT is a transformation that converts a signal from its original domain (e.g., time or space) into the frequeny domain.

This article introduces the theory of the DFT and FFT and gives some examples in MATLAB.

## Discrete Fourier transform

### Roots of unity

**Definition** (Roots of unity)**.** Let *n* be a positive integer. A complex number *z* is said to be an *n*-th root of unity if and only if *z ^{n}* = 1.

**Corollary.** If *z* is an *n*-th root of unity, so too is *z ^{w}* for any complex number

*w*.

*Proof*. (*z ^{w}*)

^{n}= (

*z*)

^{n}^{w}= 1

^{w}= 1. ∎

**Corollary.** Let *n* be a positive integer and 𝜔_{n} = *e*^{-2𝜋i/n}. Then, 𝜔_{0}, …, 𝜔_{n-1} are the only *n*-th roots of unity.

*Proof*. It’s easily verified that 𝜔_{0}, …, 𝜔_{n-1} define *n* distinct *n*-th roots of unity. Uniqueness follows from the fundamental theorem of algebra applied to the polynomial *z ^{n}* - 1. ∎

**Definition** (Kronecker delta)**.** The Kronecker delta 𝛿_{jj’} is defined to be 1 if and only if *j* = *j’* and 0 otherwise.

**Lemma.** Let *n* be a positive integer and *z* ≠ 1 be an *n*-th root of unity. Then, *z*^{0} + … + *z*^{n-1} = 0.

*Proof*. The sum in question is a geometric series, with closed form (1 - *z ^{n}*) / (1 -

*z*). Since

*z*is an

*n*-th root of unity, this is identically zero. ∎

**Corollary.** Let *n* be a positive integer and *j* and *j’* be nonnegative integers strictly smaller than *n*. Then,

*Proof*. If *j* = *j’*, each summand is 1. Otherwise, $\omega_n^{(j-j^\prime)}$ is an *n*-th root of unity and the desired result follows from the previous lemma. ∎

### Forward and inverse transforms

**Definition.** The conjugate transpose of a complex number *z* = *a* + *i* *b* is *z ^{*}* = a - i b. The conjugate transpose

*z*^{*}of a complex vector (or matrix)

**is obtained by taking the transpose of the vector (or matrix) and conjugating each entry.**

*z***Theorem.** Let

Then, the vectors *u*^{(0)}, …, *u*^{(n-1)} form an orthonormal basis for ℂ^{n}.

*Proof*. Using the previous corollary,

∎

The above establishes that the matrix *F*

(whose rows are the basis elements in the theorem) is unitary. The matrix *F* is called the *forward transform*. Because this matrix is unitary, its conjugate transpose *F*^{*} is its inverse (i.e., *F ^{*}F* =

*I*). As such, we call the conjugate transpose the

*inverse transform*.

*Remark*. In the above, we used the scaling factor 1/√*n* to ensure that the matrix *F* was unitary, simplifying mathematical discussion. MATLAB’s definitions of forward and inverse transforms `fft`

and `ifft`

do not use the same scaling factor (they are *F√n* and *F ^{*}/√n*, respectively) so as to avoid the cost of an extra vector-scalar multiply in the forward transform. Software libraries have their own conventions when it comes to the scaling factors for forward and inverse transforms, so it’s best to proceed carefully!

### Some immediate results

**Theorem** (Plancherel theorem)**.** Let ** x** and

**be vectors in ℂ**

*y*^{n}and denote by

**=**

*X**F*

**and**

*x***=**

*Y**F*

**their forward transforms. Then, ⟨**

*y***,**

*x***⟩ = ⟨**

*y***,**

*X***⟩.**

*Y**Proof*. ⟨** X**,

**⟩ = ⟨**

*Y**F*

**,**

*x**F*

**⟩ =**

*y*

*y*^{*}

*F*

^{*}

*F*

**=**

*x*

*y*^{*}

**= ⟨**

*x***,**

*x***⟩. ∎**

*y***Corollary** (Parseval’s theorem)**.** Let ** x** be a vector in ℂ

^{n}and denote by

**=**

*X**F*

**its forward transform. Then, ⟨**

*x***,**

*x***⟩ = ⟨**

*x***,**

*X***⟩ (i.e., the original vector and its transform have the same Euclidean norm).**

*X*### Convolution

For integers *a* and *b* with *b* nonnegative, denote by *a % b* the least nonnegative residue of *a* modulo *b* (i.e., the usual definition of the `%`

operator provided by most programming languages). For a vector ** x** with

*n*entries, we introduce the indexing convention

*x*=

_{j}*x*

_{j % n}whenever

*j*is negative or at least as large as

*n*. Subject to this convention, we have the following results:

**Definition** (Circular convolution)**.** The circular convolution of vectors ** x** and

**in ℂ**

*y*^{n}is a vector

**✳︎**

*x***with entries [**

*y***✳︎**

*x***]**

*y*_{j}=

*x*

_{0}y

_{j-0}+ … +

*x*

_{n-1}y

_{j-(n-1)}.

**Theorem** (Convolution theorem)**.** Let ⨂ denote the element-wise product, ** x** and

**be vectors in ℂ**

*y*^{n}. Then,

*F*(

**✳︎**

*x***) =**

*y**F*

**⨂**

*x**F*

**and**

*y**F*(

**⨂**

*x***) =**

*y**F*

**✳︎**

*x**F*

**.**

*y*The above theorem tells us that, for example, the convolution of vectors can be computed by

- taking their discrete Fourier transforms
and*X**Y* - multiplying these element-wise to get
=*Z*⨂*X*, and*Y* - taking the inverse Fourier transform of
.*Z*

In a subsequent section, we will prove that the discrete Fourier transform can be computed using only *O*(*n* lg *n*) FLOPs, suggesting that the above procedure is superior to naively computing the circular convolution from the formula, which requires *O*(*n*^{2}) FLOPs.

Similar results can be established for the circular cross-correlation.

### Real signals

In practice, the input to the forward transform is often a real signal (i.e., a vector in ℝ^{n}). It is useful to derive some facts about such signals.

**Theorem** (Conjugate symmetry)**.** Let ** x** be a vector in ℝ

^{n}and

**=**

*X**F*

**. Then, (**

*x**X*)

_{k}^{*}=

*X*

_{n-k}.

*Proof*. First, note that *X*^{*} = (*F* ** x**)

^{*}=

*x*^{⊺}

*F*

^{*}= ((

*F*

^{*})

^{⊺}

**)**

*x*^{⊺}where ⊺ denotes the ordinary transpose operation. Now, using the definition of

*F*and the conjugate transpose operation, it is straightforward to establish that

from which the desired result follows.

**Theorem.** Let ** x** and

**be vectors in ℝ**

*y*^{n}and

**=**

*z***+**

*x**i*

**. Further let**

*y***,**

*X***, and**

*Y***be their corresponding forward transforms. Then, 2**

*Z**X*=

_{k}*Z*+ (

_{k}*Z*

_{n-k})

^{*}and 2

*Y*=

_{k}*Z*- (

_{k}*Z*

_{n-k})

^{*}.

In other words, the DFT of two real signals ** x** and

**can be computed by**

*y*- creating a complex signal
whose real part is*z*and whose imaginary part is*x*,*y* - computing the DFT of
, and*z* - retrieving the DFTs of
and*x*using the above formulas.*y*

The proof is left as an exercise. MATLAB code for this procedure is provided below:

```
function [X, Y] = fft_real_signals (x, y)
% FFT_REAL_SIGNALS Computes the DFT of two real signals using one FFT.
z = x + 1.i * y;
Z = fft (z);
% Create Zr_conj, whose (j+1)-th component (MATLAB indices start at 1)
% is the conjugate of the (n-j+1)-th component of Z.
Zr_conj = conj ([Z(1) fliplr (Z(2:end))]);
X = (Z + Zr_conj) / 2.;
Y = (Z - Zr_conj) / 2.;
end
```

## Fast Fourier transform

The Cooley-Tukey algorithm is the most common FFT algorithm. It re-expresses the DFT of an arbitrary composite size *n* = *ab* in terms of smaller DFTs of sizes *a* and *b*.

### Radix-2 algorithm

Assuming *n* is even, the radix-2 algorithm corresponds to taking *a* = *b* = *n* / 2. The motivation comes from noting that

and hence

where ** E** and

**are the DFTs of the even and odd parts of the original signal**

*O***, respectively.**

*x*By the periodicity of the DFT and

the above can be equivalently expressed as the two equations

for nonnegative integers *k* strictly less than *n*/2.

## Image compression (with MATLAB code)

Image compression is one possible application of the FFT. It can be achieved by

- transforming the image into the frequency domain,
- dropping high-frequency components, and
- saving the result.

To view the image, one must inverse transform back to the original domain.

```
function [dft_R, dft_G, dft_B] = compress_image (image, ratio)
% COMPRESS_IMAGE Compresses an image.
%
% Inputs
% ------
% image The original image.
% ratio The compression ratio in (0, 1].
%
% Outputs
% -------
% dft_R The DFT of the red channel.
% dft_G The DFT of the green channel.
% dft_B The DFT of the blue channel.
% Normalize.
image_d = double (image) / double (max (max (max (image))));
% FFT
fft_image_d = fft2 (image_d);
dft_R = sparse (fft_image_d(:, :, 1));
dft_G = sparse (fft_image_d(:, :, 2));
dft_B = sparse (fft_image_d(:, :, 3));
% Size of image.
[m, n, ~] = size (image);
p = ceil (m / 2);
q = ceil (n / 2);
% Number of frequencies to remove first and second dimensions.
f = (1. - sqrt (ratio)) / 2;
dm = round (f * m);
dn = round (f * n);
% The new image will only have r * m * n frequencies.
i = p - dm + 1 : p + dm;
j = q - dn + 1 : q + dn;
dft_R(i, :) = 0;
dft_R(:, j) = 0;
dft_G(i, :) = 0;
dft_G(:, j) = 0;
dft_B(i, :) = 0;
dft_B(:, j) = 0;
end
```

```
function [image] = decompress_image (dft_R, dft_G, dft_B)
% DECOMPRESS_IMAGE Returns a full representation of the image.
%
% Inputs
% ------
% dft_R The DFT of the red channel.
% dft_G The DFT of the green channel.
% dft_B The DFT of the blue channel.
%
% Outputs
% -------
% image The image.
[m, n] = size (dft_R);
image = zeros (m, n, 3);
image(:, :, 1) = real (ifft2 (full (dft_R)));
image(:, :, 2) = real (ifft2 (full (dft_G)));
image(:, :, 3) = real (ifft2 (full (dft_B)));
end
```

```
% Code used to create the example.
earth = imread ('earth.jpg')
ratio = 0.15;
[dft_R, dft_G, dft_B] = compress_image (earth, ratio);
imshow (decompress_image (dft_R, dft_G, dft_B));
```