Applications of the Fourier Transform

In the previous lecture notebook, we looked into detail about how the 1D FFT works in Python, and saw an example of using the FFT to detect a weak sinusoidal signal in a noisy dataset.

In this lecture notebook, you will explore the application of the 1D FFT for filtering signals, and also learn about the 2D FFT and and application of it in calculating diffraction patterns.

Learning objectives: After finishing this notebook, you should be able to:

  1. Use the FFT to filter numerical data
  2. Interpret artifacts in your filtered data arising from the FFT
  3. Calculate the 2D Fourier transform of 2D data
  4. Construct 2D images of mask patterns and calculate the far-field diffraction pattern
In [ ]:
# Notebook code 
import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append('resource/asnlib/public/')
sys.path.append('../answer_checking/') 
from validate_answers import *
In [ ]:
# Notebook code
plt.rcParams['figure.dpi'] = 100

Background information on Filtering Signals

A common way to remove noise from your data is to apply a low-pass filter:

https://en.wikipedia.org/wiki/Low-pass_filter

This is equivalent to passing an electrical signal through the following RC circuit:

To predict what the output voltage is for a given input voltage, we need to know the infinite impulse response, which is the response of the circuit to an input an input voltage that is a delta-function $\delta(0)$. For the circuit above, the impulse response is given by (see here):

\begin{equation*} h(t) = \frac{1}{\tau} e^{-\frac{-t}{\tau}} \theta(t) \end{equation*}

where $\tau = RC$ is the RC time constant of the circuit, and $\theta(t)$ is a step function. We can see what this looks like using this code:

In [ ]:
# Notebook code 
tau = 1
plt.figure(figsize=(12,3))
t = np.linspace(-5,5,1000)
h = np.exp(-t/tau)*(t>=0)
plt.plot(t,h)
plt.ylabel("Impulse response h(t)")
plt.xlabel("Time");

Once we know the impulse function of a filter, we can then calculate output $V_{out}(t)$ of the filter circuit for any given input voltage $V_{in}(t)$ using a convolution:

\begin{equation} V_{out}(t) = \int_{-\infty}^t V_{in}(t-\tau) h(\tau) d\tau \end{equation}

A nice way to understand what this does is to make a plot. The output voltage $V_{out}$ at time $t_0$ is a "weighted average" of all the input voltages $V_{in}$ at times $t < t_0$ where the "weighing factor" is given by the function $h(-\tau)$ that is decaying exponentially "backwards" in time:

In [ ]:
# Notebook code 
fix, ax = plt.subplots(figsize=(12,3))
vin = np.random.normal(size=len(t))*0.05
plt.plot(t,vin,'.', color='grey', label=r"$V_{in}$")
plt.plot(t,np.exp(t)*(t<0),'r', label=r"Impulse response")
plt.fill_between(t,np.exp(t)*(t<0), alpha=0.1, color='r')
plt.xlabel("Relative Time")
ax.get_yaxis().set_visible(False)
plt.legend()
plt.title("Illustration of the convolution integral");

We say that this filter has a "memory" because the value of $V_{out}$ depends on the values of $V_{in}$ in the past. We also say that it is "causal" because the value of $V_{out}$ at time $t$ does not depend on values of $V_{in}$ in the future.

(Aside: Any real filter you make in the lab will of course be "causal", but in software you can easily make a non-causal filter if you want...)

While it is possible to apply a low pass filter by convolving $V_{in}(t)$ with $h(\tau)$ using, for example, two nested for loops, it is far more common in practice to perform filtering in the frequency domain.

To see how this works, we first look at an important mathematical result, known as the convolution theorem. Consider now a completely general function $c(x)$ given by the convolution of two function $f(x)$ and $g(x)$:

\begin{equation} c(x) = \int_{-\infty}^{\infty} f(y) g(x-y) dy \end{equation}

If $\tilde f(\omega)$ and $\tilde g(\omega)$ are the Fourier transforms of $f(x)$ and $g(x)$ respectively, then the convolution theorem states that the Fourier transform $\tilde c(\omega)$ of $c(x)$ is given by:

\begin{equation} \tilde c(\omega) = \tilde f(\omega) \ \tilde g(\omega) \end{equation}

