00001
00006 #ifndef FFT_DATA_TYPE
00007 #include <stdlib.h>
00008 #include <malloc.h>
00009 #include <stdio.h>
00010 #include <math.h>
00011 #include <time.h>
00012
00013 #include "fft.h"
00014
00015 #ifndef NULL
00016 #define NULL 0
00017 #endif
00018 #ifndef TRUE
00019 #define TRUE 1
00020 #endif
00021 #ifndef FALSE
00022 #define FALSE 0
00023 #endif
00024 #ifndef PI
00025 #ifdef M_PI
00026 #define PI M_PI
00027 #else
00028 #define PI 3.1415926
00029 #endif
00030 #endif
00031
00032 static double *sintab = NULL;
00033 static int *bitrevtab = NULL;
00034 static int lastlog = (-1);
00035 static int lastalloc = 0;
00036
00037
00038
00039
00040 static int
00041 setup_tables(int log2size)
00042 {
00043 int offset, index1, index2;
00044 int buffln = 1 << log2size;
00045 int halfln = buffln >> 1;
00046 int quadln = halfln >> 1;
00047
00048 if (lastalloc < buffln) {
00049 if (sintab)
00050 free(sintab);
00051 if (bitrevtab)
00052 free(bitrevtab);
00053 sintab = (double *)malloc(sizeof(double) * (quadln + 1));
00054 bitrevtab = (int *)malloc(sizeof(int) * buffln);
00055 if (!sintab || !bitrevtab) {
00056 lastalloc = 0;
00057 return (FALSE);
00058 }
00059 lastalloc = buffln;
00060 }
00061 for (index1 = 0; index1 <= quadln; index1++) {
00062 sintab[index1] = sin((double)index1 * 2.0 * PI / (double)buffln);
00063 }
00064 bitrevtab[0] = 0;
00065 offset = halfln;
00066 for (index1 = 1; index1 < buffln; index1 <<= 1) {
00067 for (index2 = 0; index2 < index1; index2++) {
00068 bitrevtab[index2 + index1] = bitrevtab[index2] + offset;
00069 }
00070 offset >>= 1;
00071 }
00072 lastlog = log2size;
00073 return (TRUE);
00074 }
00075
00076 #define FFT_DATA_TYPE float
00077 #define fft fft32
00078 #include "fft.c"
00079 #undef FFT_DATA_TYPE
00080 #undef fft
00081
00082 #define FFT_DATA_TYPE double
00083 #endif
00084
00085 int
00086 fft(FFT_DATA_TYPE re[], FFT_DATA_TYPE im[], int log2size, int inverse)
00087 {
00088 int buffln, halfln, quadln;
00089 int blkcnt, blkidx, sinarg;
00090 int offset, index1, index2;
00091
00092 if (!re || !im)
00093 return (FALSE);
00094 if (log2size <= 2)
00095 return (FALSE);
00096 buffln = 1 << log2size;
00097 halfln = buffln >> 1;
00098 quadln = halfln >> 1;
00099 if ((lastlog != log2size) && !setup_tables(log2size)) {
00100 return (FALSE);
00101 }
00102
00103
00104
00105 index1 = buffln;
00106 do {
00107 index1--;
00108 if (index1 < (index2 = bitrevtab[index1])) {
00109 FFT_DATA_TYPE reswap = re[index1];
00110 FFT_DATA_TYPE imswap = im[index1];
00111 re[index1] = re[index2];
00112 im[index1] = im[index2];
00113 re[index2] = reswap;
00114 im[index2] = imswap;
00115 }
00116 } while (index1 > 0);
00117
00118
00119
00120 offset = 1;
00121 blkcnt = halfln;
00122 do {
00123 sinarg = halfln - blkcnt;
00124 index1 = offset + buffln;
00125 index2 = offset + index1;
00126 offset <<= 1;
00127 do {
00128 double wsin, wcos;
00129 if (sinarg <= quadln) {
00130 wsin = sintab[sinarg];
00131 wcos = sintab[quadln - sinarg];
00132 }
00133 else {
00134 wsin = sintab[halfln - sinarg];
00135 wcos = -sintab[sinarg - quadln];
00136 }
00137 if (!inverse) {
00138 wsin = -wsin;
00139 }
00140 blkidx = blkcnt;
00141 index1 -= (buffln + 1);
00142 index2 -= (buffln + 1);
00143 do {
00144 FFT_DATA_TYPE remult,immult;
00145 remult = (FFT_DATA_TYPE)(wcos * re[index2] - wsin * im[index2]);
00146 immult = (FFT_DATA_TYPE)(wcos * im[index2] + wsin * re[index2]);
00147 re[index2] = re[index1] - remult;
00148 im[index2] = im[index1] - immult;
00149 re[index1] += remult;
00150 im[index1] += immult;
00151 index1 += offset;
00152 index2 += offset;
00153 } while (--blkidx);
00154 } while ((sinarg -= blkcnt) >= 0);
00155 blkcnt >>= 1;
00156 } while (--log2size);
00157 return (TRUE);
00158 }
