. Fourier Methods

Paul Abbott

Pre-requisites

To improve the appearance of graphics, TraditionalForm will be used as the default format type for text that appears in graphics,

$FormatType = TraditionalForm;

and the default style for text in graphics is set as follows.

$TextStyle = {FontSize -> 10};

A consistent set of options for ListDensityPlot and ListPlot3D is desirable.

SetOptions[ListDensityPlot, Mesh -> False, 
        Frame -> False, PlotRange -> All];

SetOptions[ListPlot3D, Mesh -> False, 
        PlotRange -> All];

You should have already met Fourier series (see e.g., http://mathworld.wolfram.com/FourierSeries.html) on a number of occasions. Fourier series are expansions of any periodic function, f(x), in terms of an infinite of the form

f(x) = 
      ((Subscript[a, 0])/2) +  (∑_(n = 1)^m)((Subscript[a, n]) cos(
                    n x) + (Subscript[b, n]) sin(
                    n x)).

.. Orthogonality

In this section the the basic orthogonality integrals will be examined.

(∫_0^(2  π))sin(n x) sin(
          m x) dx

In this section you will prove that

(∫_0^(2  π))sin(n x) sin(
            m x) dx = π (Subscript[δ, n, m])

You should now have the desired result.

\(TraditionalForm
    \`∫\_0\%\(2  π\)\(sin(n\ x)\)\ \(cos(m\ x)\) 
        dx\)

\(TraditionalForm
    \`∫\_0\%\(2  π\)\(cos(n\ x)\)\ \(cos(m\ x)\) 
        dx\)

\(TraditionalForm\`\(1\/L\) \(∫\_\(-L\)\%L\( 
            cos(\(m\ π\ x\)\/L)\)\ \(cos(\(n\ π\ x\)\/L)\) dx\)\)

.. Square wave

The Fourier series of a function assumes that the function is periodic outside its specified domain. Consider the function defined by

f(x)={({{1,  0 ≤ x ≤ π}, {-1, π ≤ x ≤ 2  π}}),

with the additional (periodicity) condition

f(x)=f(x mod (2 π)),

which causes the function to repeat itself outside the domain [0, 2  π].

f(
        x_ /;  0 ≤ x ≤ π) :=  ...

f(x_ /;   ... ) :=  ...

f(x_):=f(Mod[..., ...])

.. Fourier coefficients

From the orthogonality integrals computed above it can be shown that the Fourier coefficients are:

\(TraditionalForm\`a\_0 = \(1\/π\) \(∫\_0\%\(2  π\)\(f(
              x)\)\ dx\), 
a\_n = \(1\/π\) \(∫\_0\%\(2  π\)\(f(x)\)\ \(cos(
              n\ x)\) dx\), \ \ n ≥ 1, 

    b\_n = \(1\/π\) \(∫\_0\%\(2  π\)\(f(x)\)\ \(sin(
              n\ x)\) dx\), \ \ n ≥ 
      1. \)

.. Partial sum of the Fourier series to m terms

Using (Subscript[a, 0]), (Subscript[a, n]), and (Subscript[b, n]) computed in §.. the sum to m terms found in the previous section can be written as

FourierSeries[m_] := 
      ((Subscript[a, 0])/2) +  (∑_(n = 1)^m)((Subscript[a, n]) cos(n x) + (Subscript[b, n]) sin(
                  n x))

.. Lanczos σ-factors

The Gibbs phenomenon (see e.g. http://mathworld.wolfram.com/GibbsPhenomenon.html) is the overshoot that occurs at each discontinuity in Fourier Series — you should be able to see this in the graphics of §...  Note also that, in general, Fourier series are slowly convergent. The Lanczsos σ-factor (see Lanczsos 1956 and Hamming 1973),

(Subscript[σ, m_, n_]) = ((sin(((n π)/m)))/(((n π)/m)));

improves the convergence of Fourier series. Including (Subscript[σ, m, n]), the modified partial sum to m terms is

σSeries[m_] := 
      ((Subscript[a, 0])/2) +  (∑_(n = 1)^m) (Subscript[σ, m, n]) ((Subscript[a, n]) cos(
                    n x) + (Subscript[b, n]) sin(n x))

References

R W Hamming Numerical Methods for Scientists and Engineers, 2nd ed., McGraw-Hill, New York, 1973.

C Lanczsos, Applied Analysis, Prentice-Hall, Chapter 4, 1956.

.. Triangular wave

. DFT and FFT

In this section you will be introduced to the basics of the discrete Fourier transform (DFT) and fast Fourier transform (FFT).

As an aside, there are a number of other transform functions relevant including the discrete cosine transform (DCT), used in JPEG image compression, the Hartley transform, and the discrete wavelet transform (DWT). Wavelets are finitely restricted waves (wave-packets) which permit orthogonal expansion of functions which are both spatially and temporally restricted. Recently, wavelets have been adopted by the FBI as their image compression standard for fingerprints.

.. DFT and FFT

Read the discussion on fast transform algorithms to see why the FFT is faster than the DFT.

.. One-dimensional Example

Analysis of a Periodic Signal

Discrete Fourier Basis

One way of visualising the discrete Fourier basis is to take the inverse (discrete) Fourier transform of unit vectors.

Filter

𝒹 = 𝓈 + 
        Table[Random[Real, {-(1/2), (1/2)}], {x, 256}];

We can use the FT to construct a "filter" which removes the (high-frequency) random "noise" from the signal.

.. Two-dimensional Example

Centred Fourier Transform

It was pointed out in §.. that the Fourier transform, displayed as of the function of (angular) frequency, is “wrapped”, i.e., the low frequency terms appear at both ends of the plot and the high frequency terms are in the middle of these plots.  Fourier can compute the FFT in n-dimensions directly. Using list operations, it is simple to re-arrange the data so that the low frequency terms are in the middle and the high frequency terms are near the edge. Here is the centred FFT.

FT(f_) := 
        RotateLeft[
          Fourier[f], ⌊Dimensions[f]/2⌋];

and the corresponding inverse transform:

IFT(F_) := 
        InverseFourier[
          RotateRight[
            F, ⌊Dimensions[F]/2⌋]];

Point Aperture

Consider the following dataset:

𝒹 = 
        Table[If[
            15 ≤ x ≤ 18 ∧ 
              15 ≤ y ≤ 18, 1.0, 0.0], {x, 33}, {y, 
            33}];

Slit

.. Simple Image Processing

Letter E

Enter the following code for generating a rectangle:

Box(x_,y_,xw_,yw_,(({{x0_, y0_}, {x1_, y1_}}))):=1 /; 
          x0 - xw ≤ x ≤ x1 + xw ∧ 
            y0 - yw ≤ y ≤ y1 + yw

Box(x_,y_,xw_,yw_,(({{x0_, y0_}, {x1_, y1_}})))=0;

Using Box it is easy to construct many shapes including the letter ℰ:

ℰ=Transpose[N[Table[Box[x,y,1,1,(({{0, -5}, {0, 5}}))]+Box[x,y,1,1,(({{3, -5}, {5, -5}}))]+Box[x,y,1,1,(({{3, 0}, {5, 0}}))]+Box[x,y,1,1,(({{3, 5}, {5, 5}}))],{x, -9, 15},{y, -12, 12}]]];

Section of the transformed function

The FT arises in many areas of Optics. For example, holography and spatial filtering both rely on the FT. You may know that, taking a region from a hologram, the whole image can be reconstructed, albeit at a lower resolution.  In this example, you will sample the FT of ℰ  (corresponding to taking a region from a hologram) and reconstruct the image (using the inverse FT).

Select the following region of frequency space:

𝓇  =  
        Transpose[
          Take[Transpose[Take[𝒻,  {5,  20}]],  {5,  
              20}]];

Crude sampling

This image can be improved somewhat by rejecting all image (pixel) values that are less than a certain cutoff value.

Another section of the transformed function

Select another region of frequency space:

𝓇 = 
      Take[Take[𝒻, 10]^T, 10]^T;

.. Convolution

In this section a “noisy” signal will be examined using a low-frequency filter using the technique of convolution.

Noisy signal

{n, m} = Dimensions[ℰ]

𝓃 = 
        Table[Random[Real, {-0.1, 0.1}], {n}, {m}];

𝓃ℰ = 𝓃 + ℰ;

ListDensityPlot[𝓃ℰ];

ListPlot3D[𝓃ℰ];

ListPlot3D[𝓃ℰ /. {x_ /; x < 0.2 -> 0, 
            x_ /; x > 0.8 -> 1}];

Since the signal has a peak value of 1, the signal to noise ratio, S/N, is about 10. Note that the signal is no longer very clear. Nevertheless, the eye+brain is very good at recognising the pattern from the noise.

Filter

The following code produces a Gaussian (low-pass) filter suitable for smoothing the signal:

σ(k_) := 
        Table[N[(E^(-k (((((j/m)))^2) + ((((i/n)))^2))))], {i, n}, {j, m}];

Modifying the factor k leads to different (width) filters.

ListDensityPlot[σ(3)]

Convolution

A common and important application of the FT is doing convolutions. The basic idea is that integrals of the form g(
          x) = ∫f(y) k(y - x) dy can be evaluated simply by multiplying the FTs of f and k.  Convolutions are used extensively, for example, in smoothing data.

The following code takes the product of the FTs, removes the small imaginary components using Chop and then applies InverseFourier to the result:

Σ(k_) := 
        InverseFourier[
          Chop[Fourier[𝓃ℰ] Fourier[σ(
                  k)]]];

In effect this performs the integral convolution, ∫f(y) k(
              y - x) dy.

ListDensityPlot[\[LeftBracketingBar]Σ(
            1000)\[RightBracketingBar]];

ListPlot3D[\[LeftBracketingBar]Σ(
                1000)\[RightBracketingBar]];

ListPlot3D[\[LeftBracketingBar]Σ(
                10)\[RightBracketingBar]];

Convolve

.. Edge detection

A straightforward image processing problem is that of edge detection, whereby sharp intensity transitions are enhanced relative to slow transitions.  An edge can be approximated by a step-function.  Derivatives of step-functions (also called Heaviside functions) lead to the Dirac delta “function” which will clearly reveal sharp transitions.

Approximating a step function

Enter the following useful approximation to the step function:

ℋ(x_, 
        k_ : 5) = ((2 (tan^(-1))(k x))/π);

The syntax k_:  5 defines a default argument.  If no second argument is supplied to ℋ then it will default to 5.

Laplacian

The extension of this derivative technique for enhancing edges in two (or more) dimensions is to apply the (discrete) Laplacian operator to the (pixel) image function f(x, y). In our case f(x, y) corresponds to the matrix of intensity (pixel) values stored in ℰ. Note that the x and y indices refer to the x - y element of the matrix.

Recall that the (continuous) Laplacian operator is defined by Δ=(((∂^2))/(∂(x^2)))+(((∂^2))/(∂(y^2))), implemented in Mathematica as

Δ:=(((∂^2) #)/(∂(x^2)))+(((∂^2) #)/(∂(y^2)))&

At this time it might be a good idea to review section 3.3.1 on numerical derivatives. Noting that

Δ(f(x, y))-((f(x - h, y) + f(h + x, y) + f(x, y - h) + f(x, h + y) - 
              4 f(x, y))/((h^2)))+O[h]^3

we now have a matrix operator:

δ=(({{0, 1, 0}, {1, -4, 1}, {0, 1, 0}}));

which approximates the action of the continuous Laplacian operator.  For this operator, the geometrical relation in the x - y plane relating the f[x, y] values can be illustrated as follows: take the direct product (i.e., not the Dot product) between Δ and the values of f[x, y] discretized on a regular grid

(δ/(h^2)) (({{f(x - h, h + y), f(x, h + y), f(h + x, h + y)}, {f(x - h, y), f(x, y), f(h + x, y)}, {f(x - h, y - h), f(x, y - h), f(h + x, y - h)}}))

Adding all the elements together we recover the truncated series approximation:

Together[Plus @@ Plus @@ %]

To compute the action of the discrete Laplacian we need to convolve (i.e., compute the direct product of) the matrix δ with each 3×3 block of the matrix generated by discretizing f(x, y). The Partition function is useful here:

ℳ = 
      Table[(Subscript[f, i, j]), {i, 6}, {j, 5}]

𝓅 = 
      Partition[ℳ, {3, 3}, {1, 1}]

We can now compute the direct product of each block with δ as follows,

Map[δ #  &, 𝓅, {2}]

and then add up all entries in each sub-matrix.

Apply[Plus, %, {2, 3}]

In Mathematica version 4.0, this operation is built-in as ListConvolve:

ListConvolve[δ,  ℳ]

One should see that the edges of the image have been enhanced.

Since the Laplacian also arises in Poisson's equation and Schrödinger's equation, this discretized numerical approximation is very useful.

. Sound

And now for some fun …  Mathematica has the ability not only to Plot functions but to Play them as well.

.. Simple Sounds

For those that are curious, the lowest note on the piano is A sounding at 27.50 Hz and the highest is the top C sounding at 4186.009 Hz. A440 concert pitch refers to tuning the A above middle C to 440.000 Hz. See http://www.pianoworld.com/pitch.htm.

Play[sin(2 π 256(t^2)),{t, 0, 2}];

To make question  clearer, consider the following plot:

f(t_) = 2 (t^2) - t + 1;

Plot[sin(f(t)), {t, 0, 5}, 
      PlotStyle -> RGBColor[1, 0, 0]];

What is the instantaneous frequency of sin(f(t))? One way of answering this is to consider the periodic function, sin(f'((Subscript[t, 0])) t-(f'((Subscript[t, 0])) (Subscript[t, 0])-f((Subscript[t, 0]))))≡sin((Subscript[ω, (Subscript[t, 0])]) t - (Subscript[ϕ, (Subscript[t, 0])])), that has the same frequency as sin(f(t)) at t = (Subscript[t, 0]). For example, with (Subscript[t, 0]) = 3:

Plot[Evaluate[sin(f'(3) t-(f'(3) 3-f(3)))],{t, 0, 5},PlotStyle -> RGBColor[0, 0, 1]];

Show[%, %%];

From this it can be seen that the (angular) frequency of sin(f(t)) at t = (Subscript[t, 0]) is ω=2 πν=f'((Subscript[t, 0])) and the phase is ϕ=f((Subscript[t, 0]))-f'((Subscript[t, 0])) (Subscript[t, 0]). Hence the instantaneous frequency of sin(2 π f (t^2)) is

(((∂(2 π f (t^2)))/(∂t))/(2 π))

and the rate of increase of frequency is 2  f t Hz per second.

.. Touch-tone Phone

It turns out that all the 12 tones on a touch-tone phone (i.e., 1 - 9, 0, * and #) are generated by forming the superposition of two tones of differing frequencies, played for 0.1 seconds.

Forming pairs

The numbers 0-9

TouchTone(
        n_Integer /; 0 ≤ n ≤ 9) := 
      PlayTouchTone(
        tones[[n + 1]])

Multiple digits