Fast Fourier Transform :: SAS/STAT(R) 13.1 User's Guide

Fast Fourier Transform

As shown in the last subsection, kernel density estimates can be expressed as a submatrix of a certain convolution. The fast Fourier transform (FFT) is a computationally effective method for computing such convolutions. For a reference on this material, see Press et al. (1988).

The discrete Fourier transform of a complex vector $\mb {z}=(z_{0},\ldots ,z_{N-1})$ is the vector $\mb {Z}=(Z_{0},\ldots ,Z_{N-1})$ , where

$Z_{j}=\sum _{l=0}^{N-1}z_{l}e^{2\pi i l j/N}, \quad j=0,\ldots ,N-1$

and i is the square root of –1. The vector $\mb {z}$ can be recovered from $\mb {Z}$ by applying the inverse discrete Fourier transform formula

$z_{l}=N^{-1}\sum _{j=0}^{N-1}Z_{j}e^{-2\pi i l j/N}, \quad l=0,\ldots ,N-1$

Discrete Fourier transforms and their inverses can be computed quickly using the FFT algorithm, especially when N is highly composite; that is, it can be decomposed into many factors, such as a power of 2. By the discrete convolution theorem, the convolution of two vectors is the inverse Fourier transform of the element-by-element product of their Fourier transforms. This, however, requires certain periodicity assumptions, which explains why the vectors $\bK$ and $\bC$ require zero-padding. This is to avoid wrap-around effects (see Press et al.; 1988, pp. 410–411). The vector $\bK$ is actually mirror-imaged so that the convolution of $\bC$ and $\bK$ will be the vector of binned estimates. Thus, if S denotes the inverse Fourier transform of the element-by-element product of the Fourier transforms of $\bK$ and $\bC$ , then the first g elements of S are the estimates.

The bivariate Fourier transform of an $N_{1}\times N_{2}$ complex matrix having $(l_{1}+1,l_{2}+1)$ entry equal to $z_{l_{1}l_{2}}$ is the $N_{1}\times N_{2}$ matrix with $(j_{1}+1,j_{2}+1)$ entry given by

$Z_{j_{1}j_{2}}=\sum _{l_{1}=0}^{N_{1}-1} \sum _{l_{2}=0}^{N_{2}-1} z_{l_{1}l_{2}} e^{2\pi i(l_{1}j_{1}/N_{1}+l_{2}j_{2}/N_{2})}$

and the formula of the inverse is

$z_{l_{1}l_{2}}=(N_{1}N_{2})^{-1} \sum _{j_{1}=0}^{N_{1}-1} \sum _{j_{2}=0}^{N_{2}-1} Z_{j_{1}j_{2}} e^{-2\pi i(l_{1}j_{1}/N_{1}+l_{2}j_{2}/N_{2})}$

The same discrete convolution theorem applies, and zero-padding is needed for matrices $\bC$ and $\bK$ . In the case of $\bK$ , the matrix is mirror-imaged twice. Thus, if S denotes the inverse Fourier transform of the element-by-element product of the Fourier transforms of $\bK$ and $\bC$ , then the upper-left $g_{X}\times g_{Y}$ corner of S contains the estimates.

The KDE Procedure

Fast Fourier Transform