Using *real* FFT separably for multiple dimensions

Hello Comp.DSP,

I need to perform forward and inverse FFT of multi-dimensional, real data.
I am already successful to do this using the complex FFT separably.

Now I would like to do this using the real FFT for computational savings.

However, after I perform real FFT of the first dimension, the data is
complex and it seems I must use the complex FFT for subsequent
dimensions.

It makes sense, because the output of the real FFT is complex data, but is
it really impossible to use the real FFT separably for multiple
dimensions?

Thank you,

Anthony

On 2007-12-13 09:13:42 -0400, "dee.ess.pee" 
<dee.ess.pee.for.mee@gmail.com> said:

> Hello Comp.DSP,
> 
> I need to perform forward and inverse FFT of multi-dimensional, real data.
> I am already successful to do this using the complex FFT separably.
> 
> Now I would like to do this using the real FFT for computational savings.
> 
> However, after I perform real FFT of the first dimension, the data is
> complex and it seems I must use the complex FFT for subsequent
> dimensions.
> 
> It makes sense, because the output of the real FFT is complex data, but is
> it really impossible to use the real FFT separably for multiple
> dimensions?
> 
> Thank you,
> 
> Anthony

The symmetries present have nothing to do with whether you are
using a computational efficient version (ie FFT) of the DFT or
not. Separability of FTs is a very useful feature that you
are correct to use.

The symmetry that arises from the data being real is more
subtle in higher dimensions so you can not "just" use the real
FT (i.e. the symmetry suitable for one dimensional forms) in
higher dimensions.

There is symmetry in the intermediate arising after applying
the FT over one dimension as well as in the final result after
applying the FT over both dimensions. The details are not hard
to determine. Just different.

In one dimension the symmetry strongly suggests a suitable form
of the real DFT which is quite commonly used. In two or more
dimensions the suggestion is rather weaker and there is no common
form to compactly represent the results.

The issue of symmetry of FTs is the "bread and butter" of
crystalography. There are catalogues of the possible symmetries in
both two and three dimensions. In their terms you would seem to have
the simplest case of no spatial symmetry with just real data.

In article <PrmdnaQb1uGbrfzanZ2dnUVZ_q6mnZ2d@giganews.com>,
 "dee.ess.pee" <dee.ess.pee.for.mee@gmail.com> wrote:

> Hello Comp.DSP,
> 
> I need to perform forward and inverse FFT of multi-dimensional, real data.I 
> am already successful to do this using the complex FFT separably.
> 
> Now I would like to do this using the real FFT for computational savings.
> 
> However, after I perform real FFT of the first dimension, the data iscomplex 
> and it seems I must use the complex FFT for subsequentdimensions.
> 
> It makes sense, because the output of the real FFT is complex data, but isit 
> really impossible to use the real FFT separably for multipledimensions?


For the 2-D case, with M rows and N columns, take the 1-D real FT of the 
rows first.  The row data is Hermitian and you only need to keep half of 
it.  Then the 1-D complex FTs along the columns are on vectors of length 
M/2.  Thus you saved a factor of 2 in each direction for a total savings 
of 4.

KP

>For the 2-D case, with M rows and N columns, take the 1-D real FT of the 
>rows first.  The row data is Hermitian and you only need to keep half of

>it.  Then the 1-D complex FTs along the columns are on vectors of length

>M/2.  Thus you saved a factor of 2 in each direction for a total savings

>of 4.

Hello Gordon and Ken,

Thanks for your information!

I would like to consider the 2D case first. After the N length real FFT of
each row, I see a special case for the DC bin (0) and the Nyquist bin
(N/2), and I see conjugate symmetry between bins (1 ... N/2-1) and (N/2+1
... N-1).

But I'm not sure how to exploit this symmetry for the N/2 length complex
FFT of each column, and I wasn't able to find any reference for this
technique. Do you know a reference for this technique? I just need a
little more details.

Thank you,

Anthony

On 2007-12-16 23:50:11 -0400, "dee.ess.pee" 
<dee.ess.pee.for.mee@gmail.com> said:

>> For the 2-D case, with M rows and N columns, take the 1-D real FT of the
>> rows first.  The row data is Hermitian and you only need to keep half of
> 
>> it.  Then the 1-D complex FTs along the columns are on vectors of length
> 
>> M/2.  Thus you saved a factor of 2 in each direction for a total savings
> 
>> of 4.
> 
> Hello Gordon and Ken,
> 
> Thanks for your information!
> 
> I would like to consider the 2D case first. After the N length real FFT of
> each row, I see a special case for the DC bin (0) and the Nyquist bin
> (N/2), and I see conjugate symmetry between bins (1 ... N/2-1) and (N/2+1
> .. N-1).
> 
> But I'm not sure how to exploit this symmetry for the N/2 length complex
> FFT of each column, and I wasn't able to find any reference for this
> technique. Do you know a reference for this technique? I just need a
> little more details.
> 
> Thank you,
> 
> Anthony

