Loading...
 

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
from jlu.util import radialProfile

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

    img = pyfits.getdata(imageFile)

    F = fftpack.fftshift( fftpack.fft2(img) )
    psd2D = np.abs(F)**2
    psd1D = radialProfile.azimuthalAverage(psd2D)
    
    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

Page last modified on Tuesday 13 of September, 2011 16:31:46 EDT

Attached files

ID Name Comment Uploaded Size Downloads
12 image_fft_examp_psd1d.png jlu Wed 03 of Mar, 2010 18:00 EST 30.27 Kb 2297
11 image_fft_examp_psd2d.png jlu Wed 03 of Mar, 2010 18:00 EST 355.32 Kb 2241
10 image_fft_examp_image.png jlu Wed 03 of Mar, 2010 18:00 EST 175.35 Kb 2185
9 image_fft.py jlu Wed 03 of Mar, 2010 18:00 EST 961 b 2701