Re: FFT of sinusoidal grating

Posted by Gabriele Umbriaco on
URL: http://imagej.273.s1.nabble.com/Competition-ImageJ-Video-Tutorials-tp3694884p3694888.html

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]