Mercurial > hg > vamp-plugin-sdk
diff src/vamp-sdk/FFTimpl.cpp @ 337:d5c5a52e6c9f
Make the simple base-fft implementation accessible for use by plugins as well. Bump version to 2.4
author | Chris Cannam |
---|---|
date | Thu, 12 Jul 2012 11:37:31 +0100 |
parents | |
children | ab8e761f3ee9 |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/vamp-sdk/FFTimpl.cpp Thu Jul 12 11:37:31 2012 +0100 @@ -0,0 +1,108 @@ + +/* Public domain FFT implementation from Don Cross. */ + +static void +fft(unsigned int n, bool inverse, + const double *ri, const double *ii, + double *ro, double *io) +{ + if (!ri || !ro || !io) return; + + unsigned int bits; + unsigned int i, j, k, m; + unsigned int blockSize, blockEnd; + + double tr, ti; + + if (n < 2) return; + if (n & (n-1)) return; + + double angle = 2.0 * M_PI; + if (inverse) angle = -angle; + + for (i = 0; ; ++i) { + if (n & (1 << i)) { + bits = i; + break; + } + } + + int table[n]; + + for (i = 0; i < n; ++i) { + m = i; + for (j = k = 0; j < bits; ++j) { + k = (k << 1) | (m & 1); + m >>= 1; + } + table[i] = k; + } + + if (ii) { + for (i = 0; i < n; ++i) { + ro[table[i]] = ri[i]; + io[table[i]] = ii[i]; + } + } else { + for (i = 0; i < n; ++i) { + ro[table[i]] = ri[i]; + io[table[i]] = 0.0; + } + } + + blockEnd = 1; + + for (blockSize = 2; blockSize <= n; blockSize <<= 1) { + + double delta = angle / (double)blockSize; + double sm2 = -sin(-2 * delta); + double sm1 = -sin(-delta); + double cm2 = cos(-2 * delta); + double cm1 = cos(-delta); + double w = 2 * cm1; + double ar[3], ai[3]; + + for (i = 0; i < n; i += blockSize) { + + ar[2] = cm2; + ar[1] = cm1; + + ai[2] = sm2; + ai[1] = sm1; + + for (j = i, m = 0; m < blockEnd; j++, m++) { + + ar[0] = w * ar[1] - ar[2]; + ar[2] = ar[1]; + ar[1] = ar[0]; + + ai[0] = w * ai[1] - ai[2]; + ai[2] = ai[1]; + ai[1] = ai[0]; + + k = j + blockEnd; + tr = ar[0] * ro[k] - ai[0] * io[k]; + ti = ar[0] * io[k] + ai[0] * ro[k]; + + ro[k] = ro[j] - tr; + io[k] = io[j] - ti; + + ro[j] += tr; + io[j] += ti; + } + } + + blockEnd = blockSize; + } + + if (inverse) { + + double denom = (double)n; + + for (i = 0; i < n; i++) { + ro[i] /= denom; + io[i] /= denom; + } + } +} +