Hey! I'm a beginner at IDL and have problem with FFT. I'm trying to
perform 2d-FFT but the code doesn't work properly even on test images.
So I create an image, make the Fourier transform, then the inverse
Fourier transform and finally I expect to get the initial image. But
the resulting image is the initial one, reversed with respect to the
center.
My code:
nx=256L
x1=findgen(nx)-nx/2.0+1.0
x2=findgen(nx)-nx/2.0+1.0
ytest=fltarr(nx,nx)
for i=0l,nx-1 do begin
for j=0l,nx-1 do begin
if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 and x2[j] ge
0.0) then begin
ytest[i,j]=1.0
endif
endfor
endfor
; So the initial image is a squre in the right upper corner
FFTtest=FFT(ytest)
sh_FFTtest=SHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
inv_test=FFT(FFTtest,-1)
;The result is the square in the LEFT LOWER corner.
I would be very grateful for comments/advices/solutions
Thank for help,
Natalya
|
|
0
|
|
|
|
Reply
|
natalya.lyskova (1)
|
12/5/2010 5:48:54 PM |
|
Natalya Lyskova writes:
> Hey! I'm a beginner at IDL and have problem with FFT. I'm trying to
> perform 2d-FFT but the code doesn't work properly even on test images.
> So I create an image, make the Fourier transform, then the inverse
> Fourier transform and finally I expect to get the initial image. But
> the resulting image is the initial one, reversed with respect to the
> center.
>
> My code:
> nx=256L
> x1=findgen(nx)-nx/2.0+1.0
> x2=findgen(nx)-nx/2.0+1.0
>
> ytest=fltarr(nx,nx)
> for i=0l,nx-1 do begin
> for j=0l,nx-1 do begin
> if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 and x2[j] ge
> 0.0) then begin
> ytest[i,j]=1.0
> endif
> endfor
> endfor
>
> ; So the initial image is a squre in the right upper corner
>
> FFTtest=FFT(ytest)
> sh_FFTtest=SHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
> inv_test=FFT(FFTtest,-1)
>
> ;The result is the square in the LEFT LOWER corner.
>
> I would be very grateful for comments/advices/solutions
I think you need to read the on-line help for the FFT function. :-)
Your code should look like this:
FFTtest = FFT(ytest, -1)
inv_test = Real_Part(FFT(FFTtest, 1))
Now ytest and inv_test are essentially the same.
Cheers,
David
--
David Fanning, Ph.D.
Fanning Software Consulting, Inc.
Coyote's Guide to IDL Programming: http://www.dfanning.com/
Sepore ma de ni thui. ("Perhaps thou speakest truth.")
|
|
0
|
|
|
|
Reply
|
David
|
12/5/2010 7:06:30 PM
|
|
On Dec 5, 2:06=A0pm, David Fanning <n...@dfanning.com> wrote:
> Natalya Lyskova writes:
> > Hey! I'm a beginner at IDL and have problem with FFT. I'm trying to
> > perform 2d-FFT but the code doesn't work properly even on test images.
> > So I create an image, make the Fourier transform, then the inverse
> > Fourier transform and finally I expect to get the initial image. But
> > the resulting image is the initial one, reversed with respect to the
> > center.
>
> > My code:
> > nx=3D256L
> > x1=3Dfindgen(nx)-nx/2.0+1.0
> > x2=3Dfindgen(nx)-nx/2.0+1.0
>
> > ytest=3Dfltarr(nx,nx)
> > for i=3D0l,nx-1 do begin
> > for j=3D0l,nx-1 do begin
> > =A0if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 =A0and x2[j] ge
> > 0.0) then begin
> > =A0 =A0 ytest[i,j]=3D1.0
> > =A0endif
> > endfor
> > endfor
>
> > ; So the initial image is a squre in the right upper corner
>
> > FFTtest=3DFFT(ytest)
> > sh_FFTtest=3DSHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
> > inv_test=3DFFT(FFTtest,-1)
>
> > ;The result is the square in the LEFT LOWER corner.
>
> > I would be very grateful for comments/advices/solutions
>
> I think you need to read the on-line help for the FFT function. :-)
>
> Your code should look like this:
>
> =A0 =A0FFTtest =3D FFT(ytest, -1)
> =A0 =A0inv_test =3D Real_Part(FFT(FFTtest, 1))
>
> Now ytest and inv_test are essentially the same.
>
> Cheers,
>
> David
>
> --
> David Fanning, Ph.D.
> Fanning Software Consulting, Inc.
> Coyote's Guide to IDL Programming:http://www.dfanning.com/
> Sepore ma de ni thui. ("Perhaps thou speakest truth.")
Hello,
Im working with the FFT also but not directly in IDL, I work in ENVI.
I wonder how ENVI or IDL can perform a FFT on a rectangular image?
Theorically the image must be a square of dimension of coefficient of
2. (512x512, 1024 x1024 etc... ) So what is the process that make
possible to do a FFT on a rectangular image? Other Image Processing
software like PCI Geomatica cant do that.
Thank you,
Max
|
|
0
|
|
|
|
Reply
|
burton449
|
12/5/2010 8:05:30 PM
|
|
On Dec 5, 12:48=A0pm, Natalya Lyskova <natalya.lysk...@gmail.com> wrote:
> Hey! I'm a beginner at IDL and have problem with FFT. I'm trying to
> perform 2d-FFT but the code doesn't work properly even on test images.
> So I create an image, make the Fourier transform, then the inverse
> Fourier transform and finally I expect to get the initial image. But
> the resulting image is the initial one, reversed with respect to the
> center.
>
> My code:
> nx=3D256L
> x1=3Dfindgen(nx)-nx/2.0+1.0
> x2=3Dfindgen(nx)-nx/2.0+1.0
>
> ytest=3Dfltarr(nx,nx)
> for i=3D0l,nx-1 do begin
> for j=3D0l,nx-1 do begin
> =A0if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 =A0and x2[j] ge
> 0.0) then begin
> =A0 =A0 ytest[i,j]=3D1.0
> =A0endif
> endfor
> endfor
>
> ; So the initial image is a squre in the right upper corner
>
> FFTtest=3DFFT(ytest)
> sh_FFTtest=3DSHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
> inv_test=3DFFT(FFTtest,-1)
>
> ;The result is the square in the LEFT LOWER corner.
>
> I would be very grateful for comments/advices/solutions
>
> Thank for help,
> Natalya
As a side note, a better (both faster and more readable) way of
creating the original image instead of the double FOR loop is:
if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 and x2[j] ge
0.0) then begin
ytest[i,j]=3D1.0
upperright =3D where(x1 le 20.0 and x1 ge 0.0 and x2 le 20.0 and x2 ge
0.0, nupperright)
if nupperright gt 0 then ytest[upperright]=3D1.0
-Jeremy.
|
|
0
|
|
|
|
Reply
|
Jeremy
|
12/6/2010 2:16:07 PM
|
|
On 6 d=E9c, 15:16, Jeremy Bailin <astroco...@gmail.com> wrote:
> On Dec 5, 12:48=A0pm, Natalya Lyskova <natalya.lysk...@gmail.com> wrote:
>
>
>
>
>
> > Hey! I'm a beginner at IDL and have problem with FFT. I'm trying to
> > perform 2d-FFT but the code doesn't work properly even on test images.
> > So I create an image, make the Fourier transform, then the inverse
> > Fourier transform and finally I expect to get the initial image. But
> > the resulting image is the initial one, reversed with respect to the
> > center.
>
> > My code:
> > nx=3D256L
> > x1=3Dfindgen(nx)-nx/2.0+1.0
> > x2=3Dfindgen(nx)-nx/2.0+1.0
>
> > ytest=3Dfltarr(nx,nx)
> > for i=3D0l,nx-1 do begin
> > for j=3D0l,nx-1 do begin
> > =A0if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 =A0and x2[j] ge
> > 0.0) then begin
> > =A0 =A0 ytest[i,j]=3D1.0
> > =A0endif
> > endfor
> > endfor
>
> > ; So the initial image is a squre in the right upper corner
>
> > FFTtest=3DFFT(ytest)
> > sh_FFTtest=3DSHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
> > inv_test=3DFFT(FFTtest,-1)
>
> > ;The result is the square in the LEFT LOWER corner.
>
> > I would be very grateful for comments/advices/solutions
>
> > Thank for help,
> > Natalya
>
> As a side note, a better (both faster and more readable) way of
> creating the original image instead of the double FOR loop is:
>
> =A0if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 =A0and x2[j] ge
> 0.0) then begin
> =A0 =A0 ytest[i,j]=3D1.0
>
> upperright =3D where(x1 le 20.0 and x1 ge 0.0 and x2 le 20.0 and x2 ge
> 0.0, nupperright)
> if nupperright gt 0 then ytest[upperright]=3D1.0
>
> -Jeremy.- Masquer le texte des messages pr=E9c=E9dents -
>
Even simpler since IDL 8.0:
ytest[where(x1 le 20.0 and x1 ge 0.0 and x2 le 20.0 and x2 ge 0.0, /
NULL)]=3D1.0
alx.
|
|
0
|
|
|
|
Reply
|
alx
|
12/6/2010 3:10:24 PM
|
|
On Dec 5, 3:05=A0pm, burton449 <burton...@gmail.com> wrote:
> On Dec 5, 2:06=A0pm, David Fanning <n...@dfanning.com> wrote:
>
>
>
> > Natalya Lyskova writes:
> > > Hey! I'm a beginner at IDL and have problem with FFT. I'm trying to
> > > perform 2d-FFT but the code doesn't work properly even on test images=
..
> > > So I create an image, make the Fourier transform, then the inverse
> > > Fourier transform and finally I expect to get the initial image. But
> > > the resulting image is the initial one, reversed with respect to the
> > > center.
>
> > > My code:
> > > nx=3D256L
> > > x1=3Dfindgen(nx)-nx/2.0+1.0
> > > x2=3Dfindgen(nx)-nx/2.0+1.0
>
> > > ytest=3Dfltarr(nx,nx)
> > > for i=3D0l,nx-1 do begin
> > > for j=3D0l,nx-1 do begin
> > > =A0if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 =A0and x2[j] =
ge
> > > 0.0) then begin
> > > =A0 =A0 ytest[i,j]=3D1.0
> > > =A0endif
> > > endfor
> > > endfor
>
> > > ; So the initial image is a squre in the right upper corner
>
> > > FFTtest=3DFFT(ytest)
> > > sh_FFTtest=3DSHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
> > > inv_test=3DFFT(FFTtest,-1)
>
> > > ;The result is the square in the LEFT LOWER corner.
>
> > > I would be very grateful for comments/advices/solutions
>
> > I think you need to read the on-line help for the FFT function. :-)
>
> > Your code should look like this:
>
> > =A0 =A0FFTtest =3D FFT(ytest, -1)
> > =A0 =A0inv_test =3D Real_Part(FFT(FFTtest, 1))
>
> > Now ytest and inv_test are essentially the same.
>
> > Cheers,
>
> > David
>
> > --
> > David Fanning, Ph.D.
> > Fanning Software Consulting, Inc.
> > Coyote's Guide to IDL Programming:http://www.dfanning.com/
> > Sepore ma de ni thui. ("Perhaps thou speakest truth.")
>
> Hello,
> Im working with the FFT also but not directly in IDL, I work in ENVI.
> I wonder how ENVI or IDL can perform a FFT on a rectangular image?
> Theorically the image must be a square of dimension of coefficient of
> 2. (512x512, 1024 x1024 etc... ) So what is the process that make
> possible to do a FFT on a rectangular image? Other Image Processing
> software like PCI Geomatica cant do that.
Taking a DFT (Discrete Fourier Transformation) of an array is possible
for any array size. There is an algorithm called FFT (Fast Fourier
Transformation)
that happens to be very efficient if the size is in the form 2^N for
some N.
However, modern incarnations of FFT can deal with other sizes too,
albeit less
efficiently (the smaller the factors in the prime decomposition of the
size,
the better).
Please don't let what a particular software does or fail to do be your
guide to what is possible or not (from a mathematical standpoint). If
you
want to learn a bit more about the FFT, read for instance the chapter
about
it on the numerical recipes book.
Ciao,
Paolo
>
> Thank you,
> Max
|
|
0
|
|
|
|
Reply
|
Paolo
|
12/6/2010 3:28:17 PM
|
|
On Dec 6, 4:10=A0pm, alx <lecacheux.al...@wanadoo.fr> wrote:
> Even simpler since IDL 8.0:
>
> ytest[where(x1 le 20.0 and x1 ge 0.0 and x2 le 20.0 and x2 ge 0.0, /
> NULL)]=3D1.0
Or, even simpler, and without a need for v8.0:
ytest =3D (x1 le 20.0 and x1 ge 0.0 and x2 le 20.0 and x2 ge 0.0)
Cheers
Timm
|
|
0
|
|
|
|
Reply
|
Timm
|
12/6/2010 4:26:36 PM
|
|
On Dec 6, 10:28=A0am, Paolo <pgri...@gmail.com> wrote:
> On Dec 5, 3:05=A0pm, burton449 <burton...@gmail.com> wrote:
>
>
>
> > On Dec 5, 2:06=A0pm, David Fanning <n...@dfanning.com> wrote:
>
> > > Natalya Lyskova writes:
> > > > Hey! I'm a beginner at IDL and have problem with FFT. I'm trying to
> > > > perform 2d-FFT but the code doesn't work properly even on test imag=
es.
> > > > So I create an image, make the Fourier transform, then the inverse
> > > > Fourier transform and finally I expect to get the initial image. Bu=
t
> > > > the resulting image is the initial one, reversed with respect to th=
e
> > > > center.
>
> > > > My code:
> > > > nx=3D256L
> > > > x1=3Dfindgen(nx)-nx/2.0+1.0
> > > > x2=3Dfindgen(nx)-nx/2.0+1.0
>
> > > > ytest=3Dfltarr(nx,nx)
> > > > for i=3D0l,nx-1 do begin
> > > > for j=3D0l,nx-1 do begin
> > > > =A0if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 =A0and x2[j=
] ge
> > > > 0.0) then begin
> > > > =A0 =A0 ytest[i,j]=3D1.0
> > > > =A0endif
> > > > endfor
> > > > endfor
>
> > > > ; So the initial image is a squre in the right upper corner
>
> > > > FFTtest=3DFFT(ytest)
> > > > sh_FFTtest=3DSHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
> > > > inv_test=3DFFT(FFTtest,-1)
>
> > > > ;The result is the square in the LEFT LOWER corner.
>
> > > > I would be very grateful for comments/advices/solutions
>
> > > I think you need to read the on-line help for the FFT function. :-)
>
> > > Your code should look like this:
>
> > > =A0 =A0FFTtest =3D FFT(ytest, -1)
> > > =A0 =A0inv_test =3D Real_Part(FFT(FFTtest, 1))
>
> > > Now ytest and inv_test are essentially the same.
>
> > > Cheers,
>
> > > David
>
> > > --
> > > David Fanning, Ph.D.
> > > Fanning Software Consulting, Inc.
> > > Coyote's Guide to IDL Programming:http://www.dfanning.com/
> > > Sepore ma de ni thui. ("Perhaps thou speakest truth.")
>
> > Hello,
> > Im working with the FFT also but not directly in IDL, I work in ENVI.
> > I wonder how ENVI or IDL can perform a FFT on a rectangular image?
> > Theorically the image must be a square of dimension of coefficient of
> > 2. (512x512, 1024 x1024 etc... ) So what is the process that make
> > possible to do a FFT on a rectangular image? Other Image Processing
> > software like PCI Geomatica cant do that.
>
> Taking a DFT (Discrete Fourier Transformation) of an array is possible
> for any array size. There is an algorithm called FFT (Fast Fourier
> Transformation)
> that happens to be very efficient if the size is in the form 2^N for
> some N.
> However, modern incarnations of FFT can deal with other sizes too,
> albeit less
> efficiently (the smaller the factors in the prime decomposition of the
> size,
> the better).
>
> Please don't let what a particular software does or fail to do be your
> guide to what is possible or not (from a mathematical standpoint). If
> you
> want to learn a bit more about the FFT, read for instance the chapter
> about
> it on the numerical recipes book.
>
> Ciao,
> Paolo
>
>
>
> > Thank you,
> > Max
>
>
Hi Paolo,
Thank you for your comment. As a student in Remote Sensing, I have a
lot of basic things to understand. The image I would like to filter in
the frequential domain (using a Butterworth filter) is a side scan
sonar image mosaic of 9166 x 4093 pixels. Would you recommend a FFT
and if yes would you use a special algorithm?
greetings,
Max
|
|
0
|
|
|
|
Reply
|
burton449
|
12/7/2010 2:33:39 AM
|
|
On Dec 6, 9:33=A0pm, burton449 <burton...@gmail.com> wrote:
> On Dec 6, 10:28=A0am, Paolo <pgri...@gmail.com> wrote:
>
>
>
> > On Dec 5, 3:05=A0pm, burton449 <burton...@gmail.com> wrote:
>
> > > On Dec 5, 2:06=A0pm, David Fanning <n...@dfanning.com> wrote:
>
> > > > Natalya Lyskova writes:
> > > > > Hey! I'm a beginner at IDL and have problem with FFT. I'm trying =
to
> > > > > perform 2d-FFT but the code doesn't work properly even on test im=
ages.
> > > > > So I create an image, make the Fourier transform, then the invers=
e
> > > > > Fourier transform and finally I expect to get the initial image. =
But
> > > > > the resulting image is the initial one, reversed with respect to =
the
> > > > > center.
>
> > > > > My code:
> > > > > nx=3D256L
> > > > > x1=3Dfindgen(nx)-nx/2.0+1.0
> > > > > x2=3Dfindgen(nx)-nx/2.0+1.0
>
> > > > > ytest=3Dfltarr(nx,nx)
> > > > > for i=3D0l,nx-1 do begin
> > > > > for j=3D0l,nx-1 do begin
> > > > > =A0if (x1[i] le 20.0 and x1[i] ge 0.0 and x2[j] le 20.0 =A0and x2=
[j] ge
> > > > > 0.0) then begin
> > > > > =A0 =A0 ytest[i,j]=3D1.0
> > > > > =A0endif
> > > > > endfor
> > > > > endfor
>
> > > > > ; So the initial image is a squre in the right upper corner
>
> > > > > FFTtest=3DFFT(ytest)
> > > > > sh_FFTtest=3DSHIFT(FFTtest,nx/2.0-1.0,nx/2.0-1.0)
> > > > > inv_test=3DFFT(FFTtest,-1)
>
> > > > > ;The result is the square in the LEFT LOWER corner.
>
> > > > > I would be very grateful for comments/advices/solutions
>
> > > > I think you need to read the on-line help for the FFT function. :-)
>
> > > > Your code should look like this:
>
> > > > =A0 =A0FFTtest =3D FFT(ytest, -1)
> > > > =A0 =A0inv_test =3D Real_Part(FFT(FFTtest, 1))
>
> > > > Now ytest and inv_test are essentially the same.
>
> > > > Cheers,
>
> > > > David
>
> > > > --
> > > > David Fanning, Ph.D.
> > > > Fanning Software Consulting, Inc.
> > > > Coyote's Guide to IDL Programming:http://www.dfanning.com/
> > > > Sepore ma de ni thui. ("Perhaps thou speakest truth.")
>
> > > Hello,
> > > Im working with the FFT also but not directly in IDL, I work in ENVI.
> > > I wonder how ENVI or IDL can perform a FFT on a rectangular image?
> > > Theorically the image must be a square of dimension of coefficient of
> > > 2. (512x512, 1024 x1024 etc... ) So what is the process that make
> > > possible to do a FFT on a rectangular image? Other Image Processing
> > > software like PCI Geomatica cant do that.
>
> > Taking a DFT (Discrete Fourier Transformation) of an array is possible
> > for any array size. There is an algorithm called FFT (Fast Fourier
> > Transformation)
> > that happens to be very efficient if the size is in the form 2^N for
> > some N.
> > However, modern incarnations of FFT can deal with other sizes too,
> > albeit less
> > efficiently (the smaller the factors in the prime decomposition of the
> > size,
> > the better).
>
> > Please don't let what a particular software does or fail to do be your
> > guide to what is possible or not (from a mathematical standpoint). If
> > you
> > want to learn a bit more about the FFT, read for instance the chapter
> > about
> > it on the numerical recipes book.
>
> > Ciao,
> > Paolo
>
> > > Thank you,
> > > Max
>
> Hi Paolo,
>
> Thank you for your comment. As a student in Remote Sensing, I have a
> lot of basic things to understand. The image I would like to filter in
> the frequential domain (using a Butterworth filter) is a side scan
> sonar image mosaic of 9166 x 4093 pixels. Would you recommend a FFT
> and if yes would you use a special algorithm?
>
> greetings,
> Max
Well if it's just a single image, then you can certainly go ahead
and implement filtering with the FFT - the array is pretty big
and it will take a little while. If performance becomes an
issue, you can expand it to 9216x4096 and it should run a bit
faster.
Ciao,
Paolo
|
|
0
|
|
|
|
Reply
|
Paolo
|
12/9/2010 6:58:02 PM
|
|
|
8 Replies
816 Views
(page loaded in 0.137 seconds)
|