DSPRelated.com
Forums

fftw: definition of the wave vector k in fourier space

Started by paquillo March 19, 2006
Hello everybody!
I am having problems with the defintion of the k vector trying to solve a
2 and 3 dimensional differential equation in Fourier space, namely:
div v_vector(x)=c*d(x).
This can be written in Fourier space like:
v_vector(k)=i**c*d(k)*k_vector/k^2.
I am doing the complex to complex fftw transform of d(x) adding  zero
imaginary part to d(x).
Now I would expect the following defintion of the k_vector to be correct:
	   for (int i=0 ; i<Nx;i++){
	       kx[i] = i;		
	     }	  	
	   for (int i=0 ; i<Ny;i++){
	       ky[i] = i;		
	     }
	   for (int i=1 ; i<=Nx/2;i++){
		   kx[Nx-i] = -kx[i];
	     }
	   for (int i=1 ; i<=Ny/2;i++){
		   ky[Ny-i] = -ky[i];
	     }
This means, that I have first the 0 frequency, then from 1 till the
nyquist frequency, N/2, and then from -N/2 til -1.
However, this defintions give wrong results (unreasonable results),
instead the following definition seems to work (very reasonable results)
and I do not understand why , since it should be completely wrong:
	   for (int i=0 ; i<Nx;i++){
	       kx[i] = i+1;		
	     }	  	
	   for (int i=0 ; i<Ny;i++){
	       ky[i] = i+1;		
	     }
	   for (int i=0 ; i<Nx/2;i++){
		   kx[Nx-1-i] = kx[i];
	     }
	   for (int i=0 ; i<Ny/2;i++){
		   ky[Ny-1-i] = ky[i];
	     }

Does somebody know which is the correct definition of the wave vector k in
Fourier space?

Thank you in advance!

PD:
The solution of the dif eq in 2d is then:

	 for (int i=0;i<Nx;i++){
	     for (int j=0;j<Ny;j++){		     
		 re(vx[j+Ny*i]) = c* kx[i]/(kx[i]*kx[i]+ky[j]*ky[j]) im(d[j+Ny*i]);	 
		 im(vx[j+Ny*i]) = c* kx[i]/(kx[i]*kx[i]+ky[j]*ky[j])-re(d[j+Ny*i]);	 
		 re(vy[j+Ny*i]) = c* ky[j]/(kx[i]*kx[i]+ky[j]*ky[j]) im(d[j+Ny*i]);	 
		 im(vy[j+Ny*i]) = c* ky[j]/(kx[i]*kx[i]+ky[j]*ky[j])-re(d[j+Ny*i]);	 
	       }
	   }



paquillo wrote:
> Now I would expect the following defintion of the k_vector to be correct: > for (int i=0 ; i<Nx;i++){ > kx[i] = i; > } > for (int i=1 ; i<=Nx/2;i++){ > kx[Nx-i] = -kx[i]; > }
This looks mostly fine at first glance, except kx[Nx/2] and ky[Ny/2] are not defined (for even Nx and Ny): you need a "<=" on your first loop. Note also that you should really use 2*pi*i/(Nx*dx) where dx is your discretization in the x direction, and similarly for y; this matters if your computational cell has a different size in the x and y directions since otherwise your kx and ky won't have compatible units. Also, maybe you forgot to special-case the (0,0) point when you are solving the PDE, since if you are not careful you'll end up dividing by zero there. Cordially, Steven G. Johnson
Thank you for your comment Steven.
however let me say you that I am still confused.

1)"kx[Nx/2] and ky[Ny/2] are not defined (for even Nx and Ny)", why, what
can I do?

2)Now, d(x) is a Nx*Ny=Nx*Nx array: d(0) ... d(Nx*Nx-1)
The corresponding kx component should then go from kx(0) till kx(Nx-1), so
I do not understand why you tell me to put <= in the "first loop".

3)I imposed Nx=Ny (always even) with an equidistant grid, and thus kx=ky.
The normalization 2*pi/(Nx*dx) should not affect the qualitative result.
Besides that I thought this normalization is unnecessary since it is
carried in the exponent in the definition of the fftw.


4)I fully agree with you in the final comment. I took this into account.
For the first definition of the k_vector I imposed v=0 if k=0. For the
second definition this does not happen (k=!0).

Cordially,
paco


Sorry, I misread your code, I thought your first loop went from i <
Nx/2 when in fact you used i < Nx.

Dear community, I got the bug!

The first definition I gave was indeed correct besides of a factor
cx=2pi/N, as Steven pointed out. Forget about the second definition.
The bug was in the sign of the solution of the dif eq that I wrote in the
PD. This sign appears because of the nabla operator: -ik, and has to be in
agreement with the fftw definition of the FFT and not with your own
one!!!!

Thank you everybody!
Cheers!

Correct definition:
> for (int i=0 ; i<Nx;i++){ > kx[i] = cx*i; > } > for (int i=0 ; i<Ny;i++){ > ky[i] = cy*i; > } > for (int i=1 ; i<=Nx/2;i++){ > kx[Nx-i] = -kx[i]; > } > for (int i=1 ; i<=Ny/2;i++){ > ky[Ny-i] = -ky[i]; > }
Correct solution:
>PD: >The solution of the dif eq in 2d is then: > > for (int i=0;i<Nx;i++){ > for (int j=0;j<Ny;j++){ > re(vx[j+Ny*i]) = c* kx[i]/(kx[i]*kx[i]+ky[j]*ky[j])
*(-1)*im(d[j+Ny*i]);
> im(vx[j+Ny*i]) = c* kx[i]/(kx[i]*kx[i]+ky[j]*ky[j])*re(d[j+Ny*i]); > re(vy[j+Ny*i]) = c* ky[j]/(kx[i]*kx[i]+ky[j]*ky[j])
*(-1)*im(d[j+Ny*i]);
> im(vy[j+Ny*i]) = c* ky[j]/(kx[i]*kx[i]+ky[j]*ky[j])*re(d[j+Ny*i]); > } > } > > > >