A thing or two about the Fourier Transform

Recently I have been reading up on frequency domain image processing. I am still just beginning to understand how it works. Over the last few weeks I have been trying to understand the ** Fourier Transform **. Although the gist of  Fourier Series is easy to understand from its formula, that of the Fourier Transform isn’t (to me at least). This blog post describes the Fourier Transform as I understood it, along with some helpful diagrams. I am using NumPy and matplotlib for plotting. As a reference for all the equations mentioned in this post, I am using Digital Image Processing (3rd Edition).

The Fourier Series

I haven’t done any digging into its theory, but the formula for Fourier Series is pretty straight forward. Given f has a period of T

f(t) = \sum\limits_{n=-\infty}^\infty c_{n} e^{i\frac{2 \pi n}{T}}

where

c_n = \frac{1}{T} \int\limits_{\frac{-T}{2}}^{\frac{T}{2}} f(t)e^{-i\frac{2 \pi n}{T}}

for
n = 0,\pm 1, \pm2, .....

To simplify this, we can replace \omega=\frac{2\pi}{T} . f now becomes.

f(t) = \sum\limits_{n=-\infty}^\infty  c_n e^{i \omega n}

Using Euler’s Formula

f(t) = \sum\limits_{n=-\infty}^\infty  c_n (cos(n \omega t) + i. sin(n \omega t))

We can say that, f is now exactly equal to an infinite sum of certain sines and cosines. The Wikipedia Page has some pretty nice visualizations of the concept.

The Discrete Fourier Transform

This is what took me a while to understand. The discrete fourier transform, is derived from the Fourier Transform. Given a discrete function f whose M samples have been collected, we have.

F(u) = \sum\limits_{x=0}^{M-1} f(x)e^{-i\frac{2\pi u x}{M}}

for

u = 0, 1, 2, ..., M -1

As you can see, the function f of x is transformed into the function F of u .
We can get f back by the inverse Fourier Transform

f(x) = \frac{1}{M} \sum\limits_{u=0}^{M-1} F(u)e^{i\frac{2\pi u x}{M}}
for
x = 0, 1, 2, ... , M -1

Although the relationship of f and F is mathematically proven, the underlying meaning was not immediately apparent to me. What puzzled me was, why F is computed at all, and how can it tell us anything about f . The way I looked at it to understand it better is as follows.

Let
\omega = \frac{2\pi}{M}

f(x) = \frac{1}{M} \sum\limits_{u=0}^{M-1} F(u)e^{i\omega u x}

Using Euler’s formula we have

f(x) = \frac{1}{M}[{F(0).(cos(0) + i.sin(0)) + F(1).(cos(1\omega x) + i.sin(1\omega x)) + F(2).(cos(2\omega x) + i.sin(2\omega x)) + ...} ]

From the above equation, it is clear why the Fourier Transform is important. The function F , has given us coefficients to express f as a sum of sines and cosines. I will try to demonstrate this with some code. I wrote a small helper function to help me input curves. The following function taken an image of a curve drawn in black on white and saves it into a .npy file. You will need the scikit-image library to use this.

from skimage import io, color, util
from matplotlib import pyplot as plt

pylab.rcParams['figure.figsize'] = (8.0, 8.0)

def img2curve(imgfile, arrayfile):
    """ Converts an image of a curve into numpy array and saves it"""

    img = io.imread(imgfile)
    img = color.rgb2gray(img)
    img = util.img_as_ubyte(img)
    array = np.zeros((img.shape[1], ))

    for r in range(img.shape[0]):
        for c in range(img.shape[1]):
            if img[r, c] < 128:
                array[c] = img.shape[0] - r
    np.save(arrayfile, array)

def show_curve(array):
    """ Plots real and imaginary parts of an array """

    ax = plt.subplot(111)
    ax.set_xlim(0, len(array))
    ax.set_ylim(-10, len(array))
    plt.plot(np.real(array), color='blue')
    plt.plot(np.imag(array), color='red')

Original Image

img = io.imread("curve.png")
io.imshow(img)

image

Derived Curve

img2curve("curve.png", "curve.npy")
curve = np.load("curve.npy")
show_curve(curve)

curve

If you’d like to, you can download the original image and the curve.npy file.

Computing the DFT

We will compute the discrete fourier transform using NumPy’s np.fft.fft method and observe the first 5 terms.

curve_fft = np.fft.fft(curve)
print curve_fft[:5]

curve_fft is our function F as described above. Applying the inverse transform to curve_fft should give us back curve exactly. But let’s see what happens by considering an element in F at a time. First we take only F(0) and apply the inverse transform by using np.fft.ifft.

Contribution of F(0)

t_0(x) = \frac{1}{M}[F(0).(cos(0) + i.sin(0))]

tmp = curve_fft.copy()
tmp[1:] = 0

array = np.fft.ifft(tmp)
show_curve(array)

F0

As you can see, it is a real constant shown in blue.

Contribution of F(1)

t_1(x) = \frac{1}{M}[F(1).(cos(1\omega x) + i.sin(1\omega x))]

tmp = curve_fft.copy()
tmp[2:] = 0
tmp[0:1] = 0

array = np.fft.ifft(tmp)
show_curve(array)

F1

Contribution of F(2)

t_2(x) = \frac{1}{M}[F(2).(cos(2\omega x) + i.sin(2\omega x))]

tmp = curve_fft.copy()
tmp[3:] = 0
tmp[0:2] = 0

array = np.fft.ifft(tmp)
show_curve(array)

F2

Contribution of F(0), F(1) and F(2)

t_0(x) + t_1(x) + t_2(x)

tmp = curve_fft.copy()
tmp[3:] = 0

array = np.fft.ifft(tmp)
show_curve(array)

F123

As you can see the resultant curve is trying to match the original curve. The imaginary part in red, will eventually get zeroed out when all the terms of F are considered.

Fourier Transform in action

The code snippet below generates successive images of contribution of each term in F . In nth still, the blue and red lines show the real and imaginary part of the curve so far, and the dotted lines show the contribution of F(n) .

import numpy as np
from matplotlib import pyplot as plt

curve = np.load("curve.npy")
curve_fft = np.fft.fft(curve)

result = np.zeros_like(curve, dtype=np.complex)

for i in range(curve_fft.shape[0]):


    tmp = np.zeros_like(curve, dtype=np.complex)
    tmp[i] = curve_fft[i]

    tmp_curve = np.fft.ifft(tmp)
    result += tmp_curve

    fig = plt.figure()
    ax = plt.subplot(111)
    ax.set_xlim(0, len(curve))
    ax.set_ylim(-10, len(curve))

    plt.plot(np.real(result), color='blue')
    plt.plot(np.imag(result), color='red')

    plt.title("Adding F(%d)" % i)
    plt.plot(np.real(tmp_curve), 'r--', color='blue')
    plt.plot(np.imag(tmp_curve), 'r--', color='red')

    fig.savefig('out/' + str(i) + ".png")

I ran the command

avconv -f image2 -r 1 -i %d.png -r 20 fft.mp4

to generate the final video.

As I understand DFT more, I’ll keep posting. Till then, thanks for reading.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s