4.3. Fourier Analysis¶
4.3.1. Fourier Transform and Discrete Fourier Transform¶
Consider the Fourier transform pair [1]:
\(x\) denotes the temporal or spatial coordinate, and \(\omega\) denotes the frequency coordinate. Equation (1) defines the forward Fourier transform from \(f(x)\) to \(F(\omega)\). Equation (2) defines the backward (inverse) Fourier transform from \(F(\omega)\) to \(f(x)\).
Suppose the function \(f(x)\) can be sampled in an interval \([0, A]\) with \(N\) discrete points of the same sub-interval \(\Delta x = A/N\) as:
The forward discrete Fourier transform (DFT) can be defined to be:
There is a relationship between \(F(\omega)\) (in Eq. (1)) and \(\tilde{F}(k/A)\) (in Eq. (3)), which will be derived in what follows.
Assume \(f(x) = 0\) for \(x < 0, x > A\). Equation (1) can then be rewritten as:
To facilitate the derivation, the integrand in Eq. (4) be defined as:
Aided by the trapezoid rule and Eq. (5), the integration of Eq. (4) can be approximated as:
Assume
then Eq. (6) can be written as:
Because the longest wave length that the sampling interval allows is \(A\), the frequency of the fundamental mode is
which is the spacing of the frequency-domain (\(\omega\)) grid that covers the frequency interval \([-\Omega/2, \Omega/2]\) with \(N\) points. Aided by using Eq. (8), it can be obtained that
and thus
Because
it can be shown that
Equations (9) and (10) are the reciprocity relations.
To proceed, write
Equation (7) becomes
Substituting Eq. (3) into the previous equation gives:
which defines the scaling relation between the Fourier transform (Eq. (1)) and the discrete Fourier transform (Eq. (3)).
4.3.2. Example Code¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 |
class Transform(object):
def __init__(self, ngrid, extent, average=False):
from numpy import arange, empty
from fftw3 import Plan
self.ngrid = ngrid
self.extent = extent
self.interval = interval = extent[1] - extent[0]
# calculate xgrid.
self.xgrid = xgrid = arange(ngrid, dtype='float64')
xgrid /= ngrid-1
xgrid *= interval
xgrid += extent[0]
self.dx = dx = xgrid[1] - xgrid[0]
# calculate bandwidth, kgrid, and kscale.
self.bw = bw = 1.0 / dx
self.kgrid = kgrid = arange(ngrid, dtype='float64')
kgrid /= ngrid
kgrid *= bw
kgrid -= bw/2
self.kscale = 1.0 if average else interval/2
self.kscale /= ngrid/2
# make x-/k-arrays.
self.xarrw = empty(ngrid, dtype='complex128')
self.karr = empty(ngrid, dtype='complex128')
self.karrw = empty(ngrid, dtype='complex128')
# make fftw plans.
self.wforward = Plan(self.xarrw, self.karrw,
direction='forward', flags=['estimate'])
self.wbackward = Plan(self.karrw, self.xarrw,
direction='backward', flags=['estimate'])
def forward(self):
from numpy.fft import fft, fftshift
self.karr[:] = fftshift(fft(self.xarrw))
self.wforward()
self.karrw[:] = fftshift(self.karrw)
self.karr *= self.kscale
self.karrw *= self.kscale
def report(self):
import sys
sys.stdout.write('ngrid: %d; ' % self.ngrid)
sys.stdout.write('extent: %g, %g; ' % tuple(self.extent))
sys.stdout.write('interval: %g; ' % self.interval)
sys.stdout.write('dx: %g; ' % self.dx)
sys.stdout.write('bandwidth: %g; ' % self.bw)
sys.stdout.write('krange: %g, %g ' % (self.kgrid[0], self.kgrid[-1]))
sys.stdout.write('\n')
class SineTransform(Transform):
def __init__(self, ngrid, extent, freq, **kw):
from numpy import sin, pi
super(SineTransform, self).__init__(ngrid, extent, **kw)
# remember the frequency.
self.freq = freq
# initialize x/t data.
self.xarrw[:] = sin(2*pi * freq * self.xgrid)
# for plotting.
self.fig = None
self.xax = None
self.kax = None
def plot(self, figsize=(12, 6)):
from numpy import absolute
from matplotlib import pyplot as plt
# create the figure.
self.fig = fig = plt.figure(figsize=figsize)
# plot in t/x-space.
self.xax = xax = fig.add_subplot(1, 2, 1)
xax.plot(self.xgrid, self.xarrw.real)
xax.set_title('$N$ = %d' % self.ngrid)
xax.set_xlim(self.xgrid[0], self.xgrid[-1])
xax.set_ylim(-1.1, 1.1)
xax.set_xlabel('$t$/$x$ (s/m)')
xax.grid()
# plot in f/k-space.
self.kax = kax = fig.add_subplot(1, 2, 2)
kax.plot(self.kgrid, absolute(self.karr), label='numpy.fft.fft')
kax.plot(self.kgrid, absolute(self.karrw), label='fftw3.plan')
kax.set_xlim(self.kgrid[0], self.kgrid[-1])
kax.set_xlabel('$f$/$k$ (Hz/$\\frac{1}{\mathrm{m}}$')
kax.grid()
kax.legend()
class RectTransform(Transform):
def __init__(self, ngrid, extent, **kw):
from numpy import absolute, sinc
super(RectTransform, self).__init__(ngrid, extent, **kw)
# initialize x/t data.
self.xarrw.fill(0)
self.xarrw[absolute(self.xgrid) < 0.5] = 1
self.kana = sinc(self.kgrid)
# for plotting.
self.fig = None
self.xax = None
self.kax = None
def plot(self, figsize=(12, 6)):
from numpy import absolute
from matplotlib import pyplot as plt
# create the figure.
self.fig = fig = plt.figure(figsize=figsize)
# plot in t/x-space.
self.xax = xax = fig.add_subplot(1, 2, 1)
xax.plot(self.xgrid, self.xarrw.real)
xax.set_title('$N$ = %d' % self.ngrid)
xax.set_xlim(self.xgrid[0], self.xgrid[-1])
xax.set_ylim(-0.1, 1.1)
xax.set_xlabel('$t$/$x$ (s/m)')
xax.grid()
# plot in f/k-space.
self.kax = kax = fig.add_subplot(1, 2, 2)
kax.plot(self.kgrid, absolute(self.karr), label='numpy.fft.fft')
kax.plot(self.kgrid, absolute(self.karrw), label='fftw3.Plan')
kax.plot(self.kgrid, absolute(self.kana), label='analytical')
kax.set_xlim(self.kgrid[0], self.kgrid[-1])
kax.set_xlabel('$f$/$k$ (Hz/$\\frac{1}{\mathrm{m}}$')
kax.grid()
kax.legend()
def main():
from matplotlib import pyplot as plt
stfm = SineTransform(2**7, (-1.5, 1.5), 1.0, average=True)
stfm.report()
stfm.forward()
stfm.plot()
rtfm1 = RectTransform(2**5, (-1., 1.), average=True)
rtfm1.report()
rtfm1.forward()
rtfm1.plot()
rtfm2 = RectTransform(100, (-5., 5.))
rtfm2.report()
rtfm2.forward()
rtfm2.plot()
plt.show()
if __name__ == '__main__':
main()
|
-
class
pyengr.fourier.
Fourier
(ngrid, extent, average=False)¶ Fourier transform pair that supports both
numpy.fft
andfftw3.Plan
.
[1] | William L. Briggs and Van Emden Henson, The DFT: An Owners’ Manual for the Discrete Fourier Transform, SIAM, 1987. http://www.amazon.com/gp/product/0898713420/ |