diff options
Diffstat (limited to 'common')
-rw-r--r-- | common/alcomplex.c | 62 | ||||
-rw-r--r-- | common/alcomplex.h | 62 |
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 */ |