Listing One /* PCFFT.C -- by J.G.G. Dobbe -- Performs an FFT on two arrays (Re, Im) of type float (can be changed). This unit is written in C and doesn't call assembler routines. */ /* --------------------- Include directive ------------------------ */ #include "pcfft.h" /* --------------------- Local variables -------------------------- */ static float CosArray[28] = { /* cos{-2pi/N} for N = 2, 4, 8, ... 16384 */ -1.00000000000000, 0.00000000000000, 0.70710678118655, 0.92387953251129, 0.98078528040323, 0.99518472667220, 0.99879545620517, 0.99969881869620, 0.99992470183914, 0.99998117528260, 0.99999529380958, 0.99999882345170, 0.99999970586288, 0.99999992646572, /* cos{2pi/N} for N = 2, 4, 8, ... 16384 */ -1.00000000000000, 0.00000000000000, 0.70710678118655, 0.92387953251129, 0.98078528040323, 0.99518472667220, 0.99879545620517, 0.99969881869620, 0.99992470183914, 0.99998117528260, 0.99999529380958, 0.99999882345170, 0.99999970586288, 0.99999992646572 }; static float SinArray[28] = { /* sin{-2pi/N} for N = 2, 4, 8, ... 16384 */ 0.00000000000000, -1.00000000000000, -0.70710678118655, -0.38268343236509, -0.19509032201613, -0.09801714032956, -0.04906767432742, -0.02454122852291, -0.01227153828572, -0.00613588464915, -0.00306795676297, -0.00153398018628, -0.00076699031874, -0.00038349518757, /* sin{2pi/N} for N = 2, 4, 8, ... 16384 */ 0.00000000000000, 1.00000000000000, 0.70710678118655, 0.38268343236509, 0.19509032201613, 0.09801714032956, 0.04906767432742, 0.02454122852291, 0.01227153828572, 0.00613588464915, 0.00306795676297, 0.00153398018628, 0.00076699031874, 0.00038349518757 }; /* --------------------- Function implementations ----------------- */ /* --------------------- ShuffleIndex ----------------------------- */ static unsigned int ShuffleIndex(unsigned int i, int WordLength) /* Function : Finds the shuffle index of array elements. The array length must be a power of two; The power is stored in "WordLength". Return value : With "i" the source array index, "ShuffleIndex" returns the destination index for shuffling. Comment : - */ { unsigned int NewIndex; unsigned char BitNr; NewIndex = 0; for (BitNr = 0; BitNr <= WordLength - 1; BitNr++) { NewIndex = NewIndex << 1; if ((i & 1) != 0) NewIndex = NewIndex + 1; i = i >> 1; } return NewIndex; } /* --------------------- Shuffle2Arr ------------------------------ */ static void Shuffle2Arr(float *a, float *b, int bitlength) /* Function : Shuffles both arrays "a" and "b". This function is called before performing the actual FFT so the array elements are in the right order after FFT. Return value : - Comment : - */ { unsigned int IndexOld, IndexNew; float temp; unsigned int N; int bitlengthtemp; bitlengthtemp = bitlength; /* Save for later use */ N = 1; /* Find array-length */ do { N = N * 2; bitlength = bitlength - 1; } while (bitlength > 0) ; /* Shuffle all elements */ for (IndexOld = 0; IndexOld <= N - 1; IndexOld++) { /* Find index to exchange elements */ IndexNew = ShuffleIndex(IndexOld, bitlengthtemp); if (IndexNew > IndexOld) { /* Exchange elements: */ temp = a[IndexOld]; /* Of array a */ a[IndexOld] = a[IndexNew]; a[IndexNew] = temp; temp = b[IndexOld]; /* Of array a */ b[IndexOld] = b[IndexNew]; b[IndexNew] = temp; } } } /* --------------------- Fft -------------------------------------- */ void Fft(float *Re, float *Im, int Pwr, int Dir) /* Function : Actual FFT algorithm. "Re" and "Im" point to start of real and imaginary arrays of numbers, "Pwr" holds the array sizes as a power of 2 while "Dir" indicates whether an FFT (Dir>=1) or an inverse FFT must be performed (Dir<=0). Return value : The transformed information is returned by "Re" and "Im" (real and imaginary part respectively). Comment : - */ { int pwrhelp; int N; int Section; int AngleCounter; int FlyDistance; int FlyCount; int index1; int index2; float tempr, tempi; float Re1, Re2, Im1, Im2; float c, s; float scale; float sqrtn; float temp; float Qr, Qi; Shuffle2Arr(Re, Im, Pwr); /* Shuffle before (i)FFT */ pwrhelp = Pwr; /* Determine size of arrs */ N = 1; do { N = N * 2; pwrhelp--; } while (pwrhelp > 0) ; if (Dir >= 1) AngleCounter = 0; /* FFT */ else AngleCounter = 14; /* Inverse FFT */ Section = 1; while (Section < N) { FlyDistance = 2 * Section; c = CosArray[AngleCounter]; s = SinArray[AngleCounter]; Qr = 1; Qi = 0; for (FlyCount = 0; FlyCount <= Section - 1; FlyCount++) { index1 = FlyCount; do { index2 = index1 + Section; /* Perform 2-Point DFT */ tempr = 1.0 * Qr * Re[index2] - 1.0 * Qi * Im[index2]; tempi = 1.0 * Qr * Im[index2] + 1.0 * Qi * Re[index2]; Re[index2] = Re[index1] - tempr; /* For Re-part */ Re[index1] = Re[index1] + tempr; Im[index2] = Im[index1] - tempi; /* For Im-part */ Im[index1] = Im[index1] + tempi; index1 = index1 + FlyDistance; } while (index1 <= (N - 1)); /* k */ /* Calculate new Q = cos(ak) + j*sin(ak) = Qr + j*Qi */ /* -2*pi */ /* with: a = ----- */ /* N */ temp = Qr; Qr = Qr*c - Qi*s; Qi = Qi*c + temp*s; } Section = Section * 2; AngleCounter = AngleCounter + 1; } if (Dir <= 0) /* Normalize for */ { /* inverse FFT only */ scale = 1.0/N; for (index1 = 0; index1 <= N - 1; index1++) { Re[index1] = scale * Re[index1]; Im[index1] = scale * Im[index1]; } } } /* ---------------------------------------------------------------- */ Listing Two /* PCFFT.H -- by J.G.G. Dobbe -- Headers for PCFFT.C */ #ifndef PCFFT_H /* If not defined yet, use this file */ #define PCFFT_H /* --------------------- External function ------------------------ */ void Fft(float *Re, float *Im, int Pwr, int Dir); /* ---------------------------------------------------------------- */ #endif Listing Three /* PCFFTest.C -- by J.G.G. Dobbe -- Test program written in Turbo C/Borland C that uses the C or Assembler version of FFT. It depends on the type of FFT object file whether the C- or ASM-version is linked in. In both cases, the same FFT.H file is used. */ /* --------------------- Include directives ----------------------- */ #include #include #include "pcfft.h" /* --------------------- Constant definition ---------------------- */ #define SIZE 16 /* Size of data arrays (re, im) */ /* --------------------- Variable definitions --------------------- */ float re[SIZE]; /* Array holding Real part of data */ float im[SIZE]; /* Array holding Imaginary part of data */ /* --------------------- Function implementations ----------------- */ /* --------------------- DispArr ---------------------------------- */ void DispArr(char *Txt) { int i; /* Loop counter */ clrscr(); /* Clear screen */ printf("\n%s:\n", Txt); /* Display header */ for (i = 0; i < SIZE; i++) /* Display data points */ printf("i = %4d: Re = %8.2f, Im = %8.2f\n", i, re[i], im[i]); printf("Press to continue\n"); /* Display message */ getch(); /* Wait for keystroke */ } /* --------------------- main ------------------------------------- */ int main() { int i; for (i = 0; i < SIZE; i++) /* Clear data arrays */ { re[i] = 0.0; im[i] = 0.0; } re[0] = 100.0; /* Fill array with pulse signal */ DispArr("Input Time Data"); /* Display time data (pulse signal) */ Fft(re, im, 4, 1); /* FFT on data points */ DispArr("Frequency Response (after FFT)"); /* Display freq data */ Fft(re, im, 4, -1); /* Inverse FFT on freq data */ DispArr("Time Response (after inverse FFT)"); /* Display time data (pulse signal) */ return 0; } /* ---------------------------------------------------------------- */