> for (tr = 0; tr < 40; tr++)
>    for (ch = 0; ch < 64; ch++)
>        for (col = 0; col < 20; col++)
>            for (n=0; n < 512; n++)
>                x = data[tr, ch, col, n];
>            y = fft(x);

I'm assuming that your array is stored in row-major order as
40x64x20x512.  Then the three loops you want correspond to:

howmany_rank = 3;
howmany_dims[0].n = 40;
howmany_dims[0].is = howmany_dims[0].os = 64*20*512;
howmany_dims[1].n = 64;
howmany_dims[1].is = howmany_dims[0].os = 20*512;
howmany_dims[2].n = 20;
howmany_dims[2].is = howmany_dims[0].os = 512;

Your transform is then 1d (rank=1) with dims[0].n = 512 and dims[0].is
= dims[0].os = 1.

Actually, in this particular case where you are transforming the *last*
dimension, the three loops can be coalesced into a single loop, but
FFTW will do this automatically internally.

See also the section in the FFTW manual on multi-dimensional array
format for more explanation:
    http://www.fftw.org/doc/Multi_002ddimensional-Array-Format.html

Steven G. Johnson

PS. It's more complicated if you have real data and want to use FFTW's
r2c plans.  In that case the strides are in units of real numbers
instead of complex numbes so you have to multiply by 2, and furthermore
the output is slightly bigger than the input so 512 is replaced by 514
to account for this as described in the manual (assuming you are
operating in-place).  See the manual.

I found another error; I'm using the r2c with real data, so the output
dimensions (the y array) should be:
[40, 64, 20, 257]
Also, the last line of pseudocode above should probably be (in matlab
syntax):

y(tr, ch, col, :) = fft(x);

So, the plan I'm using is:

plan = fftw_plan_guru_split_dft_r2c(
    4, dims, // rank and dims
    4, ???, // don't know how to setup howmany_dims
    x, y_re, y_im,
    FLAGS);

Hopefully this makes a little more sense now.
Adam

Just found a small error in the above pseudocode; here is what it
should have been:

data = double[40*64*20*512];
y = double[40*64*20*512];
...
for (tr = 0; tr < 40; tr++)
    for (ch = 0; ch < 64; ch++)
        for (col = 0; col < 20; col++)
            for (n=0; n < 512; n++)
                x = data[tr, ch, col, n];
            y = fft(x); 

Adam

Ok, I think I'm almost there. If I have a 4D array, with dimensions of:
[40, 64, 20, 512], and I want to calculate the FFT of the last
dimension, how do I set up the howmany_rank structure to do this?
The pseudocode for what I want to do would be:

data = double[40*64*20*512];
y = double[40*64*20*512];
...
for (tr = 0; tr < 40; tr++)
    for (ch = 0; ch < 64; ch++)
        for (col = 0; col < 20; col++)
            for (n=0; n < 512; n++)
                // this is obviously not how to access the data, but
this is for clarity
                x = data[tr, ch, col, n];
                y = fft(x);

How is the howmany_dims array configured to do this series of for
loops?

Thanks!

Adam

Steven,

I'll give your suggestions a closer look.
Thanks for the quick reply,

Adam

jawilson wrote:
> plan every time, but how is that accomplished? The way I'm envisioning
> it, I create the plan using the FFTW_MEASURE flag, then use
> fftw_execute repeatedly, changing the contents of the input array (not
> the pointer value, but the contents of the array) in between every
> execute.

Yes, that's how we envision FFTW being normally used when repeated FFTs
are required.

> Also, does the input array in the create plan function need to
> have valid data (ie the first chunk of data), or is it finding the
> optimal algorithm based on the size of the FFT alone?

The input to create_plan does not have to have any valid data, but it
does get overwritten during planning in FFTW_MEASURE or FFTW_PATIENT
mode.  This is described in the manual.

> I think I'm
> primarily getting confused by the wording of the FFTW manual, where it
> says things like "To apply a given plan to a different array, use can
> use the guru interface." By different array, does that mean that the
> original contents of x have changed, or that the pointer x itself has
> changed?

If you just change the contents of x, then you can use the same plan.
You can use the guru interface if you want to apply the plan to a
different array (different pointer) of the same size etc.

> I use the fft function like: fft(X, 512, 2), the function calculates
> the fft along the 2nd dimension, returning an array of size 100x512.
> (If I do fft(X, 512, 1) or fft(X,512), the returned answer is 512x50).
> Is this possible with FFTW?

