The function $f(t)$ is transformed to a new function $g(\alpha)$ through the definite integral. A similarity to the dot product of functions is evident in this form and this operation can be thought of as a mapping or projection of $f(t)$ into a different independent variable $\alpha$. Existence, integrability and inversion of integral transform operations are important in the study of this topic, although not covered in these notes.
Two examples of integral transforms, the Laplace and Fourier, are discussed in this lecture. It is typical to use the Laplace transform to remove the time dependence from Fick's second law in diffusion problems. The Fourier transform is used in the study of diffraction under certain conditions.
To assist in this lecture some special symbols in Python and sympy are reviewed:
In [ ]:
import sympy as sp
sp.init_printing(use_latex=True)
In [ ]:
# symbols we will need below
x,y,z,t,c = sp.symbols('x y z t c')
# note the special declaration that omega is a positive number
omega = sp.symbols('omega', positive=True)
In [ ]:
sp.I**2
The natural logarithm of $e$ is $1$:
In [ ]:
sp.log(sp.E)
In SymPy there are two ways to deal with integration. If you would like to represent an unevaluated integral, you can use the Integral function. If you want to compute the integration of an expression you can use the integrate function.
In [ ]:
sp.Integral(sp.E**(sp.I*omega*t),t)
In [ ]:
# 'omega', positive=True
sp.integrate(sp.E**(sp.I*omega*t),t)
Where we assume there is no zero frequency (as we are dividing by $\omega$) - hence the assumption positive=True in the symbol definition above. (Try replacing $\omega$ with $y$ and inspect the expression returned by integrate.)
As the domain of the periodicity increases, the frequency spectrum required required to represent the function becomes more finely divided. Recall the argument of the trigonometric terms in the functions of the Fourier series:
$$ \frac{n \pi (\omega +c)}{d} $$where n is the order of the frequency component, c the offset relative to the origin, and d the domain width. If we let the domain width go to infinity (implying that the function is not periodic) then an integral sum is required rather than a discrete summation. The, infinte, non-periodic function and its frequency spectrum are related by the Fourier transform defined by:
$$ \hat{f}(\omega) = \sqrt{\frac{1}{2\pi}} \int^{+\infty}_{-\infty} f(t) \exp[-i \omega t] dt $$This results in a mapping of the function f(t) into frequency space.
The real or complex and even or odd nature of the function $f(t)$ determines if the transformed function is even, odd, real, or complex. For the purposes of materials crystal structures in this lecture we will be using even and real functions.
A useful physical problem requiring use of the Fourier transform is diffraction. In this problem we will use a top-hat function to represent the location of an infinity of wave sources from an aperture. We use the sp.Piecewise function to generate a "tophat" function for the Fourier transform.
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = 8, 4
p = sp.Piecewise((0,x<-1),(1,x<1),(0,True))
sp.plot(p);
At some distance from the aperture we place a detector that measures the combined intensity of all the wave sources, however due to the finite width of the slit each wave travels a different distance to the detector. The phase difference between the waves at the detector is given by the Fourier transform of the aperture function when the Fraunhofer approximation is valid.
This aperture function is even and real so we expect our transformed function to also be even and real. We use the definition of the integral transform above to write an explicit integral statement of the Fourier transform of the top-hat function above. The integral is $1$ between $c$ and $-c$ and zero elsewhere - so we can integrate just the non-zero part. This is integrated as:
In [ ]:
sp.Integral(1*sp.exp(-sp.I*2*omega*x),(x,-c,c))
Calling explicitly for the integration and assigning the result to a:
In [ ]:
a = sp.sqrt(1/(2*sp.pi))*sp.integrate(1*sp.exp(-sp.I*2*omega*x),(x,-c,c))
a
This does not (at first glance) appear to be a real function due to the two exponential terms, but we can use some of the algebraic methods built into SymPy to help. We can ask for this form using sines and cosines with the rewrite method. Furthermore - we can simplify it further with the expand function. Trial and error may be required to determine the best combination and ordering of algebraic manipulations.
In [ ]:
solution = sp.expand(a.rewrite(sp.sin))
solution
Here we can use the subs (substitution) method to set the value of c. I plotted the square of the function since the intensity of a diffracted wave is related to the time averaged energy transferred by the wave. This is proportional to the amplitude squared. As our function is real valued, we can take a shortcut and just plot the square.
In [ ]:
sp.plot(solution.subs(c,1));
In [ ]:
sp.plot(solution.subs(c,1)**2);
In [ ]:
compositeIntegral = sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),(x,1,2)) + \
sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),(x,-2,-1))
compositeIntegral
In [ ]:
om = compositeIntegral.doit()
om
The diffracted intensity from this pair of slits would appear as:
In [ ]:
sp.plot(om.rewrite(sp.sin).expand()**2)
Or we could functionalize our function to explore other parameters:
In [ ]:
def diffractionFunction(d=4.0, w=1.0):
result = sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),\
(x,-(d+w),-(d-w))) + \
sp.sqrt(1/(2*sp.pi))*sp.Integral(1*sp.exp(-sp.I*2*omega*x),\
(x,(d-w),(d+w)))
return result.doit()
In [ ]:
sp.expand(diffractionFunction(10.,2.).rewrite(sp.sin))
Perform the Fourier transformation on an odd or complex valued function. Plot the real and imaginary parts of both the target function and the transformed functions.
Solve for the diffracted intensity in two dimensions from a circular aperture. It may be easier to do this as a discrete problem using the DFT below.
The discrete Fourier Transform is defined here and is regarded as one of the most important advances in computing science in the 20th century. Other resources such as Numerical Recipes, the Python help files and many other websites detail the calculation and implementations.
It is often instructive to review other implementations of the DFT to help you gain experience. I will be modeling this implementation after Jake Vanderplas' blog article here. Following the notion in the blog article:
Forward DFT:
$$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~k~n~/~N}$$Inverse DFT:
$$x_n = \frac{1}{N}\sum_{k=0}^{N-1} X_k e^{i~2\pi~k~n~/~N}$$In this section of the notebook, we use Vanderplas' description and implementation.
For simplicity, we'll concern ourself only with the forward transform, as the inverse transform can be implemented in a very similar manner. Taking a look at the DFT expression above, we see that it is nothing more than a straightforward linear operation: a matrix-vector multiplication of $\vec{x}$,
$$\vec{X} = M \cdot \vec{x}$$with the matrix $M$ given by
$$M_{kn} = e^{-i~2\pi~k~n~/~N}$$With this in mind, we can compute the DFT using simple matrix multiplication as follows:
In [ ]:
import numpy as np
def DFT_slow(x):
"""Compute the discrete Fourier Transform of the 1D array x"""
x = np.asarray(x, dtype=float)
N = x.shape[0]
n = np.arange(N)
k = n.reshape((N, 1))
M = np.exp(-2j * np.pi * k * n / N)
return np.dot(M, x)
We can use the "all close" function to check if the result from DFT_slow and Numpy are close:
In [ ]:
x_signal = np.random.random(1024)
np.allclose(DFT_slow(x_signal), np.fft.fft(x_signal))
I think it would be instructive to symbolically expand the matrix above so that it is clear how n*k leads to a two dimensional matrix. Switching to sympy symbols to expose the details we can do the following:
In [ ]:
import sympy as sp
from sympy import Matrix
import numpy as np
sp.init_printing()
x is the input vector.k is the wavenumber or frequency.n is the component of the input vector.
In [ ]:
x = sp.Matrix(sp.symbols('x0:5'))
n = sp.Matrix(sp.symbols('n0:5')).T
k = sp.Matrix(sp.symbols('k0:5'))
N = sp.symbols('N')
M = (-sp.I*2*sp.pi*k*n/N).applyfunc(sp.exp)
In [ ]:
M*x
Each frequency element is projected into each point of the input vector - the matrix links k and n. So - the contribution at each point is a sum of each frequency contribution, similar to the dot product of functions.
In [ ]:
?np.fft # This gives us information on the conventions used in the return values of the functions.
In [ ]:
?np.fft.fft # This is the main DFT function we will use.
In [ ]:
?np.fft.fftfreq # This is a helper function to prepare a vector of frequencies.
In [ ]:
?np.arange # Points in an evenly spaced interval.
This approach is derived from a nice discussion on FFT found on the blog Glowing Python.
First we will divide up time into samplingInterval sized chunks between 0 and 1. This will aid in getting the x-axis scaled correctly so that frequency can be read directly off the DFT result. You can take samplingInterval in seconds putting samplingRate in Hz. Notice the approach here - we could have done this all in one line, but, by intelligently naming our variables and exposing the details of our thoughts the code is more readable:
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
samplingRate = 150.0
samplingInterval = 1.0/samplingRate
timeVector = np.arange(0, 1, samplingInterval)
# Print out the first few elements so you can see what is going on:
timeVector[0:10:]
Next we decide on the frequency of our signal and create a list to have a signal to work with.
In [ ]:
signalFrequency = 10.0;
ourSignal = np.sin(2*np.pi*signalFrequency*timeVector) + 0.5*np.sin(2*np.pi*(2*signalFrequency)*timeVector)
Plotting the input function for clarity:
In [ ]:
fig = plt.figure()
axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
axes.plot(timeVector, ourSignal, 'r')
axes.set_xlabel('Time')
axes.set_ylabel('Signal')
axes.set_title('Our Modest Signal');
Using numpy to compute the DFT:
In [ ]:
n = ourSignal.size
frequencies = np.fft.fftfreq(n, d=1.0/samplingRate)
spectrum = np.abs(np.fft.fft(ourSignal))
fig = plt.figure()
axes = fig.add_axes([0.1, 0.1, 0.8, 0.8])
axes.scatter(frequencies, spectrum, c='r', marker='s', alpha=0.4)
axes.set_xlabel('Frequency')
axes.set_ylabel('Amplitude')
axes.set_title('Our Amplitude Spectrum');
Transmission electron microscopy utilizes diffraction to determine crystal structures and develop contrast in images. In this section of the lecture we will simulate the diffraction pattern of an atomic structure. Using this diffraction pattern we will simulate using a diffraction aperture to reconstruct a phase contrast image.
In [ ]:
%matplotlib inline
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
from numpy.fft import *
def atomic_func(x,y):
param = 64.0
return (1+np.sin(4*(x+y)*2*np.pi/param))*(1+np.sin(2*(x-2*y)*2*np.pi/param))/4
def aperture(X, Y, xoffset, yoffset, size):
return (X-xoffset)**2+(Y-yoffset)**2 > size**2
We define two functions above:
atomic_func is used to provide an image function periodic in two dimensions from which the diffraction pattern will be constructed. This can be thought of as the density of electrons in a solid that is used to approximate a crystal structure.aperture returns a Boolean array that will be used to mask the diffraction pattern so that individual frequencies can be selected for image reconstruction. aperture will return True or False.
In [ ]:
x = np.arange(0.0,256.0,1.0)
y = np.arange(0.0,256.0,1.0)
X,Y = np.meshgrid(x, y)
Z = atomic_func(X,Y)
The Z array holds the atomic image function.
In [ ]:
P = np.zeros(Z.shape,dtype=complex)
K = np.zeros(Z.shape,dtype=complex)
K = fftshift(fft2(Z, norm='ortho'))
P = np.copy(K)
P[np.where(aperture(X, Y, 128, 128, 3) & aperture(X, Y, 150, 128, 3))] = 0
The P array holds the processed Fourier spectrum. The values of P are set to zero when they are outside the aperture. We use the K array to hold a opy of the image
In this cell we create two more numpy arrays (there are other ways to do this) that have the same shape as Z. The P array we use to hold the processed Fourier spectrum. The processing uses numpy's Boolean indexing to set values in P equal to zero when they are "outside" the aperture. When we get to the images below you'll see what is meant.
Because Python passes by reference we need to call for a copy of K so that we can modify one without changing the other.
From this processed spectrum we will create an image. The K array holds the whole Fourier spectrum.
In [ ]:
Im = fftshift(ifft2(P))
Above we reprocess P into the image Im.
In [ ]:
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(30,9))
axes[0].imshow(Z, origin='lower')
axes[1].imshow(abs(K),origin='lower', cmap=plt.get_cmap('pink'))
aperture1 = plt.Circle((128,128),3**2,color='r', fill = False)
aperture2 = plt.Circle((150,128),3**2,color='y', fill = False)
axes[1].add_artist(aperture1)
axes[1].add_artist(aperture2)
axes[2].imshow(abs(Im)**2, origin='lower')
plt.show()
Apply the DFT to an image of your choosing. Select the low frequency part of the DFT and regenerate the image (i.e. take the inverse FFT) from only these selected frequencies. Use a Boolean selection to zero out parts of the frequency spectrum before you convert back. To read an image in from disk, use the ndimage function from SciPy:
from scipy.ndimage import imread
img = imread('./images/pattern2.jpg', mode='L')
checking the data type of img will prove helpful.