Page 462 - Numerical Methods for Chemical Engineering

Fourier transforms in multiple dimensions 451

We assume periodicity outside of the domain,
f (x + l x 2P x , y + l y 2P y ) = f (x, y) l x,y = 0, ±1, ±2,... (9.87)

and thus determine F (q)at q [m,n] , where m = 1, 2,..., N x , n = 1, 2,..., N y :
(m − 1)π [n] (n − 1)π
[m]
[m]
[m,n]
[n]
q = q x e x + q y e y q x = q y = (9.88)
P x P y
The maximum resolvable frequencies and the frequency resolutions are
N j π N j π
q j,max = = ( q j ) q j = j = x, y (9.89)
2P j 2 P j
The computed F(q [m,n] )have q-periodicity,
F(q x + l x 2q x,max , q y + l y 2q y,max ) = F(q x , q y ) l x,y = 0, ±1, ±2,... (9.90)

The F(q [m,n] ) are obtained by quadrature of (9.85),
( x)( y) N y [m] [n]
N x
[m,n]
F(q ) ≈ f (x j , y k )e −i q x x j +q y y k (9.91)
(2π)
j=1 k=1
In terms of the 2-D discrete Fourier transform values F mn , this becomes
N y
( x)( y) [m] [n]
$ % N x
[m,n]
F q ≈ F mn F mn = f jk e −i q x x j +q y y k (9.92)
(2π)
j=1 k=1
Defining W x = e −i( q x )( x) , W y = e −i( q y )( y) ,wehave
(
N x N y
(n−1)(k−1) (m−1)( j−1)
F mn = f jk W W (9.93)
y x
j=1 k=1
Thus, we obtain the 2-D discrete Fourier transform by first computing the 1-D discrete
Fourier transforms over y with x fixed at each x j ,
N y
(n−1)(k−1)
F k (x j ) = f (x j , y k )W (9.94)
y
k=1
and then performing a 1-D discrete Fourier transform over x,
N x
(m−1)( j−1)
F mn = F k (x j )W (9.95)
x
j=1
Thus, multidimensional discrete Fourier transforms can be obtained from recursive appli-
cation of the 1-D FFT algorithm.

FFT in multiple dimensions in MATLAB

In MATLAB, the multidimensional discrete Fourier transform and its inverse are computed
using the definitions above by fftn and ifftn respectively. A separate routine fft2 is provided
for the 2-D case. The code provided below computes the 2-D power spectrum of the function
f (x, y) = cos(x) + cos(3x) + 2 sin(x) cos(2y) + cos(4y) (9.96)