/*
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);
};
};