FFT & DFT

http://paulbourke.net/miscellaneous/dft/

http://math.ewha.ac.kr/~jylee/SciComp/sc-crandall/fft.c

http://covil.sdu.dk/publications/MDobroczynski2DFFT2006.pdf

http://www.nayuki.io/res/free-small-fft-in-multiple-languages/fft.c

/* 

 * Free FFT and convolution (C)

 * 

 * Copyright (c) 2014 Project Nayuki

 * http://www.nayuki.io/page/free-small-fft-in-multiple-languages

 * 

 * (MIT License)

 * Permission is hereby granted, free of charge, to any person obtaining a copy of

 * this software and associated documentation files (the "Software"), to deal in

 * the Software without restriction, including without limitation the rights to

 * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of

 * the Software, and to permit persons to whom the Software is furnished to do so,

 * subject to the following conditions:

 * - The above copyright notice and this permission notice shall be included in

 *   all copies or substantial portions of the Software.

 * - The Software is provided "as is", without warranty of any kind, express or

 *   implied, including but not limited to the warranties of merchantability,

 *   fitness for a particular purpose and noninfringement. In no event shall the

 *   authors or copyright holders be liable for any claim, damages or other

 *   liability, whether in an action of contract, tort or otherwise, arising from,

 *   out of or in connection with the Software or the use or other dealings in the

 *   Software.

 */

#include <math.h>

#include <stdlib.h>

#include <string.h>

#include "fft.h"

// Private function prototypes

static size_t reverse_bits(size_t x, unsigned int n);

static void *memdup(const void *src, size_t n);

#define SIZE_MAX ((size_t)-1)

int transform(double real[], double imag[], size_t n) {

if (n == 0)

return 1;

else if ((n & (n - 1)) == 0)  // Is power of 2

return transform_radix2(real, imag, n);

else  // More complicated algorithm for arbitrary sizes

return transform_bluestein(real, imag, n);

}

int inverse_transform(double real[], double imag[], size_t n) {

return transform(imag, real, n);

}

int transform_radix2(double real[], double imag[], size_t n) {

// Variables

int status = 0;

unsigned int levels;

double *cos_table, *sin_table;

size_t size;

size_t i;

// Compute levels = floor(log2(n))

{

size_t temp = n;

levels = 0;

while (temp > 1) {

levels++;

temp >>= 1;

}

if (1u << levels != n)

return 0;  // n is not a power of 2

}

// Trignometric tables

if (SIZE_MAX / sizeof(double) < n / 2)

return 0;

size = (n / 2) * sizeof(double);

cos_table = malloc(size);

sin_table = malloc(size);

if (cos_table == NULL || sin_table == NULL)

goto cleanup;

for (i = 0; i < n / 2; i++) {

cos_table[i] = cos(2 * M_PI * i / n);

sin_table[i] = sin(2 * M_PI * i / n);

}

// Bit-reversed addressing permutation

for (i = 0; i < n; i++) {

size_t j = reverse_bits(i, levels);

if (j > i) {

double temp = real[i];

real[i] = real[j];

real[j] = temp;

temp = imag[i];

imag[i] = imag[j];

imag[j] = temp;

}

}

// Cooley-Tukey decimation-in-time radix-2 FFT

for (size = 2; size <= n; size *= 2) {

size_t halfsize = size / 2;

size_t tablestep = n / size;

for (i = 0; i < n; i += size) {

size_t j;

size_t k;

for (j = i, k = 0; j < i + halfsize; j++, k += tablestep) {

double tpre =  real[j+halfsize] * cos_table[k] + imag[j+halfsize] * sin_table[k];

double tpim = -real[j+halfsize] * sin_table[k] + imag[j+halfsize] * cos_table[k];

real[j + halfsize] = real[j] - tpre;

imag[j + halfsize] = imag[j] - tpim;

real[j] += tpre;

imag[j] += tpim;

}

}

if (size == n)  // Prevent overflow in 'size *= 2'

break;

}

status = 1;

cleanup:

free(cos_table);

free(sin_table);

return status;

}

