Number of points in FFT with dimension > 1
josh.liptonduffin
Thu, 09/27/2012 - 12:18 pm
I'm performing an FFT on a 2D set of real-valued data (and also a 3d volume). The 2D matrix has (2Nx2N) points, and the 3d volume has (2Nx2Nx2N) points.
In all cases, the first dimension (x axis) of the output FFT is truncated to N+1 points, while the remaining axes contain 2N points. It's to be expected that the algorithm should be truncating to one-sided for real-valued input data, but why is it only doing it in one of the two (or three) dimensions? The non-x axes appear to contain the mirror image of the frequency>0 data points with frequency<0.
I think the answer to your question is in the history of the FFT algorithm. Although some trace it back all the way to the early 1800 and Gauss, I think it is enough to go only as far as the 1960's and the Cooley and Tukey algorithm. The cost of computation and storage was such that they wanted to take every possible shortcut. In particular, when the input is real valued and the output is complex you need twice as much storage. You can almost get around that by using symmetry. I say "almost" because for N points at the input you need 1+N/2 complex outputs. As it happens though, the DC term has to have a zero imaginary part because it corresponds to the integral of the input and the input is real...
Your post says nothing about any symmetry in your higher dimensional data so I assume none exists. Without some symmetry you can't apply the efficiency above because the problem does not factor. To see this write your transform as a transform of a 1D function that has even and odd components. Next add a second dimension and you would see the problem right away. BTW: you can find a very nice discussion and diagrams of the 1D version in Bracewell's "The Fourier Transform and its Applications", McGraw-Hill 1965 pp. 14-15.
It is probably worth mentioning that in special cases where you have, for example, radial symmetry you are better off performing the transform differently.
A.G.
WaveMetrics, Inc.
September 27, 2012 at 01:49 pm - Permalink
In my particular case, the symmetry can be applied post-hoc via a series of mirror transforms in frequency space. Mostly what I want to understand is:
1) Why do the higher-order axes get mirrored (ie why are they also not truncated to N+1, as the x-axis is)?
2) From a practical standpoint, assuming I can spit out just the +x,+y,+z quadrant as a 3d wave with use of the duplicate with the "/R" switch, is there a usage of the "Concatenate" command that can allow me to efficiently join 3d waves together?
For example: let's say that I have a 3d wave corresponding to +x,+y,+z (dimensions NxNxN). Let's say I also know that the data are symmetric about the yz plane, so that I can duplicate this wave and call Reverse/DIM=0 on the new wave to have something corresponding to -x,+y,+z. What can I then do to efficiently join these two waves to form a new wave (dimensions 2NxNxN) that spans the full range of x? Can Concatenate do this?
September 27, 2012 at 02:02 pm - Permalink
So I suppose you tried to write out even/odd composition of a function in 2D and it did not answer the question why you can't fold the y-axis as well? Maybe we should start from something else: Lets assume that your 2D data are described by f(x,y). Can you find a pair of functions g() and h() such that: f(x,y)=g(x)*h(y)? If this condition holds you can factor out your transforms as two independent 1D transforms. If this condition does not hold you will have to look for computational shortcuts by taking advantage of any symmetry you may have at the input.
Here is an example:
Concatenate/np=2/o {ddd,eee}, zzz
A.G.
WaveMetrics, Inc.
September 27, 2012 at 03:58 pm - Permalink