comparison 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
comparison
equal deleted inserted replaced
336:50df48a51c97 337:d5c5a52e6c9f
1
2 /* Public domain FFT implementation from Don Cross. */
3
4 static void
5 fft(unsigned int n, bool inverse,
6 const double *ri, const double *ii,
7 double *ro, double *io)
8 {
9 if (!ri || !ro || !io) return;
10
11 unsigned int bits;
12 unsigned int i, j, k, m;
13 unsigned int blockSize, blockEnd;
14
15 double tr, ti;
16
17 if (n < 2) return;
18 if (n & (n-1)) return;
19
20 double angle = 2.0 * M_PI;
21 if (inverse) angle = -angle;
22
23 for (i = 0; ; ++i) {
24 if (n & (1 << i)) {
25 bits = i;
26 break;
27 }
28 }
29
30 int table[n];
31
32 for (i = 0; i < n; ++i) {
33 m = i;
34 for (j = k = 0; j < bits; ++j) {
35 k = (k << 1) | (m & 1);
36 m >>= 1;
37 }
38 table[i] = k;
39 }
40
41 if (ii) {
42 for (i = 0; i < n; ++i) {
43 ro[table[i]] = ri[i];
44 io[table[i]] = ii[i];
45 }
46 } else {
47 for (i = 0; i < n; ++i) {
48 ro[table[i]] = ri[i];
49 io[table[i]] = 0.0;
50 }
51 }
52
53 blockEnd = 1;
54
55 for (blockSize = 2; blockSize <= n; blockSize <<= 1) {
56
57 double delta = angle / (double)blockSize;
58 double sm2 = -sin(-2 * delta);
59 double sm1 = -sin(-delta);
60 double cm2 = cos(-2 * delta);
61 double cm1 = cos(-delta);
62 double w = 2 * cm1;
63 double ar[3], ai[3];
64
65 for (i = 0; i < n; i += blockSize) {
66
67 ar[2] = cm2;
68 ar[1] = cm1;
69
70 ai[2] = sm2;
71 ai[1] = sm1;
72
73 for (j = i, m = 0; m < blockEnd; j++, m++) {
74
75 ar[0] = w * ar[1] - ar[2];
76 ar[2] = ar[1];
77 ar[1] = ar[0];
78
79 ai[0] = w * ai[1] - ai[2];
80 ai[2] = ai[1];
81 ai[1] = ai[0];
82
83 k = j + blockEnd;
84 tr = ar[0] * ro[k] - ai[0] * io[k];
85 ti = ar[0] * io[k] + ai[0] * ro[k];
86
87 ro[k] = ro[j] - tr;
88 io[k] = io[j] - ti;
89
90 ro[j] += tr;
91 io[j] += ti;
92 }
93 }
94
95 blockEnd = blockSize;
96 }
97
98 if (inverse) {
99
100 double denom = (double)n;
101
102 for (i = 0; i < n; i++) {
103 ro[i] /= denom;
104 io[i] /= denom;
105 }
106 }
107 }
108