This code illustrates how to find the different momentum components in
the output of fftrn. The initial array of size (,
) is
generated with only two sinusoidal modes excited, the real part of
and the imaginary part of
. The loop at
the end shows how to find the different momenta. The excitation
only shows up once but the mode
and its complex
conjugate
are redundantly included in the output. Note that
no matter what you put in as input the output components
,
,
and
will be real.
#define N1 8 // Number of real points in the data - first dimension #define N2 4 // Second dimension void main() { int a,b,ka,kb; float f[N1][N2]; // An array of real numbers float fnyquist[2*N1]; // For storing k_2 = N2/2 frequencies int size[2]={N1,N2}; // Size of data in all dimensions float pi=2.*asin(1.),x1=2.*pi/(float)N1,x2=2.*pi/(float)N2; for(a=0;a<N1;a++) for(b=0;b<N2;b++) f[a][b]=5.*cos(3.*x1*(float)a + x2*(float)b)+sin(2.*x1*(float)a); fftrn((float *)f,(float *)fnyquist,2,size,1); for(a=0;a<N1;a++) { ka = (a<=N1/2 ? a : a-N1); // ka is the frequency stored at position a for(b=0;b<N2/2;b++) { kb = b; // kb is the frequency stored at position b printf("Frequency={%d,%d}, Transform value={%f + i %f}\n",ka,kb,f[a][2*b],f[a][2*b+1]); } printf("Frequency={%d,%d}, Transform value={%f + i %f}\n",ka,N2/2,fnyquist[2*a],fnyquist[2*a+1]); } } OUTPUT: Frequency={0,0}, Transform value={-0.000002 + i 0.000000} Frequency={0,1}, Transform value={0.000000 + i 0.000003} ... Frequency={2,0}, Transform value={0.000002 + i 15.999999} ... Frequency={3,1}, Transform value={79.999992 + i -0.000023} ... Frequency={-2,0}, Transform value={0.000002 + i -15.999999} ... Frequency={-1,1}, Transform value={0.000004 + i 0.000004} Frequency={-1,2}, Transform value={-0.000000 + i 0.000001}