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

Start from:

• cos(A+B) = cos(A) cos(B) - sin(A) sin(B)
• sin(A+B) = sin(A) cos(B) + cos(A) sin(B)
• cos(2A) = 1 - 2 sin²(A)
• e^{i θ} = cos(θ) + i sin(θ)

So:

e^{i (φ+δ)}

= cos(φ + δ) + i sin(φ + δ)

= cos(φ) cos(δ) - sin(φ) sin(δ) + i [sin(φ) cos(δ) + cos(φ) sin(δ)]

= cos(φ) [ 1 - 2 sin²(δ/2) ] + i sin(φ) [ 1 - 2 sin²(δ/2) ] + i sin(δ) [ i * sin(φ) + cos(φ) ]

= [ cos(φ) + i sin(φ) ] [ 1 - 2 sin²(δ/2) ] + [ cos(φ) + i sin(φ) ] i sin(δ)

= e^{i φ} + e^{i φ} [ - 2 sin²(δ/2) + i sin(δ)]

Edit: That was a lot of useless formatting on my part. It's actually way simpler:

y^(a+b) = y^a × y^b for any y. So:

e^{i (φ+δ)}

= e^{i φ} e^{i δ}

= e^{i φ} [ cos(δ) + i sin(δ) ]

= e^{i φ} [ 1 - 2 sin²(δ/2) + i sin(δ) ]

Answer 2

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 3

The reason is for numerical accuracy. If you closely examine the following code:

wtemp=sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi=sin(theta);

and

wr=(wtemp=wr)*wpr-wi*wpi+wr;
wi=wi*wpr+wtemp*wpi+wi;

they are designed to work together to yield the correct expected result. A naive approach would be implemented as follows:

wpr = cos(theta);
wpi = sin(theta);

and

wtemp = wr;
wr =wr*wpr - wi*wpi;
wi =wi*wpr + wtemp*wpi;

and with infinite precision would be functionally equivalent.

However when theta is close to zero (i.e. a lot of sample points or large nn), cos(theta) becomes problematic, because for small angles cos(theta) approaches 1 and has a slope approaching 0. And at some angle cos(theta) == 1. My experimentation shows that for float cos(2*PI/N) == 1 exactly for N >= 25736 for float (i.e. 32-bit precision). A 25,736 point FFT is conceivable. So to avoid this problem wpr is computed as cos(theta)-1 using the half angle formula from trigonometry. It is computed with sin which is very precise for small angles, thus for small angles both wpr and wpi are small and precise. This is then undone in the update code by adding the 1 back back on after the complex multiply. Expressed mathematically, we get:

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

and the update rule is:

w = w * (w_p + 1) 
  = w + w*w_p

Answer 4

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 5

What is confusing is that NR uses the Math/Physics version of the FFT which rotates the twiddle factors the opposite way that EEs define the FFT. So the NR forward is the EE inverse and visa versa. Notice that in NR the forward has a positive exponent instead of the EE negative exponent. The EE method transforms time to frequency where the Math and Physics version plays with angular momentum.

Answer 6

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 substitutions....tedious yes.

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

Paul

Answer 7

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...

_{(source: librow.com)}