PRO Freqfilt

   ; Read the convec.dat data file.

file = Filepath('convec.dat', SubDir=['examples', 'data'])
OpenR, lun, file, /Get_Lun
convec = BytArr(248,248)
ReadU, lun, convec
Free_Lun, lun

   ; Open a window and load a color table.

WINDOW, 0, XSIZE=248*3, YSIZE=248, Title='Low Pass Frequency Filtering'
LOADCT, 13

   ; Display the original image.

TV, convec, 0

   ; Transform the image to the frequency domain.

freqDomainImage = FFT(convec, -1)

   ; Calculate the power spectrum and display it.

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

   ; Create and apply the low-pass filter.

filter = 1.0 / (1.0d + DIST(248)/15.0)^2
filtered = FFT(freqDomainImage * filter, 1)

   ; Display the filtered image.

TV, filtered, 2

END