Re: FFT of sinusoidal grating

Posted by Ben Tupper on Oct 02, 2008; 1:27pm
URL: http://imagej.273.s1.nabble.com/Competition-ImageJ-Video-Tutorials-tp3694884p3694889.html

Hi,

On the IDL side of things, you might find David Fanning's tutorial on  
frequency filtering useful.  Check it out at ...

http://dfanning.com/ip_tips/freqfiltering.html


Cheers,
Ben

On Oct 2, 2008, at 8:55 AM, Gabriele Umbriaco wrote:

> Thanks, I have understood my error.
> I am trying to confront the result of the FFT with that Power  
> Spectrum obtained to the focus of lens with parrallel beam and  
> amplitude mask.
> Gabriele
>
>
> Stefan Heim ha scritto:
>> Hi Gabriele,
>>
>> On Thu, 2008-10-02 at 11:45 +0200, Gabriele Umbriaco wrote:
>>
>>> If I compute the FFT with ImageJ the result is not the same of IDL.
>>> I do not understand because, it is wrong ?
>>> You can help me?
>>>
>>
>> General advise: it takes a lot of caution and double-checking to  
>> compare
>> FFT results from two distinct FFT implementations.
>>
>>
>>>> IDL code:
>>>>
>>>> ; IDL Version 5.2 (Win32 x86)
>>>> ; Journal File for rockford@C3B0E6
>>>> ; Working directory: D:\RSI\IDL52
>>>> ; Date: Sat Jan 17 18:09:01 2004
>>>> cd, 'C:\immagini'
>>>> read_jpeg,'sin_grid_30.jpg',ima
>>>> image=fltarr(512,512)
>>>> image(*,*)=ima(1:512,1:512)
>>>>
>>
>> This looks suspicious. sin_grid_30.jpg is a periodic 800x800 pixel
>> image, and you pick a 512x512 subimage, (accidentally?) offset by 1x1
>> pixel (array indices in IDL are zero-based!). You don't have a  
>> periodic
>> image here anymore, which will skew the resulting power spectrum  
>> unless
>> you pad or taper. Using jpeg compression btw is also suboptimal  
>> for your
>> sinusoidal grating, as it will introduce artifacts in your pristine
>> frequency spectrum.
>>
>>
>>>> dimx=512
>>>> dimy=512
>>>> window, 1, xs=dimx, ys=dimy, title='original'
>>>> tvscl,image
>>>>
>>>> ;Power Spectrum F(u)*F^(u)
>>>> window, 2, xs=dimx, ys=dimy, title='Filtered Power Spectrum'
>>>> pspect= shift((fft(image, -1)*fft(conj(image),1)),dimx/2,dimy/2)
>>>>
>>
>> Inefficient (and strictly speaking wrong by a factor of 2) way to
>> calculate the power spectrum. Try
>>
>> power = SHIFT(ABS(FFT(image, -1)), dimx/2, dimy/2)
>>
>> instead. For any real world images, you'll also want to have a
>> logarithmic power spectrum for displaying purposes, i.e.
>>
>> spectrum = FFT(image, -1)
>> zero_dB = ABS(spectrum(0,0))
>> power_dB = SHIFT(20. * ALOG10(ABS(spectrum) / zero_dB), dimx/2,  
>> dimy/2)
>>
>> for a power spectrum in decibel with DC (mean gray value in the  
>> image)
>> as reference level.
>>
>> If you check the "Raw Power Spectrum" and uncheck the "FFT Window"
>> checkboxes under Process->FFT->FFT Options in ImageJ, the resulting
>> power spectrum will _look_ the same as your IDL calculated one as a
>> result of the DC component being so dominant that it completely  
>> messes
>> up linear scaling.
>>
>> Differences in the frequency noise are a mere artifact of the two
>> distinct implementations of the FFT algorithm.
>>
>> HTH,
>> -Stefan
>>
>>
>>
>
> --
> =================================================
> Gabriele Umbriaco
> Dipartimento di Astronomia, Universita` di Padova
> Laboratorio d'ottica Padova/Asiago
> vicolo dell'Osservatorio 3 I-35122 Padova, Italy
> office: ([phone] fax: [phone]
> laboratory PD: ([phone]
> laboratory AS: ([phone]
> email: [hidden email]