Fast Fourier Transform Source File
Submitted by Daryle on Fri, 05/18/2012 - 13:37
Filename: | fft.c |
Project Name: | Fast Fourier Transforms (FFT) |
Language: | C |
Description: | |
Code:
/******************************************** ** FFT: the subroutine to calculate and ** ** return the 1D-FF transform of a input- ** ** ted image. ** ** ** ** reverse: a subroutine to take an ** ** integer of n bits and reverse it from ** ** high-endian to low-endian order. ** ** ** ** Author: Daryle Niedermayer ** ** Date: Oct 5, 1999 ** ** Credits: Thanks to Paul Bourke for ** ** his contribution and explan- ** ** ations ** ** (http://www.swin.edu.au/ ** ** austronomy/pbourke) ** ********************************************/ /******************************************** ** REVERSE: Strategy... ** ** using a single pass, take the most sig- ** ** nificant and least significant bits and ** ** switch them. For this purpose, leftbit ** ** and rightbit are temporary values of ** ** the bits requiring switching, and left- ** ** place and rightplace are the positions ** ** of these bits in the integer. ** ********************************************/ unsigned char reverse(unsigned char input, int n) { unsigned char leftbit, rightbit, output; int i, leftplace, rightplace; leftplace = n-1, rightplace = 0; output = 0; while (leftplace >= rightplace) { /* Create bit masks for left and right places */ leftbit = 1 << leftplace; rightbit = 1 << rightplace; /* Use masks to get the value of the bits */ leftbit = input & leftbit; rightbit = input & rightbit; /* Switch the places for these bits */ leftbit = leftbit >> (leftplace - rightplace); rightbit = rightbit << (leftplace - rightplace); /* Use masks to insert values into input */ output = output | leftbit; output = output | rightbit; /* Increment/Decrement placement holders */ leftplace--; rightplace++; } /* while (leftplace > rightplace) */ return output; } /******************************************** ** DFT calculates the 1D-Discrete Fourier ** ** Transform using the definition of a ** ** Fourier Transform: ** ** F(u)=(1/N)Sum[from x=0 to N-1]f(x)exp^ ** ** (-j2*pi*ux/N) ** ********************************************/ COMPLEX *DFT(COMPLEX F[], int N) { int u, x; double cosarg = 0; double sinarg = 0; double W; COMPLEX *ft; COMPLEX tmp; for (u=0; u<N; u++) { /* Zero out the temporary array as we go */ TF[u].r = 0; TF[u].i = 0; W = 2*PI*u/N; for (x=0; x<N; x++) { tmp.r = cos(W*x); tmp.i = -1*sin(W*x); /**************************************************** ** To find the sums in complex form, we must ** cross-multiply f(x) with the cos and sin functions ** from Euler's Formula: ** TF.r*(cosarg - sinarg) + TF.i*(cosarg - sinarg) = ** TF.r*cosarg - TF.i*(sinarg) = TF.r ** TF.i*cosarg + TF.r*(sinarg) = TF.i ****************************************************/ TF[u] = complex_plus(TF[u],(complex_times(F[x],tmp))); } TF[u].r = TF[u].r/N; TF[u].i = TF[u].i/N; } ft = TF; return ft; } /******************************************** ** FFT calculates the 1D-FFT of a complex ** ** array passed into the function. The ar- ** ** ray is expected to be of size N. ** ********************************************/ COMPLEX *FFT(COMPLEX F[], int N) { int i, M, j, k, u, i1, i2; COMPLEX F2[N]; COMPLEX Half, W2M; int new_index[N]; int n=(int) log(N) / log(2); #ifdef _DEBUG_ printf("The value of n is: %d\n",n); #endif Half.r = 0.5; Half.i = 0.0; /***************************************** ** Reverse the array indexes *****************************************/ for (i=0; i<N; i++) { new_index[i] = reverse(i,n); } /***************************************** ** Rearrange array elements *****************************************/ for (i=0; i<N; i++) { TF[new_index[i]] = F[i]; } /***************************************** ** Begin decomposing the array into ** subgroups. *****************************************/ M = 1; /* The initial length of subgroups */ j = N/2; /* The number of pairs of subgroups */ /** Begin successive merges for n levels */ for (i=0; i<n; i++) { /** Merge pairs at the current level */ for (k=0; k<j; k++) { i1 = k*2*M; /* Start of first group */ i2 = (k*2+1)*M; /* Start of second group */ #ifdef _DEBUG_ printf("Cycling through k=%d... i1 and i2 = %d and %d...\n",k,i1, i2); #endif for (u=0; u<M; u++) { W2M.r = cos(PI*u/M); W2M.i = -sin(PI*u/M); #ifdef _DEBUG_ printf("M is currently: %d\n",M); printf("Cycling through u=%d...\n",u); printf("Cycling through i=%d... W2M=%f\n",i,W2M); #endif /* Calculate the values for the first set of numbers */ F2[u] = complex_times(Half,(complex_plus(TF[i1+u], \ complex_times(TF[i2+u],W2M)))); /*Calculate the values for the second set of numbers */ F2[u+M] = complex_times(Half,(complex_minus(TF[i1+u], \ complex_times(TF[i2+u],W2M)))); TF[i1+u] = F2[u]; TF[i1+u+M] = F2[u+M]; #ifdef _DEBUG_ printf("F2[%d] is: ",u); complex_print(F2[u]); printf("\n"); printf("F2[%d] is: ",u+M); complex_print(F2[u+M]); printf("\n"); printf("TF[%d] is: ",i1+u); complex_print(TF[i1+u]); printf("\n"); printf("TF[%d] is: ",i1+u+M); complex_print(TF[i1+u+M]); printf("\n"); #endif } } M = 2*M; j = j/2; } return TF; } /******************************************** ** bit_print is a little utility to print ** ** out an integer bitwise for debugging ** ** purposes. ** ********************************************/ void bit_print(int a) { int i; int n = sizeof(int) * CHAR_BIT; int mask = 1 << (n - 1); for (i=1; i<=n; ++i) { putchar(((a & mask) == 0) ? '0' : '1'); a <<= 1; if (i % CHAR_BIT == 0 && i<n) putchar(' '); } printf("\n"); }