### Contributed by Jessica R. Lu

Below is an example of calculating a 1D and 2D power spectrum from an image. The main part of the code is presented below with some example figures from one of my own astronomical images. This code uses radialProfile.py.  An alternative implementation (based on J. Lu's code below) is available at http://code.google.com/p/agpy/source/browse/trunk/AG_fft_tools.

```import math
from scipy import fftpack
import pyfits
import numpy as np
import pylab as py

def run():
imageFile = 'myimage.fits'

img = pyfits.getdata(imageFile)

F = fftpack.fftshift( fftpack.fft2(img) )
psd2D = np.abs(F)**2

py.figure(1)
py.clf()
py.imshow(np.log10(img), cmap=py.cm.Greys,
vmin=math.log10(20), vmax=math.log10(3500))
py.title('Image')
py.savefig('image_fft_examp_image.png')

py.figure(2)
py.clf()
py.imshow(np.log10(psd2D), cmap=py.cm.jet)
py.title('2D Power Spectrum')
py.savefig('image_fft_examp_psd2d.png')

py.figure(3)
py.clf()
py.semilogy(psd1D)
py.title('1D Power Spectrum')
py.xlabel('Spatial Frequency')
py.ylabel('Power Spectrum')
py.savefig('image_fft_examp_psd1d.png')

py.show()```

The resulting figures are:

And the code is here: image_fft.py