next up previous
Next: How FFTEASY works Up: Storage: What Goes Where Previous: Output of Fourier Transforms

Output of Fourier Transforms of Real Data

For real arrays there is a subtlety in the arrangement of the output. Of course a real function is just a special case of a complex function and you could always take the Fourier transform of an array of real numbers by putting it in an array twice as big with every other component set to $0$ and then submitting it to fftc1 or fftcn. This method is obviously quite wasteful. The point of a Fourier transform routine for real data is to take advantage of the fact that N independent input points only produce N independent output points, so even though the result $F(k)$ consists of complex numbers they can all be stored in the original array. The reason this can be done is because the Fourier transform of a real set of data obeys the symmetry $F(-k)=F^*(k)$, so it's redundant to store the negative frequency output. You still need the values $F(0)$ and $F(N/2)$ but these are both real so they only take up half as much storage as the other frequencies. The net result is that $N$ real inputs produce $N$ real outputs in the form of $N/2-1$ complex numbers and two real ones.

For fftr1 this makes the output scheme very simple. The first two numbers in the output array are the real values $F(0)$ and $F(N/2)$. The rest of the values contain alternate real and imaginary parts for $F(1)$ through $F(N/2-1)$.

For fftrn the story becomes more complicated. If I let $k_{ndims}$ stand for the momentum in the last dimension of the array then I don't need to output any negative values of $k_{ndims}$; they can all be recovered by symmetry. On the other hand the values at $k_{ndims}=0$ and $k_{ndims}=N/2$ are in general no longer real. It is still true that the number of independent outputs equals the number of inputs, but there is no obvious scheme for packing all of the outputs into the original array. Rather than coming up with a complicated scheme that would make it difficult to find different frequency components I simply let the entire plane of values $k_{ndims}=0$ be represented in the original array and store all the values with $k_{ndims}=N/2$ in a separate array $fnyquist$. Thus if you look for the values $F(k_1,k_2,...,0)$ and $F(-k_1,-k_2,...,0)$ you will find them both in the output and they will be complex conjugates of each other.

So after all is said and done the output of fftrn is arranged with the last dimension going up from $0$ to $N/2-1$ and all the others wrapping around as described for the complex routines above. The values with $k_{ndims}=N/2$ are stored in the array $fnyquist[N_1][N_2]...[2 N_{ndims-1}]$.

If you have generated an array using fftr1 or fftrn then you can simply give that array back to the same function to do an inverse transform and recover your original data. If, however, you want to generate data in frequency space and then transform it to position space you have to be careful to ensure that your input obeys the symmetry $F(-k)=F^*(k)$. For fftr1 this simply means generating real numbers for $F_0$ and $F_{N/2}$ and putting them in the first two slots of the array and then filling all the subsequent slots with complex numbers. For fftrn, however, you have to make sure that for all the values with $k_n=0$ (the first plane of values in your lattice) this symmetry is obeyed. You then have to fill fnyquist with the values with $k_n=N/2$, obeying the same symmetry conditions.


next up previous
Next: How FFTEASY works Up: Storage: What Goes Where Previous: Output of Fourier Transforms

Go to The FFTEASY Home Page
Go to Gary Felder's Home Page
Send email to Gary at gfelder@email.smith.edu

This documentation was generated on 2003-09-30