534 Chapter 12. Fast Fourier Transform jk=nn [++nv] k8>>=1; if (ks KBF){ for(j12=1;j12<=2;j12++)[ for (kr=1;kr<=ns;kr+=ks/KBF){ for (k=1;k<=ks:k+-KBF){ cc=fread(afa[1],sizeof(float),KBF,file[na]); if (cc !KBF)nrerror("read error in fourfs"); cc=fwrite(&afa[1],sizeof(float),KBF,file[nc]); if (cc !KBF)nrerror("write error in fourfs"); http://www.nr. readable files Permission is nc=mate[nc]; na=mate [na]; -800472 co (including this one) granted fori fourew(file,&na,&nb,&nc,&nd); else if (ks ==KBF)nb=na; else break; j=1; The second phase of the transform starts here.Now,the remaining permutations are suf- 7423 ficiently local to be done in place. from NUMERICAL RECIPES IN C: for (;i){ theta=isign*3.141592653589793/(n/mm); wtemp=sin(0.5*theta); (North America to any server computer,is strictly prohibited. tusers to make one paper 1988-1992 by Cambridge University Press.Programs THE wpr =-2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; w1=0.0; only),or mm>>=1; ks=kd; copy for their kd>>=1: for(j12=1:j12<=2;j12++){ for (kr=1;kr<=ns;kr++) cc=fread(afc[1],sizeof(float),KBF,file[na]); if (cc !KBF)nrerror("read error in fourfs"); kk=1; k-ks+1; for(;;){ ART OF SCIENTIFIC COMPUTING(ISBN 0-521- tempr=((float)wr)*afc[kk+ks]-((float)wi)*afc[kk+ks+1]; tempi=((float)wi)*afc[kk+ks]+((float)wr)*afc[kk+ks+1]; afa[j]=afc[kk]+tempr; afb[i]=afc[kk]-tempr; afa[++i]afc[++kk]+tempi; afb[j++]=afc[kk++]-tempi; To order Numerical Recipes books or personaluse.Further reproduction,or 1988-1992 by Numerical Recipes if (kk k)continue; -431085 kc +kd; if (kc ==mm){ kc=0; wr=(wtemp=wr)*wpr-wi*wpi+wr; Software. wi=wi*wpr+wtemp*wpi+wi; @cambridge.org(outside North America). kk +ks; if (kk KBF)break; else k=kk+ks; if (j>KBF){ cc=fwrite(&afa[1],sizeof(float),KBF,file [nc]); if (cc !KBF)nrerror("write error in fourfs"); cc=fwrite(&afb[1],sizeof(float),KBF,file [nd]); if (cc !KBF)nrerror("write error in fourfs"); j=1;534 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). jk=nn[++nv]; } ks >>= 1; if (ks > KBF) { for (j12=1;j12<=2;j12++) { for (kr=1;kr<=ns;kr+=ks/KBF) { for (k=1;k<=ks;k+=KBF) { cc=fread(&afa[1],sizeof(float),KBF,file[na]); if (cc != KBF) nrerror("read error in fourfs"); cc=fwrite(&afa[1],sizeof(float),KBF,file[nc]); if (cc != KBF) nrerror("write error in fourfs"); } nc=mate[nc]; } na=mate[na]; } fourew(file,&na,&nb,&nc,&nd); } else if (ks == KBF) nb=na; else break; } j=1; The second phase of the transform starts here. Now, the remaining permutations are suf- ficiently local to be done in place. for (;;) { theta=isign*3.141592653589793/(n/mm); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); wr=1.0; wi=0.0; mm >>= 1; ks=kd; kd >>= 1; for (j12=1;j12<=2;j12++) { for (kr=1;kr<=ns;kr++) { cc=fread(&afc[1],sizeof(float),KBF,file[na]); if (cc != KBF) nrerror("read error in fourfs"); kk=1; k=ks+1; for (;;) { tempr=((float)wr)*afc[kk+ks]-((float)wi)*afc[kk+ks+1]; tempi=((float)wi)*afc[kk+ks]+((float)wr)*afc[kk+ks+1]; afa[j]=afc[kk]+tempr; afb[j]=afc[kk]-tempr; afa[++j]=afc[++kk]+tempi; afb[j++]=afc[kk++]-tempi; if (kk < k) continue; kc += kd; if (kc == mm) { kc=0; wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; } kk += ks; if (kk > KBF) break; else k=kk+ks; } if (j > KBF) { cc=fwrite(&afa[1],sizeof(float),KBF,file[nc]); if (cc != KBF) nrerror("write error in fourfs"); cc=fwrite(&afb[1],sizeof(float),KBF,file[nd]); if (cc != KBF) nrerror("write error in fourfs"); j=1; }