--- a/Makefile.gen	Sat May 04 19:08:01 2013 +0300
+++ b/Makefile.gen	Sat May 04 23:29:22 2013 +0300
@@ -270,7 +270,7 @@
 DMLIB_OBJS += \
 	dmfile.o dmbstr.o dmlib.o dmlerp.o dmstring.o \
 	dmargs.o dmvecmat.o dmperlin.o dmimage.o \
-	dmwav.o	dmengine.o dmq3d.o
+	dmwav.o	dmengine.o dmq3d.o dmfft.o
 
 
 TESTS_TARGETS = $(addprefix $(TESTS_BINPATH),$(addsuffix $(EXEEXT),$(TESTS_BINARIES)))

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dmfft.c	Sat May 04 23:29:22 2013 +0300
@@ -0,0 +1,151 @@
+/*
+ * Realfftf
+ * Originally (C) Copyright 1993 Philib Van Baren
+ * Cleaned up and refactored to be re-entrant by
+ *     Matti 'ccr' Hamalainen <ccr@tnsp.org>, 2013
+ *
+ * Note: Output is BIT-REVERSED! so you must use the BitReversed to
+ *       get legible output, (i.e. Real_i = buffer[ BitReversed[i] ]
+ *                                 Imag_i = buffer[ BitReversed[i]+1 ] )
+ *       Input is in normal order.
+ */
+#include "dmfft.h"
+#include <math.h>
+
+
+int dmInitializeFFT(DMFFTContext * ctx, const int fftlen)
+{
+    int i;
+
+    if (ctx == NULL)
+        return DMERR_NULLPTR;
+
+    if (fftlen <= 16)
+        return DMERR_INVALID_ARGS;
+
+    // Allocate necessary tables
+    ctx->npoints = fftlen / 2;
+    ctx->sinTable = (DMFFTType *) dmMalloc(fftlen * sizeof(DMFFTType));
+    if (ctx->sinTable == NULL)
+        return DMERR_MALLOC;
+
+    // Bitreversion buffer
+    ctx->breversed = (int *) dmMalloc(ctx->npoints * sizeof(int));
+    if (ctx->breversed == NULL)
+    {
+        dmFree(ctx->sinTable);
+        return DMERR_MALLOC;
+    }
+
+    // Calculate data
+    for (i = 0; i < ctx->npoints; i++)
+    {
+        int mask, tmp;
+        for (tmp = 0, mask = ctx->npoints / 2; mask > 0; mask >>= 1)
+        {
+            tmp = (tmp >> 1) + ((i & mask) ? ctx->npoints : 0);
+        }
+
+        ctx->breversed[i] = tmp;
+    }
+
+    for (i = 0; i < ctx->npoints; i++)
+    {
+        ctx->sinTable[ctx->breversed[i]] = -sin((2.0f * M_PI * i) / fftlen);
+        ctx->sinTable[ctx->breversed[i] + 1] =
+            -cos((2.0f * M_PI * i) / fftlen);
+    }
+
+    return DMERR_OK;
+}
+
+
+void dmEndFFT(DMFFTContext * ctx)
+{
+    if (ctx != NULL)
+    {
+        dmFree(ctx->breversed);
+        dmFree(ctx->sinTable);
+        memset(ctx, 0, sizeof(DMFFTContext));
+    }
+}
+
+
+int dmRealFFT(DMFFTContext * ctx, DMFFTType * buffer)
+{
+    if (ctx == NULL || buffer == NULL)
+        return DMERR_NULLPTR;
+
+    int rounds = ctx->npoints / 2;
+    DMFFTType *endptr1 = buffer + ctx->npoints * 2;
+
+    while (rounds > 0)
+    {
+        DMFFTType *sptr = ctx->sinTable;
+        DMFFTType *A = buffer,
+                  *B = buffer + rounds * 2;
+
+        while (A < endptr1)
+        {
+            DMFFTType qsin = *sptr;
+            DMFFTType qcos = *(sptr + 1);
+            DMFFTType *endptr2 = B;
+
+            while (A < endptr2)
+            {
+                DMFFTType v1 = (*B) * qcos + (*(B + 1)) * qsin;
+                DMFFTType v2 = (*B) * qsin - (*(B + 1)) * qcos;
+                *B = (*A + v1) * 0.5;
+                *(A++) = *(B++) - v1;
+                *B = (*A - v2) * 0.5;
+                *(A++) = *(B++) + v2;
+            }
+            A = B;
+            B += rounds * 2;
+            sptr += 2;
+        }
+
+        rounds >>= 1;
+    }
+
+    /*
+     *  Massage output to get the output for a real input sequence.
+     */
+    int *br1 = ctx->breversed + 1;
+    int *br2 = ctx->breversed + ctx->npoints - 1;
+
+    while (br1 <= br2)
+    {
+        const DMFFTType vsin = ctx->sinTable[*br1];
+        const DMFFTType vcos = ctx->sinTable[*br1 + 1];
+        DMFFTType *A = buffer + *br1,
+                  *B = buffer + *br2,
+                  HRplus, HRminus,
+                  HIplus, HIminus,
+                  tmp1, tmp2;
+
+        HRminus  = *A - *B;
+        HRplus   = HRminus + (*B * 2);
+
+        HIminus  = *(A + 1) - *(B + 1);
+        HIplus   = HIminus + (*(B + 1) * 2);
+
+        tmp1     = (vsin * HRminus) - (vcos * HIplus);
+        tmp2     = (vcos * HRminus) + (vsin * HIplus);
+
+        *A       = (HRplus + tmp1) * 0.5f;
+        *B       = *A - tmp1;
+
+        *(A + 1) = (HIminus + tmp2) * 0.5f;
+        *(B + 1) = (*(A + 1)) - HIminus;
+
+        br1++;
+        br2--;
+    }
+
+    // Handle DC bin separately
+    buffer[0] += buffer[1];
+    buffer[1] = 0;
+
+    return DMERR_OK;
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dmfft.h	Sat May 04 23:29:22 2013 +0300
@@ -0,0 +1,23 @@
+#ifndef DMFFT_H
+#define DMFFT_H 1
+
+#include "dmlib.h"
+
+
+typedef double DMFFTType;
+
+
+typedef struct
+{
+    int npoints;
+    DMFFTType *sinTable;
+    int *breversed;
+} DMFFTContext;
+
+
+int dmInitializeFFT(DMFFTContext *, int);
+void dmEndFFT(DMFFTContext *);
+int dmRealFFT(DMFFTContext *, DMFFTType *);
+
+
+#endif // DMFFT_H