FFTW input/output ordering troubles

[Kazik]

Hello,
I'm using FFTW3 library to compute Fourier transform of my data.
Yes, I have read FFTW3 manual. Still I have (probably humiliating)
question/problem.

Firstly, I should probably mention that I am not very familiar with
what you probably call "signal processing". I look at Fourier
transform more from mathematical or physical point of view. So, if I
have a function defined on some interval [-z_max, z_max] (position
space or time domain, whatever), it's natural for me that after
computation of Fourier transform I will have function on interval [-
k_z_max, k_z_max] (momentum space or frequency domain, whatever). Not
from 0 to k_z_max and then suddenly from -k_z_max to 0. Nevertheless,
I thought I have understood the ordering of FFTW. Until I wrote my
code.

So, to make my story short. I would like to compute FFT of one
dimensional array.
Firstly I've defined a grid:
########################################
/* Position grid-z-axis: */
double* z = new double[N];
double zmax = 30;
double dz = (2*zmax)/N;
for(int i=0; i<N; i++){
z = -zmax+(i*dz);
}

/* Momentum grid-z-axis: */
double* pz = new double[N];
double dpz = M_PI/zmax;
double pzmax = (N/2)*dpz;
for(int i=0; i<=N/2; i++){
pz = i*dpz;
}
for(int i=1; i<N/2; i++){
pz[i+N/2] = -pzmax+i*dpz;
}
########################################

Then, arrays and plans:
########################################
fftw_complex *psiX = (fftw_complex*) fftw_malloc(sizeof
(fftw_complex) * N);
fftw_complex *psiK = (fftw_complex*) fftw_malloc(sizeof
(fftw_complex) * N);
fftw_complex *psiX2 = (fftw_complex*) fftw_malloc(sizeof
(fftw_complex) * N);
fftw_plan XtoK_Plan = fftw_plan_dft_1d(N, psiX, psiK, FFTW_FORWARD,
FFTW_PATIENT);
fftw_plan KtoX2_Plan = fftw_plan_dft_1d(N, psiK, psiX2,
FFTW_BACKWARD, FFTW_PATIENT);
########################################

Subsequently I defined two routines:
########################################
/*------------------------------------------------------
* Routine FFTW_1d_Forward computes one-dimensional
* Fourier transform of input array ArrayX (X refers to
* position space). Output: ArrayXintoK (K refers to
* momentum space).
------------------------------------------------------*/
void FFTW_1d_Forward(fftw_complex* ArrayX, fftw_complex* ArrayXintoK,
fftw_plan XtoK_Plan, const double dz, const int N){
for(int i=0; i<N; i++){
ArrayX[0] = ArrayX[0];
ArrayX[1] = ArrayX[1];
}
fftw_execute(XtoK_Plan);
for(int i=0; i<N; i++){
ArrayXintoK[0] *= (dz*pow(-1,i));
ArrayXintoK[1] *= (dz*pow(-1,i));
}
}

/*------------------------------------------------------
* Routine FFTW_1d_Backward computes one-dimensional
* inverse Fourier transform of input array ArrayK (K
* refers to momentum space). Output: ArrayKintoX (X
* refers to position space).
------------------------------------------------------*/
void FFTW_1d_Backward(fftw_complex* ArrayK, fftw_complex* ArrayKintoX,
fftw_plan KtoX_Plan, const double dz, const int N){
for(int i=0; i<N; i++){
ArrayK[0] = ArrayK[0]*pow(-1,i);
ArrayK[1] = ArrayK[1]*pow(-1,i);
}
fftw_execute(KtoX_Plan);
for(int i=0; i<N; i++){
ArrayKintoX[0] /= (1.*N*dz);
ArrayKintoX[1] /= (1.*N*dz);
}
}
########################################

which I used in main() like that:
########################################
FFTW_1d_Forward(psiX, psiK, XtoK_Plan, dz, N);
FFTW_1d_Backward(psiK, psiX2, KtoX2_Plan, dz, N);
########################################

The multipilcations by "dz" and "N*dz" in routines FFTW_1d_Forward and
FFTW_1d_Backward are due to the normalisation of the FFT. The
multiplications by pow(-1,i) follow my experience with Numerical
Recipes (NR) code for FFT. In principle this "-1"s should govern that
the final result is ordered like [-k_z_max, ..., 0, ..., k_z_max], as
far as I understand what is going on.

The problem is that the result of the above code is wrong [in order to
learn FFTW I compare every FFTW-result with NR-result, and this NR-
code gives the correct result for 100%]. In particular I've checked
such a function:
########################################
for(int i=0; i<N; i++){
psiX[0] = exp(-z*z)*(z + 3.*z*z-7.*z*z*z
)*sin(z);
psiX[1] = 0.;
}
########################################