Knuth quotes Hoare as "Premature optimization is the root of all evil".

Try working through things without the tricks and then work out what
the tricks are once you understand things.

You start with N point data. Where does N/2 come from? Sounds like a
trick that you are using but which is now sevring to block your progress.

Sometimes the symmetry is within a column as when it has real data.
Sometimes the symmetry is between two columns so only one needs to be
processed to determine the result for both. Or you may be dealing with
rows instead of columns.

In article <ApKdnaIXZbZubPjanZ2dnUVZ_vumnZ2d@giganews.com>,
 "dee.ess.pee" <dee.ess.pee.for.mee@gmail.com> wrote:

> >For the 2-D case, with M rows and N columns, take the 1-D real FT of the 
> >rows first.  The row data is Hermitian and you only need to keep half of
> 
> >it.  Then the 1-D complex FTs along the columns are on vectors of length
> 
> >M/2.  Thus you saved a factor of 2 in each direction for a total savings
> 
> >of 4.
> 
> Hello Gordon and Ken,
> 
> Thanks for your information!
> 
> I would like to consider the 2D case first. After the N length real FFT of
> each row, I see a special case for the DC bin (0) and the Nyquist bin
> (N/2), and I see conjugate symmetry between bins (1 ... N/2-1) and (N/2+1
> ... N-1).
> 
> But I'm not sure how to exploit this symmetry for the N/2 length complex
> FFT of each column, and I wasn't able to find any reference for this
> technique. Do you know a reference for this technique? I just need a
> little more details.
> 
> Thank you,
> 
> Anthony


For each row, put bin 0 at location 0, bin N/2 at location 1, bin 1 
re+im at locations 2 and 3, ..., bin N/2-1 re+im at locations N-2 and 
N-1.

Now your 2-D array has all reals in column 0, reals in column 1, real in 
column 2, imaginary in column 3, etc..

Next take real-FFTs of length M on columns 0 and 1.  Store the result 
along each column just like the rows: bin 0 at location 0, bin N/2 at 
location 1, etc..

Take complex-FFTs of length M on pairs of columns (2,3), (4, 5), ..., 
(N-2, N-1).  Store these as before, with the real part in one colums and 
the imaginary part in another.

The best thing for you to do is to sketch this out.  You will see that 
it works out.

Ken P.

On Dec 14, 8:06 am, Ken Prager <pra...@ieee.org> wrote:
> In article <PrmdnaQb1uGbrfzanZ2dnUVZ_q6mn...@giganews.com>,
>
>  "dee.ess.pee" <dee.ess.pee.for....@gmail.com> wrote:
> > Hello Comp.DSP,
>
> > I need to perform forward and inverse FFT of multi-dimensional, real data.I
> > am already successful to do this using the complex FFT separably.
>
> > Now I would like to do this using the real FFT for computational savings.
>
> > However, after I perform real FFT of the first dimension, the data iscomplex
> > and it seems I must use the complex FFT for subsequentdimensions.
>
> > It makes sense, because the output of the real FFT is complex data, but isit
> > really impossible to use the real FFT separably for multipledimensions?
>
> For the 2-D case, with M rows and N columns, take the 1-D real FT of the
> rows first.  The row data is Hermitian and you only need to keep half of
> it.  Then the 1-D complex FTs along the columns are on vectors of length
> M/2.  Thus you saved a factor of 2 in each direction for a total savings
> of 4.
>
> KP

Your conclusion is wrong: The total savings are only a factor of 2
(roughly) in time and storage compared to the complex 2d DFT.  Also,
your description of the column FFTs is wrong: they are of length M,
not M/2, but there are only about N/2 columns to transform.

You save a factor of 2 on the real row DFTs compared to complex row
DFTs, and you save a factor of 2 on the column FFTs (which are complex
of length M, but there only about N/2 of them because you only needed
to store half the outputs from each row FFT).  These factors of 2 are
*not* multiplicative---the overall savings is still a factor of about
2.

Regards,
Steven G. Johnson

