4.3. Fourier Analysis

4.3.1. Fourier Transform and Discrete Fourier Transform

Consider the Fourier transform pair [1]:

(1)\[\newcommand{\defeq}{\buildrel {\text{def}}\over{=}}\]\[\begin{split}F(\omega) = \mathcal{F}\left\{f(t)\right\}\defeq \int_{-\infty}^{\infty} f(x) e^{-i2\pi\omega x} dx \\\end{split}\]
(2)\[f(x) = \mathcal{F}^{-1}\left\{F(\omega)\right\} \defeq \int_{-\infty}^{\infty} F(\omega) e^{i2\pi\omega x} d\omega\]

\(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:

\[f_n \defeq f(x_n) = f(n\Delta x), \quad n = 0, \ldots, N-1\]

The forward discrete Fourier transform (DFT) can be defined to be:

(3)\[\tilde{F}(\frac{k}{A}) \defeq \sum_{n=0}^{N-1} f(n\Delta x) e^{-i2\pi\frac{nk}{N}}\]

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:

(4)\[F(\omega) = \int_{0}^{A} f(x) e^{-i2\pi\omega x} dx\]

To facilitate the derivation, the integrand in Eq. (4) be defined as:

(5)\[g(x) \defeq f(x) e^{-i2\pi\omega x}\]

Aided by the trapezoid rule and Eq. (5), the integration of Eq. (4) can be approximated as:

(6)\[F(\omega) \cong \frac{\Delta x}{2}\left[ g(0) + 2\sum_{n=1}^{N-2}g(x_n) + g(A) \right]\]

Assume

\[g(0) = g(A)\]

then Eq. (6) can be written as:

(7)\[F(\omega) \cong \Delta x \sum_{n=0}^{N-1}g(x_n) = \frac{A}{N}\sum_{n=0}^{N-1} f(x_n) e^{-i2\pi\omega x_n}\]

Because the longest wave length that the sampling interval allows is \(A\), the frequency of the fundamental mode is

(8)\[\Delta\omega = \frac{1}{A}\]

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

\[\Omega = N\Delta\omega = \frac{N}{A}\]

and thus

(9)\[A\Omega = N\]

Because

\[\Delta x = \frac{A}{N}, \quad \Delta\omega = \frac{1}{A}\]

it can be shown that

(10)\[\Delta x\Delta\omega = \frac{1}{N}\]

Equations (9) and (10) are the reciprocity relations.

To proceed, write

\[x_n\omega_k = (n\Delta x)(k\Delta\omega) = \frac{nA}{N}\frac{k}{A} = \frac{nk}{N}\]

Equation (7) becomes

\[F(\frac{k}{A}) \cong \frac{A}{N} \sum_{n=0}^{N-1} f(n\Delta x) e^{-i2\pi\frac{nk}{N}}\]

Substituting Eq. (3) into the previous equation gives:

(11)\[F(\frac{k}{A}) \cong \frac{A}{N}\tilde{F}(\frac{k}{A})\]

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 and fftw3.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/