FFTW won't do zero-padding for you, but it can transform just one
dimension of an array using the advanced or guru interface.

e.g. if you have an M x N array, stored in row-major order, and you
want to only FFT the columns (the N dimension), then you would call the
advanced interface to do an FFT of size N with howmany=M, stride=1, and
dist=N.  See the manual for more information.

> (3) Again, for the multidemensional functions, would it be possible to
> do the fft of an entire 4D array, only along a specified dimension?
> This would be similar to the matlab spectrogram, but not limited to a
> 2D array.

Yes.  In general if you have a 4D array and want to transform only an
"interior" dimension (i.e. not the first or last dimension), then you
need to use the guru interface.  The reason is that transforming an
interior dimension corresponds to a double loop of transforms, not a
single loop.  (Think about it.)

For an answer to a similar question, where someone wanted to transform
the middle dimension of a 3D array, see:

http://groups.google.com/group/comp.lang.fortran/msg/463235d14226650d

Cordially,
Steven G. Johnson

Hello,

I am new to using the FFTW algorithms in C++, and need some help
getting started. I have used matlab (and done a good deal of C++) for
many years, and I am rewriting some of my matlab code in C++ to see if
it helps with speed/memory issues. I'll give an overview of what I'm
trying to do, and see if anyone has suggestions.

So, I am reading in a MAT file, and the data I want is stored in a 3D
array. Matlab arrays are represented in fortran (column-first) format,
which is not what FFTW uses, but I think that can be overcome by
specifying the fftw plan slightly differently.

The data contains epoched EEG data, where the dimensions are [samples,
channels, trials]. For example, the array might consist of 40 trials,
each with 64 channels, and each of these with 3000 samples, therefore
with dimensions of [3000, 64, 40], for a total of 7,680,000 double
elements.

I want to calculate the spectrogram along the samples dimension; I am
using the matlab definition of spectrogram, where I specify the window
length, overlap, and size of the FFT. A moving window is applied,
taking the DFT for each window (padding zeros when necessary). The
dimensions of the output are then [freqs, samples, channels, trials],
where freqs is nfft/2+1, the times are the number of samples (dependent
on the window length and overlap), and channels and trials are
unchanged.

There are a few methods of doing this (I hope...)

(1) The initial thought on how to do this is to basically pound it out
in several for loops. For each trial and channel, I take winLen samples
(say 100), pad to nfft (say 512), do the DFT via FFTW, get the
magnitude and remove redundant data (so I'm left with 257 samples),
move ahead overlap samples, and wash-rinse-repeat. This is intuitively
the simplest, and probably where I want to start, but I do not know how
do accomplish this using the fftw interfaces. With FFTW, I know I
create a plan, in this case using fftw_plan_dft_r2c_1d, passing the fft
length, input array pointer, output array pointer, and flags. My
confusion is about the fact that I will be doing this thousands of
times on different data, and I think that I do not need to recreate the
plan every time, but how is that accomplished? The way I'm envisioning
it, I create the plan using the FFTW_MEASURE flag, then use
fftw_execute repeatedly, changing the contents of the input array (not
the pointer value, but the contents of the array) in between every
execute. Also, does the input array in the create plan function need to
have valid data (ie the first chunk of data), or is it finding the
optimal algorithm based on the size of the FFT alone? I think I'm
primarily getting confused by the wording of the FFTW manual, where it
says things like "To apply a given plan to a different array, use can
use the guru interface." By different array, does that mean that the
original contents of x have changed, or that the pointer x itself has
changed? If the size of the data is not changing, it seems pretty
inefficient to recalculate the plan for every step. (Again, I'm not
familiar with FFTW, so any example code that shows EXACTLY how to do
this would be very appreciated.)

(2) How do the multidimensional functions work? If you are familiar
with matlab, if I have a 2D array X, of size 100 x 50, for example, and
I use the fft function like: fft(X, 512, 2), the function calculates
the fft along the 2nd dimension, returning an array of size 100x512.
(If I do fft(X, 512, 1) or fft(X,512), the returned answer is 512x50).
Is this possible with FFTW? The spectrogram function works this way in
matlab, where it takes a 1D vector of data, creates a 2D array based on
the winLen and overlap, and takes the FFT down the winLen dimension.

(3) Again, for the multidemensional functions, would it be possible to
do the fft of an entire 4D array, only along a specified dimension?
This would be similar to the matlab spectrogram, but not limited to a
2D array.

I really appreciate anyone who takes the time to wade through this
series of questions. Any help you can give me is obviously of
tremendous help.

Thanks!
Adam