524 Chapter 12. Fast Fourier Transform int idim; unsigned long i1,12,13,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2; unsigned long ibit,k1,k2,n,nprev,nrem,ntot; float tempi,tempr; double theta,wi,wpi,wpr,wr,wtemp; Double precision for trigonometric recur- rences for (ntot=1,idim=1;idim<=ndim;idim++) Compute total number of complex val- ntot *nn[idim]; ues. nprev=1; for (idim=ndim;idim>=1;idim--){ Main loop over the dimensions. n=nn[idim]; nrem=ntot/(n*nprev); ipi=nprev <<1; ip2=1p1*n ip3=ip2*nrem; 188891992 12rev=1; for(12=1;12<=1p2;12+=1p1){ This is the bit-reversal section of the 1.800 (including this one) if (12 i2rev){ routine. for(11=12;11<=12+1p1-2;11+=2)[ 87 for(13=i1;13<=ip3;i3+=ip2) 13rev=12rev+i3-12; 7423 to to tusers to make one paper by Cambridge University Press.Programs from NUMERICAL RECIPES IN SWAP(data[i3],data[i3rev]); SWAP(data[i3+1],data[i3rev+1]); (North America THE ibit-ip2 >1; while (ibit >ip1 &i2rev ibit){ i2rev -ibit; 是 1b1t>>=1; send i2rev +ibit; ifpl=ipl; Here begins the Danielson-Lanczos sec- Copyright (C) while (ifp1 ip2){ tion of the routine. to dir ifp2=ifpl <1; theta=1s1g即*6.28318530717959/(1fp2/1p1); Initialize for the trig.recur- wtemp=sin(0.5*theta); rence. wpr =-2.0*wtemp*wtemp; ectcustser ART OF SCIENTIFIC COMPUTING(ISBN 0-521 wpi=sin(theta); r=1.0; w1=0.0; for(13=1;13<=1fp1;13+=1p1){ for(11=13;11<=13+ip1-2;11+=2){ for(12=i1;i2<=ip3;i2+=1fp2) Numerical Recipes books or 1988-1992 by Numerical Recipes k1=12: Danielson-Lanczos formula: -431085 k2=k1+1fp1; tempr=(float)wr*data[k2]-(float)wi*data[k2+1]; tempi=(float)wr*data[k2+1]+(float)wi*data[k2]; data [k2]=data[ki]-tempr; data[k2+1]=data[k1+1]-tempi; data[k1]+tempr; (outside North America) Software. data[k1+1]+tempi; wr=(wtemp=wr)*wpr-wi*wpi+wr; Trigonometric recurrence wi=wi*wpr+wtemp*wpi+wi; ifpi=ifp2; nprev *n;524 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). { int idim; unsigned long i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2; unsigned long ibit,k1,k2,n,nprev,nrem,ntot; float tempi,tempr; double theta,wi,wpi,wpr,wr,wtemp; Double precision for trigonometric recur￾rences. for (ntot=1,idim=1;idim<=ndim;idim++) Compute total number of complex val￾ntot *= nn[idim]; ues. nprev=1; for (idim=ndim;idim>=1;idim--) { Main loop over the dimensions. n=nn[idim]; nrem=ntot/(n*nprev); ip1=nprev << 1; ip2=ip1*n; ip3=ip2*nrem; i2rev=1; for (i2=1;i2<=ip2;i2+=ip1) { This is the bit-reversal section of the if (i2 < i2rev) { routine. for (i1=i2;i1<=i2+ip1-2;i1+=2) { for (i3=i1;i3<=ip3;i3+=ip2) { i3rev=i2rev+i3-i2; SWAP(data[i3],data[i3rev]); SWAP(data[i3+1],data[i3rev+1]); } } } ibit=ip2 >> 1; while (ibit >= ip1 && i2rev > ibit) { i2rev -= ibit; ibit >>= 1; } i2rev += ibit; } ifp1=ip1; Here begins the Danielson-Lanczos sec￾while (ifp1 < ip2) { tion of the routine. ifp2=ifp1 << 1; theta=isign*6.28318530717959/(ifp2/ip1); Initialize for the trig. recur￾wtemp=sin(0.5*theta); rence. wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; for (i3=1;i3<=ifp1;i3+=ip1) { for (i1=i3;i1<=i3+ip1-2;i1+=2) { for (i2=i1;i2<=ip3;i2+=ifp2) { k1=i2; Danielson-Lanczos formula: k2=k1+ifp1; tempr=(float)wr*data[k2]-(float)wi*data[k2+1]; tempi=(float)wr*data[k2+1]+(float)wi*data[k2]; data[k2]=data[k1]-tempr; data[k2+1]=data[k1+1]-tempi; data[k1] += tempr; data[k1+1] += tempi; } } wr=(wtemp=wr)*wpr-wi*wpi+wr; Trigonometric recurrence. wi=wi*wpr+wtemp*wpi+wi; } ifp1=ifp2; } nprev *= n; } }