I have chosen such an ugly example on purpose. The nice examples (e.g.
Gaussians) are often misleading.
What I've checked/compared wit NR-result are:
1) FFTW_BACKWARD[ FFTW_FORWARD[ psiX ] ] = psiX
So, this is ok, as it should be.
2) The Plancherel theorem (the conservation of a norm of psiX and
FFTW_FORWARD[ psiX ]) also works perfectly
3) The real part of the FFT_FORWARD[psiX] is exactly the same as I
obtain form NR-program
4) !!! HERE COMES THE PROBLEM !!! :
the imaginary part of FFT_FORWARD[psiX] is mirror-reflected with
respect to "y"-axis (of course I know this by comaprison with NR-
result). So, if only I swapped the "right" and "left" part of the
imaginary part of the array FFT_FORWARD[psiX] (it's called psiK in the
code above), everything would be perfect. HOWEVER, if I did that,
point 1) crashes, i.e. I don't get identity when subsequently apply
FORWARD and BACKWARD.

I feel really lost with this. I've tried "experimentally" many other
ways of ordering, but every time some of the 4 points I mentioned
(actually points 1, 3 and 4) crashes. So, in general the result looks
"not-so-bad" (youi can definitely see that something works) but there
is always some mirror-reflection with respect to "y"-axis or the
result looks like it was globally multiplied by -1, and so on.

I would really appreciate any hints and/or advices. I spent all day
writing the code above but it seems that I failed.
Thank you in advance.

 

[Victor Bazarov]

Hello,
[..]
void FFTW_1d_Forward(fftw_complex* ArrayX, fftw_complex* ArrayXintoK,
fftw_plan XtoK_Plan, const double dz, const int N){
for(int i=0; i<N; i++){
ArrayX[0] = ArrayX[0];
ArrayX[1] = ArrayX[1];
}

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Uh... What exactly are you accomplishing by this?

fftw_execute(XtoK_Plan);
for(int i=0; i<N; i++){
ArrayXintoK[0] *= (dz*pow(-1,i));
ArrayXintoK[1] *= (dz*pow(-1,i));
}
}

[..]

V

 

[Kazik]

Hello,
[..]
void FFTW_1d_Forward(fftw_complex* ArrayX, fftw_complex* ArrayXintoK,
fftw_plan XtoK_Plan, const double dz, const int N){
  for(int i=0; i<N; i++){
    ArrayX[0] =  ArrayX[0];
    ArrayX[1] =  ArrayX[1];
  }

    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Uh...  What exactly are you accomplishing by this?

Gee, I'm sorry. It's a part of the old code. Of course it does
nothing, neither good nor bad.
What is important in routine FFTW_1d_Forward is this part:

for(int i=0; i<N; i++){
ArrayXintoK[0] *= (dz*pow(-1,i));
ArrayXintoK[1] *= (dz*pow(-1,i));
}

And one more correction to the previous post. In point 4) I have
written:

FFT_FORWARD[psiX] is mirror-reflected with
respect to "y"-axis (of course I know this by comaprison with NR-
result). So, if only I swapped the "right" and "left" part of the
imaginary part of the array FFT_FORWARD[psiX] (it's called psiK in the
code above), everything would be perfect.

Which is partially wrong. The resulting imaginary part of FFT_FORWARD
[psiX] is indeed mirror-reflected
with respect to "y"-axis, but then the "swapping" would not help for
that. What would help is a mirror
reflection with respect to "y"-axis.

I'm still hopeful for some suggestions. I really need help with this.

 

[Victor Bazarov]

Hello,
[..]
void FFTW_1d_Forward(fftw_complex* ArrayX, fftw_complex* ArrayXintoK,
fftw_plan XtoK_Plan, const double dz, const int N){
for(int i=0; i<N; i++){
ArrayX[0] = ArrayX[0];
ArrayX[1] = ArrayX[1];
}

^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Uh... What exactly are you accomplishing by this?

Gee, I'm sorry. It's a part of the old code. Of course it does
nothing, neither good nor bad.

By itself it does nothing. But _if_ it was *supposed* to be something
else, the existing code *pretends* to do something and hides the fact
that something else isn't actually performed. I don't know enough of
your algorithm (perhaps you should ask in 'sci.math' somewhere) to
verify, and that's why I pointed out that code to you.

What is important in routine FFTW_1d_Forward is this part:

for(int i=0; i<N; i++){
ArrayXintoK[0] *= (dz*pow(-1,i));
ArrayXintoK[1] *= (dz*pow(-1,i));
}

OK. Whatever.

And one more correction to the previous post. In point 4) I have
written:

