Mercurial > hg > vamp-plugin-sdk
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 |