aboutsummaryrefslogtreecommitdiffstats
path: root/common
diff options
context:
space:
mode:
Diffstat (limited to 'common')
-rw-r--r--common/alcomplex.c62
-rw-r--r--common/alcomplex.h62
2 files changed, 124 insertions, 0 deletions
diff --git a/common/alcomplex.c b/common/alcomplex.c
new file mode 100644
index 00000000..c1a312e3
--- /dev/null
+++ b/common/alcomplex.c
@@ -0,0 +1,62 @@
+
+#include "config.h"
+
+#include "alcomplex.h"
+
+#include <math.h>
+
+
+extern inline ALcomplex complex_add(ALcomplex a, ALcomplex b);
+extern inline ALcomplex complex_sub(ALcomplex a, ALcomplex b);
+extern inline ALcomplex complex_mult(ALcomplex a, ALcomplex b);
+
+void complex_fft(ALcomplex *FFTBuffer, ALsizei FFTSize, ALdouble Sign)
+{
+ ALsizei i, j, k, mask, step, step2;
+ ALcomplex temp, u, w;
+ ALdouble arg;
+
+ /* Bit-reversal permutation applied to a sequence of FFTSize items */
+ for(i = 1;i < FFTSize-1;i++)
+ {
+ for(mask = 0x1, j = 0;mask < FFTSize;mask <<= 1)
+ {
+ if((i&mask) != 0)
+ j++;
+ j <<= 1;
+ }
+ j >>= 1;
+
+ if(i < j)
+ {
+ temp = FFTBuffer[i];
+ FFTBuffer[i] = FFTBuffer[j];
+ FFTBuffer[j] = temp;
+ }
+ }
+
+ /* Iterative form of Danielson–Lanczos lemma */
+ for(i = 1, step = 2;i < FFTSize;i<<=1, step<<=1)
+ {
+ step2 = step >> 1;
+ arg = M_PI / step2;
+
+ w.Real = cos(arg);
+ w.Imag = sin(arg) * Sign;
+
+ u.Real = 1.0;
+ u.Imag = 0.0;
+
+ for(j = 0;j < step2;j++)
+ {
+ for(k = j;k < FFTSize;k+=step)
+ {
+ temp = complex_mult(FFTBuffer[k+step2], u);
+ FFTBuffer[k+step2] = complex_sub(FFTBuffer[k], temp);
+ FFTBuffer[k] = complex_add(FFTBuffer[k], temp);
+ }
+
+ u = complex_mult(u, w);
+ }
+ }
+}
diff --git a/common/alcomplex.h b/common/alcomplex.h
new file mode 100644
index 00000000..cf4683fa
--- /dev/null
+++ b/common/alcomplex.h
@@ -0,0 +1,62 @@
+#ifndef ALCOMPLEX_H
+#define ALCOMPLEX_H
+
+#include "AL/al.h"
+
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct ALcomplex {
+ ALdouble Real;
+ ALdouble Imag;
+} ALcomplex;
+
+/** Addition of two complex numbers. */
+inline ALcomplex complex_add(ALcomplex a, ALcomplex b)
+{
+ ALcomplex result;
+
+ result.Real = a.Real + b.Real;
+ result.Imag = a.Imag + b.Imag;
+
+ return result;
+}
+
+/** Subtraction of two complex numbers. */
+inline ALcomplex complex_sub(ALcomplex a, ALcomplex b)
+{
+ ALcomplex result;
+
+ result.Real = a.Real - b.Real;
+ result.Imag = a.Imag - b.Imag;
+
+ return result;
+}
+
+/** Multiplication of two complex numbers. */
+inline ALcomplex complex_mult(ALcomplex a, ALcomplex b)
+{
+ ALcomplex result;
+
+ result.Real = a.Real*b.Real - a.Imag*b.Imag;
+ result.Imag = a.Imag*b.Real + a.Real*b.Imag;
+
+ return result;
+}
+
+/**
+ * Iterative implementation of 2-radix FFT (In-place algorithm). Sign = -1 is
+ * FFT and 1 is iFFT (inverse). Fills FFTBuffer[0...FFTSize-1] with the
+ * Discrete Fourier Transform (DFT) of the time domain data stored in
+ * FFTBuffer[0...FFTSize-1]. FFTBuffer is an array of complex numbers, FFTSize
+ * MUST BE power of two.
+ */
+void complex_fft(ALcomplex *FFTBuffer, ALsizei FFTSize, ALdouble Sign);
+
+#ifdef __cplusplus
+} // extern "C"
+#endif
+
+#endif /* ALCOMPLEX_H */