PRO PowerSpectrum, image

   ; Open an image data file if image not passed in.

IF N_PARAMS() EQ 0 THEN BEGIN
   file = Filepath('convec.dat', SubDir=['examples','data'])
   OpenR, lun, file, /Get_Lun
   image = BytArr(248,248)
   ReadU, lun, image
   Free_Lun, lun
ENDIF

   ; Get image dimensions.

s = SIZE(image)
xsize = s(1)
ysize = s(2)

   ; Open a window and load a color table.

WINDOW, Title='Power Spectrum Example', $
   XSIZE=xsize*2, YSIZE=ysize*3
LOADCT, 0

   ; Display the original image.

TV, image, 0

   ; Transform the image into the frequency domain.

freqDomainImage = FFT(image, -1)

   ; Calculate and display the power spectrum.

power = SHIFT(ALOG(ABS(freqDomainImage)), xsize/2, ysize/2)
power = power - Min(power)
power = power * (255.0/Max(power))
TV, power, 1

END