FFT_FORWARD[psiX] is mirror-reflected with
respect to "y"-axis (of course I know this by comaprison with NR-
result). So, if only I swapped the "right" and "left" part of the
imaginary part of the array FFT_FORWARD[psiX] (it's called psiK in the
code above), everything would be perfect.

Which is partially wrong. The resulting imaginary part of FFT_FORWARD
[psiX] is indeed mirror-reflected
with respect to "y"-axis, but then the "swapping" would not help for
that. What would help is a mirror
reflection with respect to "y"-axis.

Let me understand... You're asking in a *C++ language forum* how to
perform a transform on a set of points with complex coordinates, yes?
Isn't this a math problem?

On a real plane to get (x2,y2) as a y-axis mirror image of (x1,y1) you'd do

x2 = -x1;
y2 = y1;

I'm still hopeful for some suggestions. I really need help with this.

Please consider asking in a math forum, like 'sci.math.num-analysis'.

V

 

[Kazik]

Ok, I'll post my question there as well. Thank you for the suggestion.

I'm asking here because it's purely programming problem. I know how
to
perform mirror reflection... I also calculate some Fourier transform
analytically on the piece of paper that lies just in front of me...
The
trouble is the order that FFTW uses. I would be perfectly happy if
could
prepare my data in human-readable way, execute the FFTW plan and then
read
the result. To make it clear: Imagine you have a function f(t) and you
plot it
on your paper sheet. Then somebody is asking you question "How does
the Fourier
transform of your function f(t) look?". What I would do is I would
draw (just
underneath) f(\omega) in such a way that the left end of \omega axis
correspond to
-Infinity and the right to +Infinity (yes, I know this is discrete FFT
and there
are no infinities here). But here it's not so easy in case of
numerical
FFT. Somehow I understand the way the routine from Numerical Recipes
works- there
you have to also make some additional "customization" if you want to
have the order
[-\omega_max,...,0, +\omega_max]. Since i thought the FFTW algorithm
works similarly,
I've written the code from my first post in a NR-like way. It doesn't
work.
Just as I've described.

Hoping for some solution,

 

[Kazik]

Sorry for the layout of the above...

 

[Victor Bazarov]

Ok, I'll post my question there as well. Thank you for the suggestion.

I'm asking here because it's purely programming problem.

Have you tried 'comp.programming'? Here we discuss purely C++ language
problems, not purely programming problems. Just so you know.

> I know how
to
perform mirror reflection... I also calculate some Fourier transform
analytically on the piece of paper that lies just in front of me...
The
trouble is the order that FFTW uses.

I take it that FFTW is some kind of library, yes? Doesn't it have its
own technical support or mailing list or online discussion forum?
Really, if we invite everybody who has a problem with a product they are
using while writing programs in C++ to 'comp.lang.c++', we would have no
room left for real C++ language discussions.

> [..]

V

 

[Jerry Coffin]

[ ... ]

What is important in routine FFTW_1d_Forward is this part:

for(int i=0; i<N; i++){
ArrayXintoK[0] *= (dz*pow(-1,i));
ArrayXintoK[1] *= (dz*pow(-1,i));
}

This looks pretty inefficient to me. At least normally, pow() will be
a very general purpose library function, but in this case you're
using a very special case: -1 raised to an even power is 1 and to an
odd power is -1. I'd try something like this instead:

for (int i=0; i<N-1; i+=2) {
ArrayXintoK[0] *= dz;
ArrayXintoK[1] *= dz;
ArrayXintoK[i+1][0] *= -dz;
ArrayXintoK[i+1][1] *= -dz;
}

// if N will never be odd, delete the following:
if (N&1) {
ArrayXintoK[N-1][0] *= dz;
ArrayXintoK[N-1][1] *= dz;
}

Of course, you'll need to do some testing to see whether this really
improves the speed, but I think there's a pretty fair chance it will.

 

[Steven G. Johnson]

The OP has posted to three different newsgroups, none of which is
really an appropriate one (if any, comp.dsp would be most appropriate,
as this is really not a problem with FFTW per se, but rather a
misunderstanding about the discrete Fourier transform). I won't
repeat my response here, except to say that this is just a difference
in sign convention between FFTW and Numerical Recipes, explained in
the FFTW FAQ, that leads to a flip in the imaginary part between the
two.

Regards,
Steven G. Johnson

 

[Steven G. Johnson]

The OP has posted to three different newsgroups, none of which is
really an appropriate one (if any, comp.dsp would be most appropriate,

I take that back; he also posted to comp.dsp too yesterday. The
irritating thing about this is the sense that, by responding, we are
rewarding the OP for spamming enough Usenet groups with separate
duplicate posts.