40 likes | 75 Views
Fast Fourier Transform. Sung-Ju Kang Department of Physics Kangwon National University. Fast Fourier transform is the algorism for faster compute than discrete Fourier transform. The discrete Fourier transform appears to be an process. These appearances are deceiving!
E N D
Fast Fourier Transform. Sung-Ju Kang Department of Physics Kangwon National University Fast Fourier transform is the algorism for faster compute than discrete Fourier transform. The discrete Fourier transform appears to be an process. These appearances are deceiving! The discrete Fourier transform can, in fact, be computed in operations with an algorithm of the fast Fourier transform, or FFT. The difference between and is immense.
That is involved computation of the discrete Fourier transform of N points. The standard answer was this : Define W as the complex number Then above the numerical formula can be written as One of the two is formed from the even-numbered points of the original N, the odd-numbered points. The proof is simply this : In the last line,W is the same constant as in (1).
Programming with C language #include <stdio.h> #include <math.h> #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr void four1(float data[], unsigned long nn, int isign); void main() { float data[2050],r,p; double x,y,z,w,t; unsigned long nn; int isign,i; char FINN[20],FOUTN[20]; FILE *FIN,*FOUT; printf("\nINPUT THE NAME OF THE INPUT FILE\n"); scanf("%s",FINN); if((FIN=fopen(FINN,"r"))==NULL){ printf("FILE OPEN ERROR...\n"); exit(-1);} printf("\nINPUT THE NAME OF THE OUTPUT FILE\n"); scanf("%s",FOUTN); if((FOUT=fopen(FOUTN,"w"))==NULL){ printf("FILE OPEN ERROR...\n"); exit(-1);} nn=1024; /* No. of data */ isign=1; for (i=1;i<=2*nn;i+=2) { fscanf(FIN,"%lf %lf %lf %lf",&z,&y,&x,&w); data[i]=(float)x;data[i+1]=0.0; } four1(data,nn,isign); for (i=1;i<=2*nn;i+=2) { r=sqrt(pow(data[i],2.0)+pow(data[i+1],2.0)); p=(pow(r,2.0)/(2.0*nn)); t=i*M_PI/nn; fprintf(FOUT,"%14.8le %14.8e\n",t,p); printf("%14.8e %14.8e\n",t,p); } fclose(FIN);fclose(FOUT); }
void four1(float data[], unsigned long nn, int isign) { unsigned long n,mmax,m,j,istep,i; double wtemp,wr,wpr,wpi,wi,theta; float 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 /* (C) Copr. 1986-92 Numerical Recipes Software 0>,m2. */