C: Numerical Recipies (FFT)

Question

This question is directed at any fans of Numerical Recipes or anyone that understands FFT well.

Can anyone explain why the real component is calculated by -2*(sin(theta/2))^2 ? I can't seem to wrap my head around it. I've seen other examples such as http://www.dspdimension.com/admin/dft-a-pied/ tutorial which simply takes cos(theta) as real and -sin(theta) as imaginary. I've also seen here in basic http://www.dspguide.com/ch12/3.htm which lists it as cos(theta) as real and -sin(theta) as imaginary. I can think of a few more resources that simply take the cos and -sin as real and imaginary.

cos(theta) = 1-2*(sin(theta/2))^2

if the above trig identity is true, then why does this not folllow?

theta=isign*(6.28318530717959/mmax);
wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);

I am assuming Numerical Recipe must be using some trig identity? I can't seem to figure it out and the book doesn't explain at all.

Code found here: http://ronispc.chem.mcgill.ca/ronis/chem593/sinfft.c.html

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(double *data,unsigned long nn,int isign)
{
        unsigned long n,mmax,m,j,istep,i;
        double wtemp,wr,wpr,wpi,wi,theta;
        double tempr,tempi;

        n=nn << 1;
        j=1;
        for (i=1;i<n;i+=2) {
                if (j > i) {
                        SWAP(data[j],data[i]);
                        SWAP(data[j+1],data[i+1]);
                }
                m=n >> 1;
                while (m >= 2 && j > m) {
                        j -= m;
                        m >>= 1;
                }
                j += m;
        }
        mmax=2;
        while (n > mmax) {
                istep=mmax << 1;
                theta=isign*(6.28318530717959/mmax);
                wtemp=sin(0.5*theta);
                wpr = -2.0*wtemp*wtemp;
                wpi=sin(theta);
                wr=1.0;
                wi=0.0;
                for (m=1;m<mmax;m+=2) {
                        for (i=m;i<=n;i+=istep) {
                                j=i+mmax;
                                tempr=wr*data[j]-wi*data[j+1];
                                tempi=wr*data[j+1]+wi*data[j];
                                data[j]=data[i]-tempr;
                                data[j+1]=data[i+1]-tempi;
                                data[i] += tempr;
                                data[i+1] += tempi;
                        }
                        wr=(wtemp=wr)*wpr-wi*wpi+wr;
                        wi=wi*wpr+wtemp*wpi+wi;
                }
                mmax=istep;
        }
}
#undef SWAP

Answer 1

One form of the half angle identity for cosines is:

cos(theta) = 1 - 2*(sin(theta/2)^2)

Not sure if that answers your question.

Answer 2

I don't know FFT well I've done one but its been a long time.

So I looked up trig identities at http://www.sosmath.com/trig/Trig5/trig5/trig5.html

and from the first "Product-To-Sum" identity for sin(u)*sin(v) we have

sin^2(theta/2) = (1/2) [cos(zero) - cos(theta)] = 0.5 - 0.5 cos(theta)

Does this help?

Answer 3

They are using trig identities to minimize the number of circular functions they need to compute. Many fast implementations just look up these functions. If you really want to know you need to just work out the details by unrolling the loop above and making the appropriate variable subsitutions....tedious yes.

BTW, the NR implementation is known to be very sloooow.

Paul

Answer 4

Ok so here is the trig identity. The reason its not the half cos(theta) trig identity is because of the reliance on euler and imaginary numbers. The math is still beyond me...