1 #define MYRECIPLN2 1.442695040888963407359924681001892137426 // 1.0/log(2) 2 3 /* some useful conversions between a number and its power of 2 */ 4 #define LOG2(a) (MYRECIPLN2*log(a)) // floating point logarithm base 2 5 #define POW2(m) ((unsigned long) 1 << (m)) // integer power of 2 for m<32 6 7 /******************************************************************* 8 lower level fft stuff called by routines in fftext.c and fft2d.c 9 *******************************************************************/ 10 11 void fftCosInit(long M, float *Utbl); 12 /* Compute Utbl, the cosine table for ffts */ 13 /* of size (pow(2,M)/4 +1) */ 14 /* INPUTS */ 15 /* M = log2 of fft size */ 16 /* OUTPUTS */ 17 /* *Utbl = cosine table */ 18 19 void fftBRInit(long M, short *BRLow); 20 /* Compute BRLow, the bit reversed table for ffts */ 21 /* of size pow(2,M/2 -1) */ 22 /* INPUTS */ 23 /* M = log2 of fft size */ 24 /* OUTPUTS */ 25 /* *BRLow = bit reversed counter table */ 26 27 void ffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow); 28 /* Compute in-place complex fft on the rows of the input array */ 29 /* INPUTS */ 30 /* *ioptr = input data array */ 31 /* M = log2 of fft size (ex M=10 for 1024 point fft) */ 32 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */ 33 /* *Utbl = cosine table */ 34 /* *BRLow = bit reversed counter table */ 35 /* OUTPUTS */ 36 /* *ioptr = output data array */ 37 38 void iffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow); 39 /* Compute in-place inverse complex fft on the rows of the input array */ 40 /* INPUTS */ 41 /* *ioptr = input data array */ 42 /* M = log2 of fft size */ 43 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */ 44 /* *Utbl = cosine table */ 45 /* *BRLow = bit reversed counter table */ 46 /* OUTPUTS */ 47 /* *ioptr = output data array */ 48 49 void rffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow); 50 /* Compute in-place real fft on the rows of the input array */ 51 /* The result is the complex spectra of the positive frequencies */ 52 /* except the location for the first complex number contains the real */ 53 /* values for DC and Nyquest */ 54 /* INPUTS */ 55 /* *ioptr = real input data array */ 56 /* M = log2 of fft size */ 57 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */ 58 /* *Utbl = cosine table */ 59 /* *BRLow = bit reversed counter table */ 60 /* OUTPUTS */ 61 /* *ioptr = output data array in the following order */ 62 /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */ 63 64 65 void riffts1(float *ioptr, long M, long Rows, float *Utbl, short *BRLow); 66 /* Compute in-place real ifft on the rows of the input array */ 67 /* data order as from rffts1 */ 68 /* INPUTS */ 69 /* *ioptr = input data array in the following order */ 70 /* M = log2 of fft size */ 71 /* Re(x[0]), Re(x[N/2]), Re(x[1]), Im(x[1]), Re(x[2]), Im(x[2]), ... Re(x[N/2-1]), Im(x[N/2-1]). */ 72 /* Rows = number of rows in ioptr array (use Rows of 1 if ioptr is a 1 dimensional array) */ 73 /* *Utbl = cosine table */ 74 /* *BRLow = bit reversed counter table */ 75 /* OUTPUTS */ 76 /* *ioptr = real output data array */ 77