On Dec 13, 8:58 am, Gordon Sande <g.sa...@worldnet.att.net> wrote:
> There is symmetry in the intermediate arising after applying
> the FT over one dimension as well as in the final result after
> applying the FT over both dimensions. The details are not hard
> to determine. Just different.
>
> In one dimension the symmetry strongly suggests a suitable form
> of the real DFT which is quite commonly used. In two or more
> dimensions the suggestion is rather weaker and there is no common
> form to compactly represent the results.

As you know, in principle you could store the results using exactly
half the storage of a corresponding complex DFT.  However, as you also
know, the storage format required to do such "packing" becomes
increasingly convoluted in higher dimensions.

Our experience for a long while has been that it just isn't worth it
to use these complicated packings to save *exactly* a factor of two in
storage.

When you transform the first (real) dimension of length N, just store
the non-redundant half as a complex array of length N/2+1, i.e.
storing the 0th and N/2th bins as separate complex numbers even though
they are real.  This wastes some space, but the wastage (a fraction 2/
N) is negligible in practice.

Then, the subsequent transform dimensions are all just ordinary
complex DFTs, but there are roughly half as many of them.  This is not
only easier to implement efficiently, it is also much easier to use,
because the result is just an ordinary complex multidimensional array
with one dimension N replaced by N/2+1, and no funny packing to
interpret.

(This is the format for multidimensional DFTs that we adopted in
FFTW.  It can be used in-place if you trivially pad each length-N row
of your real input with two extra unused storage slots.)

Regards,
Steven G. Johnson

Gordon - Re: premature optimization

This is good advice. First I have the multiple dimensions transform
working using complex FFT without optimizations, then I have the first
dimension working using the real FFT optimization, now I want to do the
later dimensions optimization. The N/2 length transform is a reference to
the later dimensions optimization.

Ken - Re: N/2 length with packing method

Thank you for the extra information, now I understand this optimization
and I am successful to implement it. Requiring exactly N/2 space, this
method is elegant for 2D transforms.

Steven - Re: N/2 + 1 length without packing method

When I try to extend the N/2 length with packing method to higher
dimensions, I see the complexity problem. Now I have implemented the N/2 +
1 length without packing method, and the code is elegant to generalize to
higher dimensions.

Thanks to all for your advice!

Anthony

In article 
<cd8df422-4130-4302-aa6c-f1a1fffc232e@e23g2000prf.googlegroups.com>,
 "Steven G. Johnson" <stevenj@alum.mit.edu> wrote:

> On Dec 14, 8:06 am, Ken Prager <pra...@ieee.org> wrote:
> > In article <PrmdnaQb1uGbrfzanZ2dnUVZ_q6mn...@giganews.com>,
> >
> >  "dee.ess.pee" <dee.ess.pee.for....@gmail.com> wrote:
> > > Hello Comp.DSP,
> >
> > > I need to perform forward and inverse FFT of multi-dimensional, real 
> > > data.I
> > > am already successful to do this using the complex FFT separably.
> >
> > > Now I would like to do this using the real FFT for computational savings.
> >
> > > However, after I perform real FFT of the first dimension, the data 
> > > iscomplex
> > > and it seems I must use the complex FFT for subsequentdimensions.
> >
> > > It makes sense, because the output of the real FFT is complex data, but 
> > > isit
> > > really impossible to use the real FFT separably for multipledimensions?
> >
> > For the 2-D case, with M rows and N columns, take the 1-D real FT of the
> > rows first.  The row data is Hermitian and you only need to keep half of
> > it.  Then the 1-D complex FTs along the columns are on vectors of length
> > M/2.  Thus you saved a factor of 2 in each direction for a total savings
> > of 4.
> >
> > KP
> 
> Your conclusion is wrong: The total savings are only a factor of 2
> (roughly) in time and storage compared to the complex 2d DFT.  Also,
> your description of the column FFTs is wrong: they are of length M,
> not M/2, but there are only about N/2 columns to transform.

The hazzard of responding too hastily--you are right and I corrected 
this in a subsequent post.


> You save a factor of 2 on the real row DFTs compared to complex row
> DFTs, and you save a factor of 2 on the column FFTs (which are complex
> of length M, but there only about N/2 of them because you only needed
> to store half the outputs from each row FFT).  These factors of 2 are
> *not* multiplicative---the overall savings is still a factor of about
> 2.

Right again.  This time because of the separability of the 2D DFT.  For 
those viewing from home, we went from...

M[5Nlog2(N)] + N[5Mlog2(M)]

to

1/2M[5Nlog2(N)] + 1/2N[5Mlog2(M)] = 1/2{M[5Nlog2(N)] + N[5Mlog2(M)]}

Regards,

Ken Prager