C: numerical repeats (FFT)

Question

This question is for any Numerical Recipes fan or anyone who understands FFT well.

Can someone explain why the real component is calculated at -2 * (sin (theta / 2)) ^ 2? I don’t seem to wrap my head around it. I saw other examples, such as a http://www.dspdimension.com/admin/dft-a-pied/ tutorial that simply accepts cos (theta) as real and -sin (theta) as imaginary. I also saw here at the base http://www.dspguide.com/ch12/3.htm , which lists it as cos (theta) as real and -sin (theta) as imaginary. I can come up with a few more resources that simply accept cos and -sin as real and imaginary.

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

if the above trigger value is true, then why doesn't it look like?

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

I assume Numericical Recipe should use some kind of trigger? I cannot understand this, and the book does not 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 with:

• 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 & theta;} = cos (? theta;) + i sin (? theta;)

So:

e ^{i (ph + del)}

= cos (? + delta) + i sin (? + delta)

= cos (?) cos (?) - sin (?) sin (?) + i [sin (?) cos (?) + cos (?) sin ( & delta;)]

= cos (?) [1 - 2 sin ² (? / 2)] + i sin (?) [1 - 2 sin ² (? / 2)] + i sin (? ) [i * sin (? phis;) + cos (? phi;)]

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

= e ^{i & phi;} + e ^{i & phi;} [- 2 sin ² (? / 2) + i sin (?);]

Edit: It was a lot of useless formatting on my part. This is actually simpler:

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

e ^{i (ph + del)}

= e ^{i & phi;} e ^i?

= e ^{i & phi;} [cos (?) + i sin (?)]

= e ^{i & phi;} [1 - 2 sin ² (? / 2) + i sin (?)]

Answer 2

One form of polygon identity for cosines:

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

Not sure if this answers your question.

Answer 3

I do not know the FFT, I did one, but it was a long time.

So, I looked at trigger 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 it help?

Answer 4

They use trigger identities to minimize the number of circular functions that they need to calculate. Many quick implementations just look for these features. If you really want to know that you just need to drill down to the details by expanding the loop above and making the appropriate variable replacements ... tedious yes.

By the way, the implementation of NR is known very slowly.

Floor

Answer 5

So here is the trigger identity. The reason this is not a half cos (theta) trig identity is due to the Euler and imaginary numbers. The math is still after me ...

_{(source: librow.com )}

Answer 6

What is misleading is that NR uses the FFT version of Math / Physics, which rotates the torsion coefficients in the opposite way that the EFTs determine the FFT. Thus, NR forward is the EE reverse and vice versa. Note that in NR, the forward has a positive indicator instead of a negative EE indicator. The EE method converts time to frequency, where the version of Math and Physics is played back with angular momentum.

Answer 7

The reason is numerical accuracy. If you carefully study 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 provide the correct expected result. The naive approach will 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 accuracy will be functionally equivalent.

However, when theta is close to zero (i.e. many sampling points or large nn ), cos(theta) becomes problematic because at small angles cos(theta) approaches 1 and has a slope approaching 0. And for some angle cos(theta) == 1 . My experiments show that for float cos(2*PI/N) == 1 exactly for N >= 25736 for float (i.e. 32-bit precision). Possible 25.736 FFT point. Therefore, in order to avoid this problem, wpr calculated as cos(theta)-1 using the trigonometry polygon formula. It is calculated using sin , which is very accurate for small angles, so for small angles both wpr and wpi are small and accurate. Then this code is canceled in the update code, adding 1 back after the completion of complex multiplication. Mathematically speaking, we get:

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

and update rule:

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