This is a big deal! Why? In the first equation above, for each point in $h(x)$ I need to calculate a full integral. If these are all vectors of size $N$, this integral takes $N$ steps, and so calculating $h(x)$ takes order $N^2$ steps. In the second equation, this is only point-wise multiplication of the two functions! So if I work in Fourier space, calculating the Fourier transform $\tilde c(\omega)$ takes only order $N$ steps. If $N$ is $10^6$, this is the difference between $10^6$ steps and $10^{12}$ steps, which is a big difference! If a step takes 1 ms, this is the difference between a calculation of 3.6 hours or 31 years!

(Of course, we still need to take the Fourier transform and inverse Fourier transform, but thanks to Gauss' genius, this scales only like $N\log N$, so only $6\times10^6$ steps.)

To calculate $V_{out}$ of our filter, we then need to find the Fourier transform of our impulse response function $h(t)$. For the low pass filter above, this is given by:

\begin{equation} \tilde h(\omega) = \frac{1}{1+i\omega \tau} \end{equation}

The Fourier transform $\tilde h(\omega)$ of the impulse response function is also referred to as the Transfer Function. In our case, it tells you how the amplitude and phase of an oscillating input signal is related to the amplitude and phase of the output signal. Using phasor notation $V(t) = \tilde V e^{-i\omega t}$, the input and output signals are related by:

\begin{equation} \tilde V_{out}(\omega) = \tilde h(\omega) \tilde V_{in}(\omega) \end{equation}

You can see a plot of the transfer function for an RC filter here:

In [ ]:
# Notebook code 
tau = 1
w = np.geomspace(1e-3,1e3, 100)
h = 1/(1+1j*w*tau)

plt.figure(figsize=(12,4))
plt.subplot(121)
plt.plot(w,np.abs(h))
plt.xlabel("$\omega$")
plt.ylabel("$|h(\omega)|$")
plt.xscale('log')
plt.yscale('log')
plt.grid()

plt.subplot(122)
plt.plot(w, np.angle(h)/np.pi*180)
plt.xscale('log')
plt.xlabel("$\omega$")
plt.ylabel(r"Phase angle of $\tilde h$ (degrees)")
plt.grid()
plt.yticks(np.arange(-90,1,15));

To get back to the output singal $V_{out}(t)$, all we have to do is then apply the inverse Fourier transform:

$$ V_{out}(t) = \int_{-\infty}^{\infty} \tilde V_{in}(\omega) \tilde h(\omega) e^{i\omega t} d\omega $$

In python, you would first use np.fft.fft() to calculate the FFT of the input signal, multiply this by the filter transfer function you want to use, and then take the inverse FFT using np.fft.ifft(). Done!

Filtering Exercises

You will start by calculating the result of filtering a square pulse using a low-pass filter.

Exercise 1(a): Make a square voltage pulse:

$$ V(t) = \begin{cases} 0 & \text{for } 0 <t <=50 \text{ s}\\ 1 & \text{for } 50 <t <=100 \text{ s} \\ 0 & \text{for } 100 <t <=150 \text{ s} \end{cases} $$

Your pulse should include 1000 points

In [ ]:
# t = np.linspace(...)
###
### Your code here
###


# If you are well versed with vectorized comparison operators, you
# can create v in one line! Otherwise, a for loop is OK too...
# v = ....
###
### Your code here
###


plt.plot(t,v)
plt.ylabel("Voltage (V)")
plt.xlabel("Time (s)")

answer_8_1a_1 = np.copy(t)
answer_8_1a_2 = np.copy(v)
In [ ]:
question = "answer_8_1a"
num = 2

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

Exercise 1(b): Write a function to use a FFT to apply a low-pass filter to the data. You filter should have a time constant $\tau = 5$ seconds. Make a plot of the original data along with the filtered data.

Note: When using np.fft.ifft(), it will return a complex-valued array even if the data has no imaginary component (as it should if you have done your filtering correctly). For plotting, you will need to use np.real() to convert your array to real values before plotting (if you forget, matplotlib will give you a warning messgae). It is therefore handy to have your low_pass function return the real-valued part of the filtered data.

In [ ]:
def low_pass(v, tau, dt):
    ###
### Your code here
###


vfilt = low_pass(v, 5, t[1]-t[0])

plt.plot(t,v, label="V_{in}")
plt.ylabel("Voltage (V)")
plt.xlabel("Time (s)")
plt.plot(t, vfilt, label="V_{out}");

answer_8_1b_1 = np.copy(vfilt)
In [ ]:
question = "answer_8_1b"
num = 1

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

Exercise 1(c): Use your function to apply a low pass filter with $\tau = 20$ seconds. As usual, your plot should have appropriate axis labels, and a legend.

In [ ]:
# vfilt = low_pass(....)
# plt.plot(....)
# .....
###
### Your code here
###


answer_8_1c_1 = np.copy(vfilt)
In [ ]:
question = "answer_8_1c"
num = 1

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

What is strange about your filtered function at $t=0$? What is going on?

Exercise 1(d): We can solve the problem we saw in Exercise 1(c) by using zero-padding of the array during the Fourier transform. Read the documentation page of the numpy fft to learn how to do this:

https://docs.scipy.org/doc/numpy/reference/generated/numpy.fft.fft.html

Write a new version of the low pass function low_pass2() function, this time padding this signal you are filtering with as many points as it already contains. If you do this correctly, it will eliminate the "boundary" effects from your FFT we saw above. Show this by making a plot of your low-pass filtered data after filtering with the new function.

In [ ]:
def low_pass2(v, tau, dt):
    ###
### Your code here
###


vfilt2 = low_pass2(v, 5, t[1]-t[0])
plt.plot(t,v, label="V_{in}")
plt.ylabel("Voltage (V)")
plt.xlabel("Time (s)")
plt.plot(t, vfilt2, label="V_{out}");

answer_8_1d_1 = np.copy(vfilt2)
In [ ]:
question = "answer_8_1d"
num = 1

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

2D Fourier transforms

So far, we have been considering the Fourier transform of one-dimensional data, and thinking in particular about the case of signals varying in time, for which the Fourier transform gives a decomposition of the signal in the frequency domain.

In 2D, it is more typical to take the Fourier transform of spatial data (for example, of an image in space). In this case, the units of the axis of your Fourier transform is not frequency, but instead a two-dimensional wave vector, usually denoted by letter $k$, with units of inverse distance m$^{-1}$. For a 2D FT, the wave vector is a two-dimensional vector $\bf{k}$ with vector components $k_x$ and $k_y$.

The 2D FT is used extensively in solid-state physics, and also in image analysis.

Here, we will explore two applications of the 2D Fourier transform for the calculation of Fraunhofer diffraction patterns.

Fraunhofer Diffraction

Fraunhofer diffraction considers the very-far field scattering of a colimated light off of a planar object that is much smaller than the distance to the screen (and smaller than the diameter of the beam). A classical example of 1D Fraunhofer diffraction is the single and double slit diffraction patterns. Here, we will explore the Fraunhofer diffraction patterns of 2D planar masks.

Calculating the Fraunhofer diffraction pattern is actually pretty simple: in the far field limit, the image of intensity the light on the screen is directly given by the square magnitude of the Fourier transform of the shape of the aperture.

Note that as long as the assumptions above hold, the image of the diffraction pattern on the screen is independent of the wavelength: changing from a green to red laser pointer, for example, changes only the size / scaling of the image (the dots on the wall will be further apart for a green laser than a red one).

Exercise 2(a): We will aim to create a "mask" for calculating a diffraction pattern. Your mask should consist of an array of $(N+1)\times(N+1) = 2001\times2001$ points, corresponding to $N=2000$.

To do this, we will use the meshgrid() command to create 2D arrays filled with X and Y coordinates. Have your X and Y coordinates run from -1000 to 1000 nanometers.

To understand exactly what meshgrid() does, the following code will make a plot of the two matrices size by side using subplot(121) and subplot(122) together with imshow().

In [ ]:
N = 2000

# Lower case letters correspond to the 1D array of coordinates in the two directions
x = np.linspace(-N/2,N/2,N+1)
y = np.linspace(-N/2,N/2,N+1)

# Upper case letters correspont to the 2D matrices containing the X or Y coordinates
X, Y = np.meshgrid(x, y)

# This we need to tell python what the range of the image plot are
extent = [x[0],x[-1], y[0], y[-1]]

# We will use the subplots command to make two plots beside each other
plt.subplots(figsize=(10,3.5))
plt.subplot(121)
plt.imshow(X,origin='lower',extent=extent)
plt.xlabel("X")
plt.ylabel("Y")
plt.colorbar()
plt.subplot(122)
plt.imshow(Y,origin='lower', extent=extent)
plt.xlabel("X")
plt.ylabel("Y")
plt.colorbar()
plt.tight_layout()

What do you see in the plots? What do your X and Y matrices look like?

Exercise 2(b): Use your X and Y matrices to make a mask matrix M in which all points should be 1, except for points with a radius of $r = 50$ points from the center of the array, which should be zero. The code will make two plots of your mask: one of the full matrix, and one zoomed into the range of a 200x200 pixel zoom at the center of the image. The mask value of 1 will indicate where the light will be transmitted, and a mask value of 0 will indicate where the light will be blocked, and so this corresponds to the diffraction through a circular aperture.

In [ ]:
# It seems like this is magic, but using the right matrices, it does just work!
# M =  _____ < 50
###
### Your code here
###


plt.subplots(figsize=(10,3.5))

# The full plot
plt.subplot(121)
exts = [x[0], x[-1], y[-1], y[0]]
plt.imshow(M, extent=exts)
plt.xlabel("x (nm)")
plt.ylabel("y (nm)")

# The zoom
plt.subplot(122)
zoom = 200
start = int(len(x)/2-zoom/2)
end = start + zoom
M_zoom = M[start:end,start:end]
exts_zoom = [x[start], x[end], y[start], y[end]]
plt.imshow(M_zoom, extent=exts_zoom)
plt.xlabel("x")
plt.ylabel("y")

answer_8_2b_1 = np.copy(M)
In [ ]:
question = "answer_8_2b"
num = 1

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

Exercise 2(c): Now calculate and plot the diffraction pattern from your mask. Plot a zoom in of the diffraction pattern around the center 200 pixels of the image. Your axis labels on your diffraction pattern should correspond to wavenumbers with the appropriate units. (Hint: make sure to use fftshift, which also works on 2D arrays!)

In [ ]:
# First, we take the FFT of the mask
Mt = np.fft.fft2(M)

# Now translate this into a diffraction pattern
# diffraction_pattern = ....
###
### Your code here
###


# Calculate the zoom 
zoom=100
# start=....
# end=....
###
### Your code here
###


# Extents: this one we have to think about carefully! 
# A good trick: what is the largest wave number k_max? 
# Your full image will then run from -k_max to k_max. From this,
# you can work out the extents 

# k_max = ....
# k_zoom = k_max * zoom / len(x)
# extent_zoom = [ , , , ]
###
### Your code here
###


# Now the 
plt.imshow(diffraction_pattern, extent=exts)
# plt.xlabel(...)
# plt.ylabel(...)
###
### Your code here
###

plt.colorbar()

answer_8_2c_1 = np.copy(diffraction_pattern)
answer_8_2c_2 = k_max
answer_8_2c_3 = k_zoom
answer_8_2c_4 = np.copy(extent_zoom)
In [ ]:
question = "answer_8_2c"
num = 4

to_check = [question + "_%d" % (n+1) for n in range(num)]
feedback = ""
passed = True
for var in to_check:
    res, msg = check_answer(eval(var), var)
    passed = passed and res
    print(msg); feedback += msg + "n"

assert passed == True, feedback

Exercise 2(d) Here, you will use the ipywidgets library to make an interactive zoom of your diffraction pattern. Fill in the code in the update function to generate the zoom as specified in the number of pixels around the center of the image.

After running the code, you will get an image with two sliders that allow you to adjust the zoom and the maximum of the colormap (in logarithmic intervals...). The update is not super-fast: it is best to select the slider and use the arrow keys to adjust them (using the mouse will sometimes result in the image lagging behind).

It is quite cool to zoom out and take the vmax to the lowest value in the array, particularly if you make your figure very big by using plt.figure(figsize=(10,10)) (or some size that fills your screen nicely).

In [ ]:
from ipywidgets import interact

# To really be able to tweak the colormap over a large range,
# we will make an array of maximum colors values that is 
# gemoetrically spaced
N = 20
vmax_array = np.geomspace(1e1,6e7,N)

def update(zoom=1000, v_index=N-1):
    vmax = vmax_array[v_index]
    ###
### Your code here
###

    
interact(update, zoom=(20,1000,10), v_index=(1,N-1,1))

Can you understand where the funky patterns you see come from?