FFT的C语言实现
void fft(COMPLEX *x, int m)
{
? COMPLEX *w;?????????? /* used to store the w complex array */
? int mstore = 0;?????? /* stores m for future reference */
? int n = 1;??????????? /* length of fft stored for future */
?COMPLEX u,temp,tm;
?COMPLEX *xi,*xip,*xj,*wptr;
?int i,j,k,l,le,windex;
?double arg,w_real,w_imag,wrecur_real,wrecur_imag,wtemp_real;
?if(m != mstore) {
??/* free previously allocated storage and set new m */
??if(mstore != 0) free(w);
??mstore = m;
??if(m == 0) return;?????? /* if m=0 then done */
??/* n = 2**m = fft length */
??n = 1 << m;
??le = n/2;
??/* allocate the storage for w */
??w = (COMPLEX *) calloc(le-1,sizeof(COMPLEX));
??if(!w) {
???printf("\nUnable to allocate complex W array\n");
???exit(1);
??}
??/* calculate the w values recursively */
??arg = 4.0*atan(1.0)/le;???????? /* PI/le calculation */
??wrecur_real = w_real = cos(arg);
??wrecur_imag = w_imag = -sin(arg);
??xj = w;
??for (j = 1 ; j < le ; j++) {
???xj->real = (float)wrecur_real;
???xj->imag = (float)wrecur_imag;
???xj++;
???wtemp_real = wrecur_real*w_real - wrecur_imag*w_imag;
???wrecur_imag = wrecur_real*w_imag + wrecur_imag*w_real;
???wrecur_real = wtemp_real;
??}
?}
?/* start fft */
?le = n;
?windex = 1;
?for (l = 0 ; l < m ; l++) {
??le = le/2;
??/* first iteration with no multiplies */
??for(i = 0 ; i < n ; i = i + 2*le) {
???xi = x + i;
???xip = xi + le;
???temp.real = xi->real + xip->real;
???temp.imag = xi->imag + xip->imag;
???xip->real = xi->real - xip->real;
???xip->imag = xi->imag - xip->imag;
???*xi = temp;
??}
??/* remaining iterations use stored w */
??wptr = w + windex - 1;
??for (j = 1 ; j < le ; j++) {
???u = *wptr;
???for (i = j ; i < n ; i = i + 2*le) {
????xi = x + i;
????xip = xi + le;
????temp.real = xi->real + xip->real;
????temp.imag = xi->imag + xip->imag;
????tm.real = xi->real - xip->real;
????tm.imag = xi->imag - xip->imag;
????xip->real = tm.real*u.real - tm.imag*u.imag;
????xip->imag = tm.real*u.imag + tm.imag*u.real;
????*xi = temp;
???}
???wptr = wptr + windex;
??}
??windex = 2*windex;
?}
?/* rearrange data by bit reversing */
?j = 0;
?for (i = 1 ; i < (n-1) ; i++) {
??k = n/2;
??while(k <= j) {
???j = j - k;
???k = k/2;
??}
??j = j + k;
??if (i < j) {
???xi = x + i;
???xj = x + j;
???temp = *xj;
???*xj = *xi;
???*xi = temp;
??}
?}
?
?free(w);
?w=NULL;
}