int transform_bluestein(double real[], double imag[], size_t n) {

// Variables

int status = 0;

double *cos_table, *sin_table;

double *areal, *aimag;

double *breal, *bimag;

double *creal, *cimag;

size_t m;

size_t size_n, size_m;

size_t i;

// Find a power-of-2 convolution length m such that m >= n * 2 + 1

{

size_t target;

if (n > (SIZE_MAX - 1) / 2)

return 0;

target = n * 2 + 1;

for (m = 1; m < target; m *= 2) {

if (SIZE_MAX / 2 < m)

return 0;

}

}

// Allocate memory

if (SIZE_MAX / sizeof(double) < n || SIZE_MAX / sizeof(double) < m)

return 0;

size_n = n * sizeof(double);

size_m = m * sizeof(double);

cos_table = malloc(size_n);

sin_table = malloc(size_n);

areal = calloc(m, sizeof(double));

aimag = calloc(m, sizeof(double));

breal = calloc(m, sizeof(double));

bimag = calloc(m, sizeof(double));

creal = malloc(size_m);

cimag = malloc(size_m);

if (cos_table == NULL || sin_table == NULL

|| areal == NULL || aimag == NULL

|| breal == NULL || bimag == NULL

|| creal == NULL || cimag == NULL)

goto cleanup;

// Trignometric tables

for (i = 0; i < n; i++) {

double temp = M_PI * (size_t)((unsigned long long)i * i % ((unsigned long long)n * 2)) / n;

// Less accurate version if long long is unavailable: double temp = M_PI * i * i / n;

cos_table[i] = cos(temp);

sin_table[i] = sin(temp);

}

// Temporary vectors and preprocessing

for (i = 0; i < n; i++) {

areal[i] =  real[i] * cos_table[i] + imag[i] * sin_table[i];

aimag[i] = -real[i] * sin_table[i] + imag[i] * cos_table[i];

}

breal[0] = cos_table[0];

bimag[0] = sin_table[0];

for (i = 1; i < n; i++) {

breal[i] = breal[m - i] = cos_table[i];

bimag[i] = bimag[m - i] = sin_table[i];

}

// Convolution

if (!convolve_complex(areal, aimag, breal, bimag, creal, cimag, m))

goto cleanup;

// Postprocessing

for (i = 0; i < n; i++) {

real[i] =  creal[i] * cos_table[i] + cimag[i] * sin_table[i];

imag[i] = -creal[i] * sin_table[i] + cimag[i] * cos_table[i];

}

status = 1;

// Deallocation

cleanup:

free(cimag);

free(creal);

free(bimag);

free(breal);

free(aimag);

free(areal);

free(sin_table);

free(cos_table);

return status;

}

int convolve_real(const double x[], const double y[], double out[], size_t n) {

double *ximag, *yimag, *zimag;

int status = 0;

ximag = calloc(n, sizeof(double));

yimag = calloc(n, sizeof(double));

zimag = calloc(n, sizeof(double));

if (ximag == NULL || yimag == NULL || zimag == NULL)

goto cleanup;

status = convolve_complex(x, ximag, y, yimag, out, zimag, n);

cleanup:

free(zimag);

free(yimag);

free(ximag);

return status;

}

int convolve_complex(const double xreal[], const double ximag[], const double yreal[], const double yimag[], double outreal[], double outimag[], size_t n) {

int status = 0;

size_t size;

size_t i;

double *xr, *xi, *yr, *yi;

if (SIZE_MAX / sizeof(double) < n)

return 0;

size = n * sizeof(double);

xr = memdup(xreal, size);

xi = memdup(ximag, size);

yr = memdup(yreal, size);

yi = memdup(yimag, size);

if (xr == NULL || xi == NULL || yr == NULL || yi == NULL)

goto cleanup;

if (!transform(xr, xi, n))

goto cleanup;

if (!transform(yr, yi, n))

goto cleanup;

for (i = 0; i < n; i++) {

double temp = xr[i] * yr[i] - xi[i] * yi[i];

xi[i] = xi[i] * yr[i] + xr[i] * yi[i];

xr[i] = temp;

}

if (!inverse_transform(xr, xi, n))

goto cleanup;

for (i = 0; i < n; i++) {  // Scaling (because this FFT implementation omits it)

outreal[i] = xr[i] / n;

outimag[i] = xi[i] / n;

}

status = 1;

cleanup:

free(yi);

free(yr);

free(xi);

free(xr);

return status;

}

static size_t reverse_bits(size_t x, unsigned int n) {

size_t result = 0;

unsigned int i;

for (i = 0; i < n; i++, x >>= 1)

result = (result << 1) | (x & 1);

return result;

}

static void *memdup(const void *src, size_t n) {

void *dest = malloc(n);

if (dest != NULL)

memcpy(dest, src, n);

return dest;

}