| Re: IDL 5.5, 2D FFT indexing confusion. [message #44821] |
Tue, 19 July 2005 13:49  |
R.G. Stockwell
Messages: 363
Registered: July 1999
|
Senior Member |
|
|
"Pitufa" <c.c.calderon@gmail.com> wrote in message
news:1121772201.952005.96070@g43g2000cwa.googlegroups.com...
> Hi,
>
> I have been trying to generate an real even function in fourier space
> that I can INVERSE FFT in order to get a function which has no
> imaginary part. I have no problems when the function is a vector, but I
> get an imaginary part when it is a two dimensional array.
...
Hello Pitufa,
I shrunk the size of the array and took a look.
I think your nyquist rows and columns should all be positive (i.e. don't
flip the signs).
(by nyquist rows/columns i mean the npix/2+1 columm and the npix/2+1 row)
Here is a npix=6 example that i fixed
f = [ $
[0.00, 0.00, 0.00, 0.00, 0.00, 0.00],$
[0.00, 1.00, 0.80, 0.60, -0.80, -1.00],$
[0.00, 0.80, 1.00, 0.92, -1.00, -0.80],$
[0.00, 0.60, 0.92, 1.00, 0.92, 0.60],$
[0.00,-0.80, -1.00, 0.92, 1.00, 0.80],$
[0.00,-1.00, -0.80, 0.60, 0.80, 1.00] $
]
Cheers,
bob
PS
Note that if your "npix" is odd, you have both a positive and negative
nyquist points
and they are both complex.
|
|
|
|
| Re: IDL 5.5, 2D FFT indexing confusion. [message #44825 is a reply to message #44821] |
Tue, 19 July 2005 09:44  |
Pitufa
Messages: 8
Registered: July 2005
|
Junior Member |
|
|
Thanks for your reply, I was wondering if there was a way to shrink the
phi derivation!
About the symmetry in fft, here is why I thought it had to be point
symmetric about the centre. The FFT of a function f_{kj} is given by
(for a square array of side dimension N):
F_{m,p} = (1/N^2) sum_{k,j} f_{k,j} exp[-2pi i (km + jp)/N] eqn.
[1]
Now, if the array in fourier space has its origin at the centre of the
array, then the point (m, p) is centrally opposite to (N-m, N-p). And
the FFT for this point is:
F_{N-m,N-p} = (1/N^2) sum_{k,j} f_{k,j} exp[-2pi i (k(N-m) + j(N-p))/N]
= (1/N^2) sum_{k,j} f_{k,j} exp[2pi i (km + jp)/N] exp[-2pi i (k+ j)]
= (1/N^2) sum_{k,j} f_{k,j} exp[2pi i (km + jp)/N] eqn. [2]
which is the complex conjugate of eqn [1] if f_{k,j} is real.
Please let me know if you don't agree.
Thanks,
Pitufa.
|
|
|
|
| Re: IDL 5.5, 2D FFT indexing confusion. [message #44830 is a reply to message #44825] |
Tue, 19 July 2005 07:29 |
Benjamin Hornberger
Messages: 258
Registered: March 2004
|
Senior Member |
|
|
Pitufa wrote:
> Hi,
>
> I have been trying to generate an real even function in fourier space
> that I can INVERSE FFT in order to get a function which has no
> imaginary part. I have no problems when the function is a vector, but I
> get an imaginary part when it is a two dimensional array.
>
> Below is the test program that shows my problem. Am I defining wrong
> the variables 'centre' or 'nshift'?
>
> I would be really grateful if someone could let me know what I am doing
> wrong.
>
> Thanks!
>
> Pitufa.
I believe your TEST array is not even in the sense which is required for
the FT to be real. I think it has to be even (symmetric) in each
dimension separately, and not point-symmetric about the center. Please
correct me if I'm wrong. E.g., the following TEST array gives a real FT:
test = sqrt((rebin(f, 100, 100))^2+$
(rebin(reform(f, 1, 100), 100, 100))^2)
which is basically the same as a shifted DIST array.
With respect to centers and shifts, I prefer to shift everything by n/2
when going from real to Fourier space and -n/2 when going back:
fourier = shift(fft(real), n/2)
real = fft(shift(fourier, -n/2), /inverse)
I this case, the zero frequency is right in the center for an odd number
of pixels and just right/above the center for an even number of pixels,
and the frequency array goes from -Nyquist to +(1 index below Nyquist).
For even number of pixels, you don't even have to worry about shifting
plus or minus any more.
As a side note, you can easily vectorize the calculation of your PHI array:
nx2d = rebin(indgen(nx)-cx, nx, ny)
ny2d = rebin(reform(indgen(ny)-cy, 1, ny), nx, ny)
phi = atan(ny2d, nx2d)
where nx, ny are the numbers of pixels and cx, cy the centers.
Benjamin
>
> pro test_index
>
> npix = 100
> centre = npix/2.d - 1.d
> nshift = npix/2 + 1
>
> ;1d example:
>
> f = abs(findgen(npix) - centre)
> ifft = fft(shift(f,nshift),1,/double)
>
> print,'1D: Imaginary maximum:', Max(Abs(Imaginary(ifft)))
> print,'1D: Real maximum :', Max(Abs(Double(ifft)))
>
> ;2d example:
>
> ;angle of position vector w.r.t x axis:
>
> phi = dblarr(NPIX,NPIX)
> FOR X = 0, NPIX-1 DO FOR Y = 0, NPIX-1 DO $
> PHI[X,Y] = ATAN(Y*1.D - centre,X*1.D - centre)
>
> TEST = sin(2.d*phi)
> IFFTTEST = FFT(SHIFT(TEST,nshift,nshift),1,/DOUBLE)
>
> print,'2D: Imaginary maximum: ', $
> MAX(ABS(imaginary(IFFTTEST))),mean(ABS(imaginary(IFFTTEST)))
> print,'2D: Real maximum : ',MAX(ABS(double(IFFTTEST)))
>
> end
>
|
|
|
|
| Re: IDL 5.5, 2D FFT indexing confusion. [message #44916 is a reply to message #44821] |
Thu, 21 July 2005 10:28 |
Pitufa
Messages: 8
Registered: July 2005
|
Junior Member |
|
|
Thanks Bob!
I tried what you told me and it seemed to work, although I am not
totally sure why..
I agree, everything becomes much simpler if N is odd, since the range
of frequencies is symmetric (taking centre as (npix - 1)/2.d and nshift
as (npix - 1)/2 +1).
I tried to understand how it can work with an even N, and got the
following conditions:
1. F_{m,p} = F_{N-m,N-p}*
2. F_{m,N} = F_{N-m,N}*
3. F_{N,p} = F_{N,N-p}*
4. F_{N,N} = F_{N/2-1,N/2-1}*
5. F_{N,N/2-1} = F_{N/2-1,N}*
but then it only works if the fourier array is real, and 2-5 don't make
any sense physically! My conclusion from this is that you are better
off using an odd N if you are doing this sort of calculation.
Agree? Disagree?
Why the nyquist rows/columns need to be positive?
Thanks alot,
Pitufa.
R.G. Stockwell wrote:
> "Pitufa" <c.c.calderon@gmail.com> wrote in message
> news:1121772201.952005.96070@g43g2000cwa.googlegroups.com...
>> Hi,
>>
>> I have been trying to generate an real even function in fourier space
>> that I can INVERSE FFT in order to get a function which has no
>> imaginary part. I have no problems when the function is a vector, but I
>> get an imaginary part when it is a two dimensional array.
> ...
>
> Hello Pitufa,
> I shrunk the size of the array and took a look.
> I think your nyquist rows and columns should all be positive (i.e. don't
> flip the signs).
> (by nyquist rows/columns i mean the npix/2+1 columm and the npix/2+1 row)
>
> Here is a npix=6 example that i fixed
> f = [ $
> [0.00, 0.00, 0.00, 0.00, 0.00, 0.00],$
> [0.00, 1.00, 0.80, 0.60, -0.80, -1.00],$
> [0.00, 0.80, 1.00, 0.92, -1.00, -0.80],$
> [0.00, 0.60, 0.92, 1.00, 0.92, 0.60],$
> [0.00,-0.80, -1.00, 0.92, 1.00, 0.80],$
> [0.00,-1.00, -0.80, 0.60, 0.80, 1.00] $
> ]
>
> Cheers,
> bob
>
> PS
> Note that if your "npix" is odd, you have both a positive and negative
> nyquist points
> and they are both complex.
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
