Initial commit

AUTHORS

@@ -0,0 +1 @@
+David Barina <[email protected]>

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2017 David Barina
+
+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.

Makefile

@@ -0,0 +1,7 @@
+CFLAGS=-std=c99 -pedantic -Wall -Wextra -Wconversion -D_XOPEN_SOURCE -march=native -O3
+LDLIBS=-lm
+
+example: libufft.a
+
+libufft.a: ufft.o
+	$(AR) -rs $@ $^

README

@@ -0,0 +1,14 @@
+uFFT is a small portable C library for computing the discrete Fourier transform
+(DFT) in one dimension.
+
+The library implements forward fast Fourier transform (FFT) algorithm. The
+inverse transform is not implemented as it can be computed using the existing
+forward transform. For example, the inverse transform is the same as the forward
+one with the real and imaginary parts swapped for both input and output.
+
+The library is written in pure C99. No compiler extensions nor assembly language
+are employed. It uses floating point data type and can handle unaligned data.
+The FFT routines can be easily modified since their source code has less than a
+hundred lines. As you might expect, the uFFT performance does not outperform the
+performance of FFTW. The library uses one temporary array as large as the input
+array.

example.c

@@ -0,0 +1,25 @@
+#include <stddef.h>
+#include <complex.h>
+#include "ufft.h"
+#include <stdio.h>
+
+int main()
+{
+	size_t N = 1<<3;
+
+	float complex vector[N];
+
+	for(size_t n = 0; n < N; n++) {
+		vector[n] = n;
+	}
+
+	fft(vector, N);
+
+	printf("DFT:\n");
+
+	for(size_t n = 0; n < N; n++) {
+		printf("%f%+fi\n", creal(vector[n]), cimag(vector[n]));
+	}
+
+	return 0;
+}

ufft.c

@@ -0,0 +1,97 @@
+#include "ufft.h"
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+static int ctz(size_t N)
+{
+	int ctz1 = 0;
+
+	while( N ) {
+		ctz1++;
+		N >>= 1;
+	}
+
+	return ctz1-1;
+}
+
+static void nop_split(const float complex *x, float complex *X, size_t N)
+{
+	for(size_t n = 0; n < N/2; n++) {
+		X[0   + n] = x[2*n+0];
+		X[N/2 + n] = x[2*n+1];
+	}
+}
+
+static void fft_split(const float complex *x, float complex *X, size_t N, float complex phi)
+{
+	for(size_t n = 0; n < N/2; n++) {
+		X[0   + n] = x[2*n+0] + x[2*n+1] * cexp(-2*(float)M_PI*I*phi);
+		X[N/2 + n] = x[2*n+0] - x[2*n+1] * cexp(-2*(float)M_PI*I*phi);
+	}
+}
+
+static size_t revbits(size_t v, int J)
+{
+	size_t r = 0;
+
+	for(int j = 0; j < J; j++) {
+		r |= ( (v>>j)&1 ) << (J-1-j);
+	}
+
+	return r;
+}
+
+static int nop_reverse(int b, float complex *buffers[2], size_t N)
+{
+	int J = ctz(N);
+
+	for(int j = 0; j < J-1; j++, b++) {
+		size_t delta = N>>j;
+
+		for(size_t n = 0; n < N; n += delta) {
+			nop_split(buffers[b&1]+n, buffers[~b&1]+n, delta);
+		}
+	}
+
+	return b;
+}
+
+static int fft_reverse(int b, float complex *buffers[2], size_t N)
+{
+	int J = ctz(N);
+
+	for(int j = 0; j < J; j++, b++) {
+		size_t delta = N>>j;
+
+		for(size_t n = 0; n < N; n += delta) {
+			float complex phi = (float)revbits( n/delta, j) / (float)(2<<j);
+			fft_split(buffers[b&1]+n, buffers[~b&1]+n, delta, phi);
+		}
+	}
+
+	return b;
+}
+
+int fft(float complex *vector, size_t N)
+{
+	if( !N ) return 0;
+
+	if( N & (N-1) ) return 1;
+
+	float complex *buffers[2] = { vector, malloc(N*sizeof(float complex)) };
+
+	if( !buffers[1] ) return -1;
+
+	int b = 0;
+
+	b = nop_reverse(b, buffers, N);
+	b = fft_reverse(b, buffers, N);
+	b = nop_reverse(b, buffers, N);
+
+	memmove(vector, buffers[b&1], N*sizeof(float complex));
+
+	free( buffers[1] );
+
+	return 0;
+}

ufft.h

@@ -0,0 +1,20 @@
+#ifndef UFFT_H
+#define UFFT_H
+
+#include <complex.h>
+#include <stddef.h>
+
+/**
+ * @brief FFT algorithm
+ *
+ * This function computes forward radix-2 fast Fourier transform (FFT) using Cooley-Tukey algorithm.
+ * The output is written in-place over the input.
+ *
+ * @param vector An array of @p N complex values in single-precision floating-point format.
+ * @param N The size of the transform must be a power of two.
+ *
+ * @return Zero for success.
+ */
+int fft(float complex *vector, size_t N);
+
+#endif