Hi,

I find that the reverse transform code in realft() does not always yield correct result. Here is my testing code:

//////////////////////
const int N = 8;
Vec_DP data1(N), data2(N);
// Prepare test data
for ( int i = 0; i < N; i ++ ) {
data[i] = sin(1.0*i);
data2[i] = 1.0*i;
}
// run FFT
NR::realft(data1, 1);
// result (it's correct)
//No. 0: 0.553733 - 0.807574 i
// No. 1: 2.39465 + 2.09701 i
// No. 2: -1.38668 - 0.91556 i
//No. 3: -0.881042 - 0.280414 i
NR::realft(data1, -1);
// result (it's wrong)
// 0, 0.262406, 1.48836, 0.14112, -0.756802,
// -0.379859,-0.858481,0.656987
NR::realft(data2, 1);
// result (it's correct)
//No. 0: 28 - 4 i
//No. 1: -4 - 9.65685 i
//No. 2: -4 - 4 i
//No. 3: -4 - 1.65685 i
NR::realft(data2, -1);
// result (it' correct for this test dataset!) 0,1,2,3,4,5,6,7

I also try different data sets such as cos, exp and so on and the reverse transform results are wrong.

I type the routine code of this chapter and am sure there is no type errors in the code as I've tested the routines in this chapter with different datasets. Everything is OK except the one I mentioned above.

I am blindly guessing that the issue may result from the special structure in which the complex data is stored after forward transform, i.e. the first element of the complex array contains the real part of first (0 frequency) and middle (Niquist frequency) values of the transformed array.
Anyway, is there anyone who's got the same experience? Any help will be greatly appreciated!

Lianqing
3/31/2006

And the convlv() in chapter 13 also suffers from this problem.

I have the same problem with realft.
The inverse ft doesn't work even if the direct ft works.

Willo

I did some modifications and now it's working for the inverse transform too.

void realft(float data[], unsigned long n, int isign)
{
void four1(float data[], unsigned long nn, int isign);
unsigned long i,i1,i2,i3,i4,np3;
float c1=0.5,c2,h1r,h1i,h2r,h2i;
double wr,wi,wpr,wpi,wtemp,theta;

theta=3.141592653589793/(double) (n>>1);
if (isign == 1) {
c2 = -0.5;
four1(data,n>>1,1);
} else {
c2=0.5;
theta = -theta;
four1(data,n>>1,-1);
}
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);
wr=1.0+wpr;
wi=wpi;
np3=n+3;
for (i=2;i<=(n>>2);i++) {
i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
h1r=c1*(data[i1]+data[i3]);
h1i=c1*(data[i2]-data[i4]);
h2r = -c2*(data[i2]+data[i4]);
h2i=c2*(data[i1]-data[i3]);
data[i1]=h1r+wr*h2r-wi*h2i;
data[i2]=h1i+wr*h2i+wi*h2r;
data[i3]=h1r-wr*h2r+wi*h2i;
data[i4] = -h1i+wr*h2i+wi*h2r;
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}
data[1] = (h1r=data[1])+data[2];
data[2] = h1r-data[2];
}

Don't forget that for the inverse transform result must be multiplied by 1/n.

Simon

Just to be clear: the NR code for realft is fine. It gives the correct results on the test problem given in the first post above. And one must multiply the inverse transform by 2/n (not 1/n) to get the original data. This is stated correctly in the comment to realft in the book.

Saul Teukolsky