The Allegro Wiki is migrating to github at https://github.com/liballeg/allegro_wiki/wiki

FFT

I made two routines to perform forward and reverse FFT. That's Fast Fourier Transform, the heart of Digital Signal Processing.

<highlightSyntax language="C">/***************************************************************************/ /* fft.c */ /***************************************************************************/ /* Forward and Reverse Fast Fourier Transform routines. */ /* ==================================================== */ /* By Tobias Dammers (c) 2005. */ /* Feel free to use these routines as you see fit. */ /***************************************************************************/

/* Fourier Transform is the process of transforming signals from the time domain

* (levels over time) into the frequency domain (levels over frequency). The dis-
* crete-time version of this is called the Discrete Fourier Transform (DFT).
* The standard approach is a rather costly operation, but for blocks of data
* whose length is a power of 2, a dramatic optimization is possible, resulting
* in an algorithm called Fast Fourier Transform, FFT. This is a C implementation
* of the algorithm, and of its reverse process.
*/
1. include <math.h>
1. ifndef PI
2. define PI 3.1415926535897932384626433832795f
3. endif

// The forward FFT. rex holds the real part of the signal, imx the imaginary part. // Both arrays must be at least 2^m samples long. The result is returned in the // same two arrays, ie the transform is in-place. After transforming, rex holds // the real part of the frequency domain, and imx the imaginary part. // To calculate the frequency spectrum of an arbitrary block of wave data, put the // input in rex[], and set imx[] to all-zero. void fft(int m, float* rex, float* imx) { int n = 1 << m; int nm1 = n - 1; int nd2 = n / 2; int j = nd2; float tr, ti; int i; int l;

for (i = 1; i < n; ++i) { int k; if (i < j) { tr = rex[j]; ti = imx[j]; rex[j] = rex[i]; imx[j] = imx[i]; rex[i] = tr; imx[i] = ti; } k = nd2; while ((k <= j) && (k > 0)) { j -= k; k = k >> 1; } j += k; }

for (l = 1; l <= m; ++l) { int le = 1 << l; int le2 = le/2; float ur = 1.0f; float ui = 0.0f; float sr = cos(PI/(float)le2); float si = -sin(PI/(float)le2); for (j = 1; j <= le2; ++j) { for (int i = j-1; i <= nm1; i += le) { int ip = i + le2; tr = rex[ip] * ur - imx[ip] * ui; ti = rex[ip] * ui + imx[ip] * ur; rex[ip] = rex[i] - tr; imx[ip] = imx[i] - ti; rex[i] += tr; imx[i] += ti; } tr = ur; ur = tr * sr - ui * si; ui = tr * si + ui * sr; } } }

// The inverse FFT. This function takes advantage of the special symmetry of // Fourier pairs; all we need to do is scale the input, apply a forward FFT // then scale the output. // As a test, you can feed a signal through the forward fft, then feed the result // back through the inverse fft, giving you the original signal, give or take // some roundoff errors. void inv_fft(int m, float* rex, float* imx) { int n = 1 << m; int i; for (i = 0; i < n; ++i) imx[i] = -imx[i];

fft(m, rex, imx);

for (i = 0; i < n; ++i) { rex[i] /= (float)n; imx[i] /= -(float)n; } } </highlightSyntax>