Simple FFT in C++

fft.cpp

/**
 * Implementation for FFT and IFFT functions
 *
 * By Josh Jeppson
 *
 * Note: This took me a grand total of 26 lines in Python
 * */

#include "fft.h"

#include <stdlib.h>
#include <math.h>


double *
fft(double * buffer, uint64_t size) {
	cplx * cplxBuffer = (cplx *) malloc(sizeof(cplx) * size);
	uint64_t i = 0;
	for (; i < size; i++) {
		cplxBuffer[i] = cplx(buffer[i], 0);
	}
	cplx * freqcplx = fft_cplx(cplxBuffer, size);
	free(cplxBuffer);
	double * fDomain = (double *) malloc(sizeof(double) * size);
	for (i = 0; i < size; i++) {
		fDomain[i] = std::abs(freqcplx[i]);
	}
	free(freqcplx);
	return fDomain;
}

double *
ifft(double * buffer, uint64_t size) {
	cplx * cplxBuffer = (cplx *) malloc(sizeof(cplx) * size);
	uint64_t i = 0;
	for (; i < size; i++) {
		cplxBuffer[i] = cplx(buffer[i], 0);
	}
	cplx * timecplx = ifft_cplx(cplxBuffer, size);
	free(cplxBuffer);
	double * tDomain = (double *) malloc(sizeof(double) * size);
	for (i = 0; i < size; i++) {
		tDomain[i] = std::abs(timecplx[i]);
	}
	free(timecplx);
	return tDomain;
}

cplx *
fft_cplx(cplx * buffer, uint64_t size) {
	// Create omega array for the first size primitive roots of unity
	cplx * omega = (cplx*) malloc(sizeof(cplx) * size);
	uint64_t i = 0;
	for (; i < size; i++) {
		omega[i] = cplx(
			cos(2.0 * M_PI * (float) i / (float) size)
			, sin(2.0 * M_PI * (float) i / (float) size)
		);
	}
	cplx * freqDomain = fft_helper(buffer, omega, size);
	free(omega);
	return freqDomain;
}

cplx *
ifft_cplx(cplx * buffer, uint64_t size) {
	// For the inverse fft, the angle in the cplx plane is inverted
	// from that of the regular fft
	cplx * omega = (cplx *) malloc(sizeof(cplx) * size);
	uint64_t i = 0;
	for (; i < size; i++) {
		omega[i] = cplx(
			cos(2.0 * M_PI * (float) i / (float) size)
			, -sin(2.0 * M_PI * (float) i / (float) size)
		);
	}
	cplx * timeDomain = fft_helper(buffer, omega, size);
	free(omega);
	// The IFFT must be reverse-scaled back down by n (the size)
	for (i = 0; i < size; i++) {
		timeDomain[i] /= size;
	}
	return timeDomain;
}

cplx *
fft_helper(cplx * buffer, cplx * omega, uint64_t size) {
	if (size == 1) {
		cplx * ret = new cplx;
		*ret = *buffer;
		return ret;
	}
	cplx * even = (cplx *) malloc(sizeof(cplx) * size / 2);
	cplx * odd  = (cplx *) malloc(sizeof(cplx) * size / 2);
	cplx * omegaHalf = (cplx *) malloc(sizeof(cplx) * size / 2);
	uint64_t i = 0;
	for (; i < size / 2; i++) {
		// Get the even and odd elements
		even[i] = buffer[i * 2];
		odd[i]  = buffer[i * 2 + 1];
		// Square the elements in the first half of the omega array
		omegaHalf[i] = omega[i] * omega[i];
	}

	cplx * solution_even = fft_helper(even, omegaHalf, size / 2);
	cplx * solution_odd  = fft_helper(odd,  omegaHalf, size / 2);
	// Free memory before we forget
	free(even);
	free(odd);
	free(omegaHalf);
	// Re-combine the even and odd solutions of the FFT
	cplx * solution = (cplx *) malloc(sizeof(cplx) * size);
	for (i = 0; i < size / 2; i++) {
		solution[i] = solution_even[i] + omega[i] * solution_odd[i];
		solution[i + (size / 2)] = solution_even[i] - omega[i] * solution_odd[i];
	}
	free(solution_even);
	free(solution_odd);
	return solution;
}

fft.h

/**
 * Tiny Fourier transform library originally for use in Nodesynth
 *
 * By Josh Jeppson
 * */

#ifndef FFT_H_INCLUDED
#define FFT_H_INCLUDED

#include <stdint.h>
#include <complex.h>

// typedef std::complex<double> cplx;
using cplx = std::complex<double>;
/**
 * Computes the DFT (Discrete Fourier Transform) using the FFT (Fast Fourier Transform) algorithm,
 * which is a divide and conquer algorithm of O(n) = n log(n). This wrapper takes a buffer in
 * double precision floating point. Converts to cplx type and calls fft_cplx. Requires
 * the buffer size to be a power of 2.
 *
 * @param buffer the buffer to transform into the frequency domain
 * @param size the number of elements in the buffer
 * @return the buffer transformed into frequency domain
 * */
double * fft(double * buffer, uint64_t size);

/**
 * Given the DFT (Discrete Fourier Transform), computes the inverse FFT, which is also an O(n) = n log(n)
 * algorithm. This wrapper takes a buffer in double precision floating point. Converts into cplx
 * type and calls ifft_cplx. Requires the buffer size to be a power of 2.
 *
 * @param buffer the buffer (in frequency domain) to transform into the time domain
 * @param size the number of elements in the buffer
 * @return the buffer transformed into time domain
 * */
double * ifft(double * buffer, uint64_t size);

// Helper functions using the cplx type

/**
 * Computes the DFT (Discrete Fourier Transform) using the FFT (Fast Fourier Transform) algorithm,
 * which is a divide and conquer algorithm of O(n) = n log(n). This takes elements of type cplx
 * and returns elements of type cplx. Requires the buffer size to be a power of 2.
 *
 * @param buffer the buffer to transform into the frequency domain
 * @param size the number of elements in the buffer
 * @return the buffer transformed into frequency domain
 * */
cplx * fft_cplx(cplx * buffer, uint64_t size);

/**
 * Given the DFT (Discrete Fourier Transform), computes the inverse FFT, which is also an O(n) = n log(n)
 * algorithm. This takes elements of type cplx and returns elements of type cplx. Requires the
 * buffer size to be a power of 2.
 *
 * @param buffer the buffer (in frequency domain) to transform into the time domain
 * @param size the number of elements in the buffer
 * @return the buffer transformed into time domain
 * */
cplx * ifft_cplx(cplx * buffer, uint64_t size);

/**
 * Helper function that creates the actual FFT/IFFT. This function should probably not be used by
 * the end user, as it should only be called by fft_cplx and ifft_cplx. Assumes that the
 * buffer size is a power of 2. This function takes an omega array, which should contain
 * the primitive roots of unity.
 *
 * @param buffer Buffer to convert either in or out of the frequency domain
 * @param omega Omega array--contains primitive roots of unity
 * @param size The size of both the buffer and omega array
 * */
cplx * fft_helper(cplx * buffer, cplx * omega, uint64_t size);

#endif // FFT_H_INCLUDED