PHY324H1S
Computational Assignment: Fast Fourier Transform applications
in Physics using Python
Due Date: March 20, 2014. Submit it to your group leader.
Introduction
Fourier transform methods (spectral methods) are widely used in research as efficient
data manipulating tools: signal processing, digital frequency filters, data analysis, etc.
An accessible introduction to Fourier analysis and its applications can be found at:
http://www.upscale.utoronto.ca/PVB/Harrison/FourierTransform.pdf
Check it out before writing the Python program.
Fourier analysis is a method expressing a function as a sum of periodic components and
recovering the output “signal” from those components. The Discrete Fourier Transform
separates the input (signal in the time domain) into components that contribute at discrete
frequencies. The output is called a spectrum or transform and exists in the frequency
domain.
A physical process can be represented by a function of time h(t) or by a function of
frequency H(f) where -∞ < f < ∞. h(t) and H(f) are two different representations of the
same function. The Fourier transform enables the switch back and forth between the two
representations through the following equations (the Fourier transform equations):

H( f ) 
(1)
 2ift
h
(
t
)
e
dt



h(t ) 
2ift
H
(
f
)
e
df

(1’)

(1) is the Fourier Transform of time-dependent function h(t) and (1’) is the Inverse
Fourier Transform.
The Discrete Fourier Transform is the numerical approximation to the Fourier transform.
DFT is defined as:
N 1
X k   xn e
2 i
kn
N
(2)
n 0
Where:
x = sequence of N numbers (data points)
Xk = k-th point in the Fourier-transformed data
The DFT is invertible, with its inverse given by:
1
xn 
N
N 1
X
k 0
k
e
2 i
kn
N
(3)
1
The k-th element in the data sequence can be recovered from the Fourier sequence by
adding up a series of complex waves Ak e it , where Ak  X k N 1 , and   2  k N  1 , so
each point Xk in the Fourier sequence corresponds to the amplitude and phase shift of a
complex wave, whose frequency depends on its position, k, in the sequence.
DFT is very important in numerical analysis in part because of the very fast algorithm for
computing it, called the Fast Fourier Transform (FFT). The algorithm was developed by
Cooley and Tukey [3].
The FFT algorithm is equivalent to equations (2) and (3), but is more computationally
efficient than the definition.
FFT and Python
The numpy module has a built-in FFT package called fft. Complete documentation can
be found at:
http://docs.scipy.org/doc/numpy/reference/routines.fft.html
Function fft simply outputs an array that represents a Fourier transform of given data
(which is also in a form of an array).
If A = fft(a, n), then A[0] contains the zero-frequency term (the mean of the
signal), which is always purely real for real inputs.
A[1:n/2] contains the positive-frequency terms, and A[n/2+1:] contains the
negative-frequency terms, in order of decreasingly negative frequency.
For an even number of input points, A[n/2] represents both positive and negative
frequencies, and is also purely real for real input.
For an odd number of input points, A[(n-1)/2] contains the largest positive
frequency, while A[(n+1)/2] contains the largest negative frequency.
Syntax: fft.fft(a, n=None, axis=-1)
This function computes the one-dimensional n-point discrete Fourier Transform (DFT)
with the efficient Fast Fourier Transform (FFT) algorithm.
Parameters:
a array_like. Input array can be complex.
n : int, optional. Length of the transformed axis of the output.
axis : int, optional. Axis over which to compute the FFT.
Returns:
out : complex ndarray
Error raised:
IndexError :if axes is larger than the last axis of a.
Function fftfreq outputs the frequency bins used to perform FFT.
2
fft.fftfreq(A) returns an array giving the frequencies of corresponding elements in
the output.
fft.fftshift(A) shifts transforms and their frequencies to put the zero-frequency
components in the middle, and fft.ifftshift(A) undoes that shift.
When the input a is a time-domain signal and A = fft(a), abs(A) is its amplitude
spectrum and abs(A)**2 is its power spectrum. The phase spectrum is obtained by
angle(A).
For example, let f be an array that represents data.
G = fft(f)
freq = fftfreq(len(f), dt)
G is the Fourier transform of f. Then, freq[i] is the frequency of the ith component
of G. dt is the increment of time at which f is defined as: f[i] = f(dt*i). If not
specified, dt is 1 by default.
In Table 1 you will find the numpy commands related to FFT:
Function
fft
ifft
fft2
ifft2
fftshift
ifftshift
fftfreq
Description
Performs FFT
Performs Inverse FFT
Perform FFT on 2D data
Perform inverse FFT on 2D data
Shifts the order of the output from the FFT
Inverse of fftshift
Generate frequencies corresponding to each point in
the Fourier-transformed data.
Table 1 FFT-related commands in numpy
Frequency analysis exercises
 Exercise 1
