Fourier Transform

Worksheet 1 focuses on using Python tasks to calculate the Fourier Transform of a few window functions. In this notebook, I will illustrate how you can generate a window function, and how you can calculate the associated Fourier Transform.

The Fourier Transform can be examined in many ways. The simplest interpretation is that the FT is a harmonic decomposition of a signal. You can calculate the FT analytically, if you’re lucky. If you’re even luckier, a finite Fourier Series will give you an accurate representation of your signal – provided that you have a simple signal. There is often an implicit assumption when calculating the analytical representation – that your input function goes from -infinity to +infinity.

In practice, astronomical signals are not perfect, simple sine-waves. Neither are they perfectly sampled, or sampled for infinite time. When you have a regularly sampled series, e.g., a time series, you can compute the Fourier Transform using the Fast Fourier Transform (FFT) Algorithm. This is what we’re doing in this worksheet – we are computing the FFT for a series of functions, and we are studying how similar or different the FT is to what we expect from the illustrations on the worksheet.

When in doubt, Google

There’s a reason why WS1 is short on details – there is a huge amount of resources on the internet that focuses on computing Fourier Transforms using python. In fact, if you Google Python Fourier Transform, the first hit will probably be the relevant SCIPY page.

Please note that there are lots of resources on how to generate window functions, such as the boxcar function. You can use a similar method to the one I’ve done below, or you can simply use a python function to define this. I leave it to you to figure it out.

Getting started

First, I’m going to import the important libraries, as suggested in the aforementioned URL. Note that %matplotlib inline command allows us to render plots within this notebook.

Also note how I have defined pylab and numpy in the Python namespace. This makes it easy to reference plotting and analytical functions from those libraries.

import numpy.fft as fft
import pylab as pl
import numpy as np
%matplotlib inline

The Shah Function

The Shah function (Comb) is a series of delta functions with some spacing dx.

The image below is from your WS1, and is one of the examples found in the ERA appendix on Fourier Transforms:

In the next cell I define the Shah function, by hand. Note how I use the np.mod function. This returns the remainder of division, so I’ve inserted an if statement that appends a 1 to the array if the index i is divisible by dx, and 0 if it isn’t. There are certainly many more methods to do this, but let’s do it the hard way for now.

I’ve also included comments on how I’ve constructed the Shah function and the associated plot.

# This is the step in the x direction.
dx = 10 
# Nx defines how many samples. The for-loop generates twice the number
# of samples defined in Nx.
Nx = 100

# Define empty arrays for x and y
x = []
y = []

# The for-loop runs from -Nx to +Nx in steps of 1.
for i in range(-Nx,Nx+1):
    x.append(i)
    # The if statement appends a 1 if (i,5)=0, and 0 if (i,5)!=0.
    if np.mod(np.abs(i), dx)<1:
        y.append(1.)
    else:
        y.append(0.)

# Python likes to do math on arrays, so I'm going to covert x and y.
x = np.array(x)
# Now y contains the Shah function, which I'll plot shortly.
y = np.array(y)
# Note that I'm using the step plotting function.
pl.step(x, y)
# Set the y-limits on the plot.
pl.ylim(0,2)

Calculating the Fourier Transform

Now that we have generated the Shah function, we can use the fft library defined above.

# I want the FFT to be smooth and well sampled, so I'm going to define 
# a variable that sets the number of points for the FFT.
nftpoints = 100

# Now, I can calculate the FFT of y.
# using norm='ortho' will normalize the FFT, and n=nftpoints will make 
# the FFT return the FT sampled with nftpoints.
F = fft.fft(y, norm='ortho', n=nftpoints)


# There is a method to calculate the associated frequency
# in the Fourier domain. 
# Now I'm going to include a fairly complicated command. 
# fft.fftfreq computes the frequency spacing, and fft.fftshift 
# shifts/sorts the frequency axis so that it plots nicely.
# The important command is fftfreq. 
freq = fft.fftshift(fft.fftfreq(nftpoints))

Plotting the FFT

Now that I’ve generated the window function and the associated FT, I can finally plot them in a way similar to the NRAO plot.

Here the subplot, step and ylim commands are important.

# First, plot the Shah function. 
pl.subplot(121)
# Plot this in a black line, this is what 'k-' is for.
pl.step(x, y, 'k-')
# Set the y-limit so that it looks pretty.
pl.ylim(0,2)
pl.title('Shah Function')


# Now, plot the FT.
pl.subplot(122)
# Note that I am plotting the absolute value of F, since F is an 
# array of complex numbers.
pl.step(freq, np.absolute(F), 'k-')
# Set the y-limit.
pl.ylim(0,2)
pl.title('FT of Shah Function')

<matplotlib.text.Text at 0x7ff2673bfc50>

We need to note something important here. Recall that I set dx when I defined the Shah function above. Now, from the NRAO plot, I know that space between the spikes in the Fourier domain is given by the following: ds = 1/dx. You will note that this is the case in the plot above.