realfft(): Odd result in comparison with Matlab?


jhpark
04-22-2006, 02:40 AM
Hi all,
I would like to post this to ask you help me out both the implementation or theory understanding.

First, Let's assume that the four1() routine has been sucessfully implemented. To check it, we deal with the following data:
Vec_IO_DP v="0 1 2 3 4 5 6 7";
// which corresponding to a complex sequence
c_v="0+1i 2+3i 4+5i 6+7i"
The output of four1 of v:
four1(v,1): 12 16 -8 0 -4 -4 0 -8
// which corresponding to a complex sequence
c_v1=12+16i -8+0i -4-4i 0-8i
And the fft(c_v) produced by Matlab is
c_v=[0+1i 2+3i 4+5i 6+7i];
>> b1=fft(c_v)

b1 =

12.0000 +16.0000i -8.0000 -4.0000 - 4.0000i 0 - 8.0000i
As comparison of c_v1 and b1, we can be asured that four1() is successfully implemented.

Now, first we check the theory for consituting the realfft(). For sake of convenience in explanation, we re-denote as following:
v(0) v(1) v(2) v(3) v(4) v(5) v(6) v(7)
After call four1(v,1) we have the sequence in v as:
Re{H0} Im{H0} Re{H1} Im{H1} Re{H2} Im{H2} Re{H3} Im{H3}
Let's examine the above real sequence v as an example. After calling four1(v,1) the sequence v will contain:
v(0)=Re{H0}=12
v(1)=Im{H0}=16
v(2)=Re{H1}=-8
v(3)=Im{H1}=0
v(4)=Re{H2}=-4
v(5)=Im{H2}=-4
v(6)=Re{H3}=0
v(7)=Re{H3}=-8
Thus, we have
H0=12+16i;
H1=-8+0i
H2=-4-4i
H3=0-8i
From equation (12.3.5) (c.f. p517 2nd edition), we have
F1=0.5*(H1+H3')-0.5i*(H1-H3')*exp(2i*pi*1/8)
=-9.6569 + 4.0000i
F2=0.5*(H2+H2')-0.5i*(H2-H2')*exp(2i*pi*2/8)
= -4.0000 - 4.0000i
F3=0.5*(H3+H1')-0.5i*(H3-H1')*exp(2i*pi*3/8)
=1.6569 - 4.0000i
F0 and F4 are easily covered by
F0=Re{H0}+Im{H0}
=28
F4=Re{H0}-Im{H0}
=-4
For easily comparing, after the above process, we have (*1):
F0=28
F1=-9.6569 + 4.0000i
F2=-4.0000 - 4.0000i
F3=1.6569 - 4.0000i
F4=-4
Or the output sequence now:
{F0} {F4} {Re{F1}} {Im{F1}} {Re{F2}} {Im{F2}} {Re{F3}} {Im{F3}}
The output of realfft() results the same F as above.

However, if we apply the above real sequence into Matlab
b2=fft(v)
Columns 1 through 6
28.0000 -4.0000 + 9.6569i -4.0000 + 4.0000i -4.0000 + 1.6569i -4.0000 -4.0000 - 1.6569i

Columns 7 through 8

-4.0000 - 4.0000i -4.0000 - 9.6569i

From F (*1) and b2, I couldnot figure out F in the positive or negative range of b2. Of course the F0 and F4 are OK.

Can anyone help me out this issue, I am confused that I couldnot inteprete the theory correctly.

Hereafter is the piece of my code reproduced from p. 518:
int i,i1,i2,i3,i4;
double c1=0.5,c2,h1r,h1i,h2r,h2i,wr,wi,wpr,wpi,wtemp,thet a;
int n=v.length();
theta=PI/(double)(n>>1); // Initialize the recurrence (2*pi/N)
if(isign==1)
{
c2=-0.5;
four1(v,1);// the forward transform is here
}
else
{
c2=0.5;
theta=-theta; // Otherwise, setup for an inverse transform
}// if(isign==1)
wtemp=std::sin(0.5*theta);
wpr=-2.0*wtemp*wtemp;
wpi=std::sin(theta);
wr=1.0+wpr;
wi=wpi;
for(i=1;i<(n>>2);i++)// Case i=0 will be done separately
{
i2=1+(i1=i+i);
i4=1+(i3=n-i1);
// The two separate transform are separated out of data
h1r=c1*(v(i1)+v(i3));
h1i=c1*(v(i2)-v(i4));
h2r=-c2*(v(i2)+v(i4));
h2i=c2*(v(i1)-v(i3));
// Here they are recombined to form the true form
// of the original real data
v(i1)=h1r+wr*h2r-wi*h2i;
v(i2)=h1i+wr*h2i+wi*h2r;
v(i3)=h1r-wr*h2r+wi*h2i;
v(i4)=-h1i+wr*h2i+wi*h2r;
// The recurrence, compute sin/cos(2i*pi*n/N) for next
wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;
}// for(i=1;i<(n>>2);i++)
// Squeeze the first and last data together to get them all
// within the original array
if(isign==1)
{
v(0)=(h1r=v(0))+v(1); // F0
v(1)=h1r-v(1);// FN/2
}
else
{
v(0)=c1*((h1r=v(0))+v(1));
v(1)=c1*(h1r-v(1));
nrfft(v,-1); // This is the inverse transform for the case isign=-1
} // if(isign==1)