Import libraries (pylab, scipy, numpy, numpy.fft) Generate data (sinusoidal
functions). Use fft to analyze:
a) A superposition of three waves with constant, but different frequencies
b) A wave whose frequency depends on time
Plot the graphs, sample frequencies, amplitude and power spectra with appropriate axes.
 Exercise 2
Yearly sunspot measurements are provided by the US National Geophysical Data
Centre/NOAA Satellite and Information Service. Data will be provided as a text file. Use
the data with fft to find the dominant period of the sunspot observations.
3
Bandpass filtering
FFT is not only useful for finding frequency components in waves. It can also be used to
remove noise from experimental data.
Often, the “white noise” (Gaussian-distributed random noise) contaminates all frequency
bands, uniformly. If we know the frequency range of the signal, we can remove the noise
spikes away from this frequency range, and then perform an inverse FFT to get a cleaned
data set. For example, with a frequency of the waves of 2 Hz, we will simply multiply the
noisy FFT spectrum with a Gaussian function:
 (| |  2 ) 2
F ( )  e
2 ( 0.5 2 )
(4)
Where: F is the Gaussian filter, ν = frequency (Hz)
The Gaussian function has a maximum height of 1, is centered at 2Hz, with a standard
deviation of 0.5Hz.
Bandpass filtering exercise
 Exercise 3
Use the 3 waves from Exercise 1a), but change them to have the same frequency, phase
and amplitude. Contaminate each of them with successively increasing amounts of
random, Gaussian-distributed noise.
1) Perform an FFT on the superposition of the three noise-contaminated waves.
Analyze and plot the output.
2) Filter the signal with a Gaussian function, plot the “clean” wave, and analyze the
result. Is the resultant wave 100% clean? Explain.
3) [CANCELLED] Use the function from Exercise 1b). Modify it to have a linear
dependency of frequency on time. Contaminate the function with random noise.
Perform the FFT. On the FFT, apply a Gaussian filter of your choice and clean the
noise. The center of the Gaussian filter should be at the midpoint of the frequency
range.
References:
1. G.B. Arfken, H.J. Weber: Mathematical Methods for Physicists, 4th ed., AP1995
2. W.H. Press, S.A.Teukolsky, W.T. Vetterling, B.P. Flannery: Numerical Recipes:
The Art of Scientific Computing, 3rd ed. Cambridge Univ. Press 2007
3. James W. Cooley, John W. Tukey: “An algorithm for the machine calculation of
complex Fourier series,” Math. Comput. 19: 297-301, 1965.
4. Numpy Reference: http://docs.scipy.org/doc/numpy/reference/index.html
(1.5.dev8088, February 05, 2010)
This exercise was written by: Ruxandra Serbanescu in 2010. Students: Frank Shu Zhao, Christopher
Stevens, Eve Lee, Dina Mistry, and Siming Zhang wrote and tested preliminary versions of this exercise
guide in 2009.
4

P型フーリエ記述子を用いた前髪形状からのマンガ登場人物識別

Appendix A Curvilinear coordinates