Chris@25
|
1 /* Copyright (c) 2012, Jens Nockert <jens@ofmlabs.org>, Jussi Kalliokoski <jussi@ofmlabs.org>
|
Chris@25
|
2 * All rights reserved.
|
Chris@25
|
3 *
|
Chris@25
|
4 * Redistribution and use in source and binary forms, with or without
|
Chris@25
|
5 * modification, are permitted provided that the following conditions are met:
|
Chris@25
|
6 *
|
Chris@25
|
7 * 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
|
Chris@25
|
8 * 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
|
Chris@25
|
9 *
|
Chris@25
|
10 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
|
Chris@25
|
11 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
|
Chris@25
|
12 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
|
Chris@25
|
13 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
|
Chris@25
|
14 * ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
|
Chris@25
|
15 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
|
Chris@25
|
16 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
|
Chris@25
|
17 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
|
Chris@25
|
18 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
|
Chris@25
|
19 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
|
Chris@25
|
20 */
|
Chris@25
|
21
|
Chris@25
|
22 if (!FFT) {
|
Chris@25
|
23 var FFT = {}
|
Chris@25
|
24 }
|
Chris@25
|
25
|
Chris@25
|
26 void function (namespace) {
|
Chris@25
|
27 "use strict"
|
Chris@25
|
28
|
Chris@25
|
29 function forwardButterfly2(output, outputOffset, outputStride, input, inputOffset, inputStride, product, n, twiddle, fStride) {
|
Chris@25
|
30 var m = n / 2, q = n / product, old = product / 2
|
Chris@25
|
31
|
Chris@25
|
32 for (var i = 0; i < q; i++) {
|
Chris@25
|
33 var a0 = old * i
|
Chris@25
|
34 var a1 = a0 + m
|
Chris@25
|
35
|
Chris@25
|
36 var s0 = input[inputOffset + inputStride * a0]
|
Chris@25
|
37 var s1 = input[inputOffset + inputStride * a1]
|
Chris@25
|
38
|
Chris@25
|
39 var r0 = s0 + s1
|
Chris@25
|
40 var r1 = s0 - s1
|
Chris@25
|
41
|
Chris@25
|
42 var a0 = product * i
|
Chris@25
|
43 var a1 = a0 + product - 1
|
Chris@25
|
44
|
Chris@25
|
45 output[outputOffset + outputStride * a0] = r0
|
Chris@25
|
46 output[outputOffset + outputStride * a1] = r1
|
Chris@25
|
47 }
|
Chris@25
|
48
|
Chris@25
|
49 if (old == 1) { return }
|
Chris@25
|
50
|
Chris@25
|
51 for (var i = 0; i < old / 2; i++) {
|
Chris@25
|
52 var t1_r = twiddle[2 * ((-1) + (i))], t1_i = twiddle[2 * ((-1) + (i)) + 1]
|
Chris@25
|
53
|
Chris@25
|
54 for (var j = 0; j < q; j++) {
|
Chris@25
|
55 var a0 = j * old + 2 * i - 1
|
Chris@25
|
56 var a1 = a0 + m
|
Chris@25
|
57
|
Chris@25
|
58 var s0_r = input[2 * ((inputOffset) + (inputStride) * (a0))], s0_i = input[2 * ((inputOffset) + (inputStride) * (a0)) + 1]
|
Chris@25
|
59
|
Chris@25
|
60 var s1_r = input[2 * ((inputOffset) + (inputStride) * (a1))], s1_i = input[2 * ((inputOffset) + (inputStride) * (a1)) + 1]
|
Chris@25
|
61 var v1_r = s1_r * t1_r - s1_i * t1_i, v1_i = s1_r * t1_i + s1_i * t1_r
|
Chris@25
|
62
|
Chris@25
|
63 var r0_r = s0_r + v1_r, r0_i = s0_i + v1_i
|
Chris@25
|
64 var r1_r = s0_r - v1_r, r1_i = s0_i - v1_i; r1_i = -r1_i
|
Chris@25
|
65
|
Chris@25
|
66 var a0 = j * product + 2 * i - 1
|
Chris@25
|
67 var a1 = (j - 1) * product - 2 * i - 1
|
Chris@25
|
68
|
Chris@25
|
69 output[2 * ((outputOffset) + (outputStride) * (a0))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (a0)) + 1] = r0_i
|
Chris@25
|
70 output[2 * ((outputOffset) + (outputStride) * (a1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (a1)) + 1] = r1_i
|
Chris@25
|
71 }
|
Chris@25
|
72 }
|
Chris@25
|
73
|
Chris@25
|
74 if (old % 2 == 1) { return }
|
Chris@25
|
75
|
Chris@25
|
76 for (var i = 0; i < q; i++) {
|
Chris@25
|
77 var a0 = (i + 1) * old - 1
|
Chris@25
|
78 var a1 = a0 + m
|
Chris@25
|
79
|
Chris@25
|
80 var r0_r = input[2 * ((inputOffset) + (inputStride) * (a0))]
|
Chris@25
|
81 var r1_i = -input[2 * ((inputOffset) + (inputStride) * (a1))]
|
Chris@25
|
82
|
Chris@25
|
83 var a0 = i * product + old - 1
|
Chris@25
|
84
|
Chris@25
|
85 output[2 * ((outputOffset) + (outputStride) * (a0))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (a0)) + 1] = r0_i
|
Chris@25
|
86 }
|
Chris@25
|
87 }
|
Chris@25
|
88
|
Chris@25
|
89 function backwardButterfly2(output, outputOffset, outputStride, input, inputOffset, inputStride, product, n, twiddle, fStride) {
|
Chris@25
|
90 var m = n / 2, q = n / product, old = product / 2
|
Chris@25
|
91
|
Chris@25
|
92 for (var i = 0; i < q; i++) {
|
Chris@25
|
93 var a0 = (2 * i) * q
|
Chris@25
|
94 var a1 = (2 * i + 2) * q - 1
|
Chris@25
|
95
|
Chris@25
|
96 var s0 = input[inputOffset + inputStride * a0]
|
Chris@25
|
97 var s1 = input[inputOffset + inputStride * a1]
|
Chris@25
|
98
|
Chris@25
|
99 var r0 = s0 + s1
|
Chris@25
|
100 var r1 = s0 - s1
|
Chris@25
|
101
|
Chris@25
|
102 var a0 = q * i
|
Chris@25
|
103 var a1 = q * i + m
|
Chris@25
|
104
|
Chris@25
|
105 output[outputOffset + outputStride * a0] = r0
|
Chris@25
|
106 output[outputOffset + outputStride * a1] = r1
|
Chris@25
|
107 }
|
Chris@25
|
108
|
Chris@25
|
109 if (q == 1) { return }
|
Chris@25
|
110
|
Chris@25
|
111 for (var i = 0; i < q / 2; i++) {
|
Chris@25
|
112 var t1_r = twiddle[2 * ((-1) + (i))], t1_i = twiddle[2 * ((-1) + (i)) + 1]
|
Chris@25
|
113
|
Chris@25
|
114 for (var j = 0; j < old; j++) {
|
Chris@25
|
115 var a0 = 2 * j * q + 2 * i - 1
|
Chris@25
|
116 var a1 = 2 * (j + 1) * q - 2 * i - 1
|
Chris@25
|
117
|
Chris@25
|
118 var s0_r = input[2 * ((inputOffset) + (inputStride) * (a0))], s0_i = input[2 * ((inputOffset) + (inputStride) * (a0)) + 1]
|
Chris@25
|
119 var s1_r = input[2 * ((inputOffset) + (inputStride) * (a1))], s1_i = input[2 * ((inputOffset) + (inputStride) * (a1)) + 1]
|
Chris@25
|
120
|
Chris@25
|
121 var r0_r = s0_r + s1_r
|
Chris@25
|
122 var r0_i = s0_i - s1_i
|
Chris@25
|
123
|
Chris@25
|
124 var v1_r = s0_r - s1_r
|
Chris@25
|
125 var v1_i = s0_i + s1_i
|
Chris@25
|
126
|
Chris@25
|
127 var r1_r = v1_r * t1_r - v1_i * t1_i, r1_i = v1_r * t1_i + v1_i * t1_r
|
Chris@25
|
128
|
Chris@25
|
129 var a0 = j * q + 2 * i - 1
|
Chris@25
|
130 var a1 = a0 + m
|
Chris@25
|
131
|
Chris@25
|
132 output[2 * ((outputOffset) + (outputStride) * (a0))] = r0_r, output[2 * ((outputOffset) + (outputStride) * (a0)) + 1] = r0_i
|
Chris@25
|
133 output[2 * ((outputOffset) + (outputStride) * (a1))] = r1_r, output[2 * ((outputOffset) + (outputStride) * (a1)) + 1] = r1_i
|
Chris@25
|
134 }
|
Chris@25
|
135 }
|
Chris@25
|
136
|
Chris@25
|
137 if (q % 2 == 1) { return }
|
Chris@25
|
138
|
Chris@25
|
139 for (var i = 0; i < q; i++) {
|
Chris@25
|
140 var a0 = 2 * (i + 1) * q - 1
|
Chris@25
|
141
|
Chris@25
|
142 var r0_r = input[2 * ((inputOffset) + (inputStride) * (a0))], r0_i = input[2 * ((inputOffset) + (inputStride) * (a0)) + 1]
|
Chris@25
|
143
|
Chris@25
|
144 input[2 * ((inputOffset) + (inputStride) * (a0))] = 2 * r0_r
|
Chris@25
|
145 input[2 * ((inputOffset) + (inputStride) * (a1)) + 1] = -2 * r0_i
|
Chris@25
|
146 }
|
Chris@25
|
147 }
|
Chris@25
|
148
|
Chris@25
|
149 function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) {
|
Chris@25
|
150 var p = factors.shift()
|
Chris@25
|
151 var m = factors.shift()
|
Chris@25
|
152
|
Chris@25
|
153 if (m == 1) {
|
Chris@25
|
154 for (var i = 0; i < p * m; i++) {
|
Chris@25
|
155 var x0_r = f[2 * ((fOffset) + (fStride * inputStride) * (i))], x0_i = f[2 * ((fOffset) + (fStride * inputStride) * (i)) + 1]
|
Chris@25
|
156 output[2 * ((outputOffset) + (outputStride) * (i))] = x0_r, output[2 * ((outputOffset) + (outputStride) * (i)) + 1] = x0_i
|
Chris@25
|
157 }
|
Chris@25
|
158 } else {
|
Chris@25
|
159 for (var i = 0; i < p; i++) {
|
Chris@25
|
160 work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state)
|
Chris@25
|
161 }
|
Chris@25
|
162 }
|
Chris@25
|
163
|
Chris@25
|
164 switch (p) {
|
Chris@25
|
165 case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break
|
Chris@25
|
166 case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break
|
Chris@25
|
167 case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break
|
Chris@25
|
168 default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break
|
Chris@25
|
169 }
|
Chris@25
|
170 }
|
Chris@25
|
171
|
Chris@25
|
172 var real = function (n, inverse) {
|
Chris@25
|
173 var n = ~~n, inverse = !!inverse
|
Chris@25
|
174
|
Chris@25
|
175 if (n < 1) {
|
Chris@25
|
176 throw new RangeError("n is outside range, should be positive integer, was `" + n + "'")
|
Chris@25
|
177 }
|
Chris@25
|
178
|
Chris@25
|
179 var state = {
|
Chris@25
|
180 n: n,
|
Chris@25
|
181 inverse: inverse,
|
Chris@25
|
182
|
Chris@25
|
183 factors: [],
|
Chris@25
|
184 twiddle: [],
|
Chris@25
|
185 scratch: new Float64Array(n)
|
Chris@25
|
186 }
|
Chris@25
|
187
|
Chris@25
|
188 var t = new Float64Array(n)
|
Chris@25
|
189
|
Chris@25
|
190 var p = 4, v = Math.floor(Math.sqrt(n))
|
Chris@25
|
191
|
Chris@25
|
192 while (n > 1) {
|
Chris@25
|
193 while (n % p) {
|
Chris@25
|
194 switch (p) {
|
Chris@25
|
195 case 4: p = 2; break
|
Chris@25
|
196 case 2: p = 3; break
|
Chris@25
|
197 default: p += 2; break
|
Chris@25
|
198 }
|
Chris@25
|
199
|
Chris@25
|
200 if (p > v) {
|
Chris@25
|
201 p = n
|
Chris@25
|
202 }
|
Chris@25
|
203 }
|
Chris@25
|
204
|
Chris@25
|
205 n /= p
|
Chris@25
|
206
|
Chris@25
|
207 state.factors.push(p)
|
Chris@25
|
208 }
|
Chris@25
|
209
|
Chris@25
|
210 var theta = 2 * Math.PI / n, product = 1, twiddle = new Float64Array(n)
|
Chris@25
|
211
|
Chris@25
|
212 for (var i = 0, t = 0; i < state.factors.length; i++) {
|
Chris@25
|
213 var phase = theta * i, factor = state.factors[i]
|
Chris@25
|
214
|
Chris@25
|
215 var old = product, product = product * factor, q = n / product
|
Chris@25
|
216
|
Chris@25
|
217 state.twiddle.push(new Float64Array(twiddle, t))
|
Chris@25
|
218
|
Chris@25
|
219 if (inverse) {
|
Chris@25
|
220 var counter = q, multiplier = old
|
Chris@25
|
221 } else {
|
Chris@25
|
222 var counter = old, multiplier = q
|
Chris@25
|
223 }
|
Chris@25
|
224
|
Chris@25
|
225 for (var j = 1; j < factor; j++) {
|
Chris@25
|
226 var m = 0
|
Chris@25
|
227
|
Chris@25
|
228 for (var k = 1; k < counter / 2; k++, t++) {
|
Chris@25
|
229 m = (m + j * multiplier) % n
|
Chris@25
|
230
|
Chris@25
|
231 var phase = theta * m
|
Chris@25
|
232
|
Chris@25
|
233 t[2 * (i)] = Math.cos(phase)
|
Chris@25
|
234 t[2 * (i) + 1] = Math.sin(phase)
|
Chris@25
|
235 }
|
Chris@25
|
236 }
|
Chris@25
|
237 }
|
Chris@25
|
238
|
Chris@25
|
239 this.state = state
|
Chris@25
|
240 }
|
Chris@25
|
241
|
Chris@25
|
242 real.prototype.process = function(output, outputStride, input, inputStride) {
|
Chris@25
|
243 var outputStride = ~~outputStride, inputStride = ~~inputStride
|
Chris@25
|
244
|
Chris@25
|
245 if (outputStride < 1) {
|
Chris@25
|
246 throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'")
|
Chris@25
|
247 }
|
Chris@25
|
248
|
Chris@25
|
249 if (inputStride < 1) {
|
Chris@25
|
250 throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'")
|
Chris@25
|
251 }
|
Chris@25
|
252
|
Chris@25
|
253 var product = 1, state = 0, inverse = this.state.inverse
|
Chris@25
|
254
|
Chris@25
|
255 var n = this.state.n, factors = this.state.factors
|
Chris@25
|
256 var twiddle = this.state.twiddle, scratch = this.state.scratch
|
Chris@25
|
257
|
Chris@25
|
258 for (var i = 0; i < factors.length; i++) {
|
Chris@25
|
259 var factor = factors[i], old = product, product = product * factor
|
Chris@25
|
260
|
Chris@25
|
261 var q = n / product, fStride = Math.ceil(old / 2) - 1
|
Chris@25
|
262
|
Chris@25
|
263 if (state == 0) {
|
Chris@25
|
264 var inBuffer = input, inStride = inputStride
|
Chris@25
|
265
|
Chris@25
|
266 if (this.state.factors.length % 2 == 0) {
|
Chris@25
|
267 var outBuffer = scratch, outStride = 1, state = 1
|
Chris@25
|
268 } else {
|
Chris@25
|
269 var outBuffer = output, outStride = outputStride, state = 2
|
Chris@25
|
270 }
|
Chris@25
|
271 } else if (state == 1) {
|
Chris@25
|
272 var inBuffer = scratch, inStride = 1, outBuffer = output, outStride = outputStride, state = 2
|
Chris@25
|
273 } else if (state == 2) {
|
Chris@25
|
274 var inBuffer = output, inStride = outputStride, outBuffer = scratch, outStride = 1, state = 1
|
Chris@25
|
275 } else {
|
Chris@25
|
276 throw new RangeError("state somehow is not in the range (0 .. 2)")
|
Chris@25
|
277 }
|
Chris@25
|
278
|
Chris@25
|
279 if (inverse) {
|
Chris@25
|
280 switch (factor) {
|
Chris@25
|
281 case 2: backwardButterfly2(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
282 case 3: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
283 case 4: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
284 case 5: backwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
285 default: backwardButterfly(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
286 }
|
Chris@25
|
287 } else {
|
Chris@25
|
288 switch (factor) {
|
Chris@25
|
289 case 2: forwardButterfly2(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
290 case 3: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
291 case 4: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
292 case 5: forwardButterfly3(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
293 default: forwardButterfly(outBuffer, 0, outStride, inBuffer, 0, inStride, product, n, twiddle[i], fStride); break
|
Chris@25
|
294 }
|
Chris@25
|
295 }
|
Chris@25
|
296 }
|
Chris@25
|
297 }
|
Chris@25
|
298
|
Chris@25
|
299 namespace.real = real
|
Chris@25
|
300 }(FFT)
|