/* Fast Fourier Transform Probably the most basic type. Being unable to comprehend anyone's description of the algorithm for years, it occured to me the other day that it can't be anything other than optimization of the non-fast version. So I wrote out all of the operations, looked for obvious optimizations, and made them. Probably the code could be optimized a lot more, but nothing like the drastic change between the slow fourier transform algorithm which is O(n*n) vs. one like this which is O(n*log(n)), and thus I don't care. Sadly, I don't think this is any easier to understand than anyone else's description of the algorigthm, but if you're trying to understand how it works, I left some printf()s in there you can uncomment to see a list of all of the operations it performs. Probably not much help, though. First call fft_initalize(), first argument is the number of samples, which must be a power of 2, second argument is 0 for a normal transform (to convert sample data to frequency data) or 1 for a reverse transform (to convert frequency data to sample data). Then call fft_fast(), first argument is a complex * for the output data, second argument is a complex * for the input data. If your input data is real data (which is probably the case if you have no idea what I'm talking about) then just create an array of the complex type (which is a standard C99 data type) and copy the values to that array, which will set the imaginary components to zero automatically. If you just want the frequency amplitudes from the return data, then just take cabs(output[i]) for each item in the output array, which converts the complex number to a real number that represents the magnitude of that frequency in the input data. Compile with "-std=gnu99" and "-lm" options. Also, don't forget to use "-O3" and "--fast-math" when compiling any code you want to run fast. Released under the antivirial license. Basically, you can do anything you want with it as long as what you want doesn't involve the GNU GPL. See http://www.ecstaticlyrics.com/antiviral/ for more information. */ #include #include #include #include #include static int samples = 0; static int inverse = 0; static complex *table = NULL; static complex *data = NULL; static complex *temp = NULL; void fft_initialize (int s, int i) { if (s & (s - 1)) { fprintf(stderr, "FFT: Sample size is not a power of two.\n"); exit(1); }; if (s == samples && !inverse == !i) return; samples = s; inverse = i; table = realloc(table, samples * sizeof(complex)); if (table == NULL) { fprintf(stderr, "FFT: Failed to allocate memory for table.\n"); exit(1); }; data = realloc(data, samples * sizeof(complex)); if (data == NULL) { fprintf(stderr, "FFT: Failed to allocate memory for data table.\n"); exit(1); }; temp = realloc(temp, samples * sizeof(complex)); if (temp == NULL) { fprintf(stderr, "FFT: Failed to allocate memory for temp table.\n"); exit(1); }; for (int i = 0; i < samples; i++) { table[i] = cexp((inverse ? +1 : -1) * 2 * M_PI * I * i / samples); }; }; static void recursive (complex *data, int offset, int step, int count, complex *temp) { if (count > 2) { recursive(temp, offset, step << 1, count >> 1, data); recursive(temp, offset + step, step << 1, count >> 1, data); }; //printf("offset=%d count=%d step=%d:\n", offset, count, step); for (int i = 0; i < count >> 1; i++) { int a = offset + step * i; int b = offset + step * (2 * i); int c = offset + step * (2 * i + 1); int d = (step * i) & (samples - 1); data[a] = temp[b] + temp[c] * table[d]; //printf(" o[%d] = i[%d] + i[%d] * a(%d)\n", a, b, c, d); //printf(" (%+0.2f, %+0.2f) = (%+0.2f, %+0.2f) + (%+0.2f, %+0.2f) * (%+0.2f, %+0.2f)\n", creal(data[a]), cimag(data[a]), creal(temp[b]), cimag(temp[b]), creal(temp[c]), cimag(temp[c]), creal(table[d]), cimag(table[d])); }; for (int i = 0; i < count >> 1; i++) { int a = offset + step * (i + (count >> 1)); int b = offset + step * (2 * i); int c = offset + step * (2 * i + 1); int d = (step * (i + (count >> 1))) & (samples - 1); data[a] = temp[b] + temp[c] * table[d]; //printf(" o[%d] = i[%d] + i[%d] * a(%d)\n", a, b, c, d); //printf(" (%+0.2f, %+0.2f) = (%+0.2f, %+0.2f) + (%+0.2f, %+0.2f) * (%+0.2f, %+0.2f)\n", creal(data[a]), cimag(data[a]), creal(temp[b]), cimag(temp[b]), creal(temp[c]), cimag(temp[c]), creal(table[d]), cimag(table[d])); }; }; void fft_fast (complex *output, complex *input) { if (samples) { if ((int) trunc(log(samples) / log(2) + 0.5) & 1) { memmove(temp, input, samples * sizeof(complex)); } else { memmove(data, input, samples * sizeof(complex)); }; recursive(data, 0, 1, samples, temp); for (int i = 0; i < samples; i++) { if (inverse) { output[i] = data[i]; } else { output[i] = data[i] / samples; }; }; } else { fprintf(stderr, "FFT: fft_initialize() has not been called.\n"); exit(1); }; };