508 Chapter 12.Fast Fourier Transform ①real ①real 1=0 r=0 imag ② imag real real 1 1=△ f= (④imag imag N△ Permission is read able files -D real f= N2-1 http://www.nr.com or imag N△ +D real (combination) N+2 imag ∫=±2 11-800 (including this one) granted fori N+3 real f=- N/2-1 N+imag N△ from NUMERICAL RECIPES IN 2N-3 real 1=(N-2)△ 2N-2 imag (North America to any server computer, tusers to make one paper 1988-1992 by Cambridge University Press. THE 2w-①real 2N-D real t=(N-1)△ 1 N ④ imag N△ 是 Programs (a) (b) send Figure 12.2.2.Input and output arrays for FFT.(a)The input array contains N (a power of 2) strictly prohibited complex time samples in a real array of length 2N,with real and imaginary parts altemating.(b)The output array contains the complex Fourier spectrum at N values of frequency.Real and imaginary parts to dir Copyright (C) again alternate.The array starts with zero frequency,works up to the most positive frequency (which is ambiguous with the most negative frequency).Negative frequencies follow,from the second-most negative up to the frequency just below zero. ectcustser ART OF SCIENTIFIC COMPUTING(ISBN wpi=sin(theta); wr=1.0; v@cam w1=0.0; 10-621 for (m=1;m<mmax;m+=2){ Here are the two nested inner loops. for (i=m;i<=n;i+=istep) 1988-1992 by Numerical Recipes j=i+mmax; This is the Danielson-Lanczos for- 43108-5 tempr-wr*data[j]-wi*data[j+1]; mula: tempi=wr*data[j+1]+wi*data[j]; data[j]=data [i]-tempr; (outside data[j+1]=data[i+1]-tempi; data[i]tempr; North Software. data[i+1]+tempi; wr=(wtemp=wr)*wpr-wi*wpi+wr; Trigonometric recurrence. America) wi=wi*wpr+wtemp*wpi+wi; mmax=istep; (A double precision version of four1,named dfour1,is used by the routine mpmul in $20.6.You can easily make the conversion,or else get the converted routine from the Numerical Recipes diskette.508 Chapter 12. Fast Fourier Transform Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 1 2 3 4 real imag real imag t = 0 t = ∆ real imag real imag t = (N − 2)∆ t = (N − 1)∆ real array of length 2 N 1 2 3 4 real imag real imag f = 0 f = N − 1 N N + 1 N + 2 N + 3 N + 4 real imag real imag real imag f = f = ± (combination) f = − real array of length 2 N 1 N∆ N/2 − 1 N∆ 1 2∆ 2N − 1 2N real imag f = − 1 N∆ 2N − 3 2N − 2 2N − 1 2N N/2 − 1 N∆ (a) (b) Figure 12.2.2. Input and output arrays for FFT. (a) The input array contains N (a power of 2) complex time samples in a real array of length 2N, with real and imaginary parts alternating. (b) The output array contains the complex Fourier spectrum at N values of frequency. Real and imaginary parts again alternate. The array starts with zero frequency, works up to the most positive frequency (which is ambiguous with the most negative frequency). Negative frequencies follow, from the second-most negative up to the frequency just below zero. wpi=sin(theta); wr=1.0; wi=0.0; for (m=1;m<mmax;m+=2) { Here are the two nested inner loops. for (i=m;i<=n;i+=istep) { j=i+mmax; This is the Danielson-Lanczos for￾tempr=wr*data[j]-wi*data[j+1]; mula: 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; Trigonometric recurrence. wi=wi*wpr+wtemp*wpi+wi; } mmax=istep; } } (A double precision version of four1, named dfour1, is used by the routine mpmul in §20.6. You can easily make the conversion, or else get the converted routine from the Numerical Recipes diskette.)