comparison fft/fft.js/src/complex.erb.js @ 25:66f9fd5ac611

Bring in some more of the third-party code
author Chris Cannam
date Wed, 07 Oct 2015 13:46:38 +0100
parents
children
comparison
equal deleted inserted replaced
24:e705de983b67 25:66f9fd5ac611
1 <%= $:.unshift('.'); require "#{File.dirname(__FILE__)}/../src/complex.rb"; File.read "#{File.dirname(__FILE__)}/../LICENSE" %>
2
3 if (!FFT) {
4 var FFT = {}
5 }
6
7 void function (namespace) {
8 "use strict"
9
10 function butterfly2(output, outputOffset, outputStride, fStride, state, m) {
11 var t = state.twiddle
12
13 for (var i = 0; i < m; i++) {
14 <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %>
15 <%= load('s1', 'output', 'outputOffset', 'i + m', 'outputStride') %>
16
17 <%= load('t1', 't', 0, 'i', 'fStride') %>
18
19 <%= cmul('v1', 's1', 't1') %>
20
21 <%= cadd('r0', 's0', 'v1') %>
22 <%= csub('r1', 's0', 'v1') %>
23
24 <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %>
25 <%= store('r1', 'output', 'outputOffset', 'i + m', 'outputStride') %>
26 }
27 }
28
29 function butterfly3(output, outputOffset, outputStride, fStride, state, m) {
30 var t = state.twiddle
31 var m1 = m, m2 = 2 * m
32 var fStride1 = fStride, fStride2 = 2 * fStride
33
34 var e = <%= imag('t', 0, 'm', 'fStride') %>
35
36 for (var i = 0; i < m; i++) {
37 <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %>
38
39 <%= load('s1', 'output', 'outputOffset', 'i + m1', 'outputStride') %>
40 <%= load('t1', 't', 0, 'i', 'fStride1') %>
41 <%= cmul('v1', 's1', 't1') %>
42
43 <%= load('s2', 'output', 'outputOffset', 'i + m2', 'outputStride') %>
44 <%= load('t2', 't', 0, 'i', 'fStride2') %>
45 <%= cmul('v2', 's2', 't2') %>
46
47 <%= cadd('i0', 'v1', 'v2') %>
48
49 <%= cadd('r0', 's0', 'i0') %>
50 <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %>
51
52 var i1_r = s0_r - i0_r * 0.5
53 var i1_i = s0_i - i0_i * 0.5
54
55 var i2_r = (v1_r - v2_r) * e
56 var i2_i = (v1_i - v2_i) * e
57
58 var r1_r = i1_r - i2_i
59 var r1_i = i1_i + i2_r
60 <%= store('r1', 'output', 'outputOffset', 'i + m1', 'outputStride') %>
61
62 var r2_r = i1_r + i2_i
63 var r2_i = i1_i - i2_r
64 <%= store('r2', 'output', 'outputOffset', 'i + m2', 'outputStride') %>
65 }
66 }
67
68 function butterfly4(output, outputOffset, outputStride, fStride, state, m) {
69 var t = state.twiddle
70 var m1 = m, m2 = 2 * m, m3 = 3 * m
71 var fStride1 = fStride, fStride2 = 2 * fStride, fStride3 = 3 * fStride
72
73 for (var i = 0; i < m; i++) {
74 <%= load('s0', 'output', 'outputOffset', 'i', 'outputStride') %>
75
76 <%= load('s1', 'output', 'outputOffset', 'i + m1', 'outputStride') %>
77 <%= load('t1', 't', 0, 'i', 'fStride1') %>
78 <%= cmul('v1', 's1', 't1') %>
79
80 <%= load('s2', 'output', 'outputOffset', 'i + m2', 'outputStride') %>
81 <%= load('t2', 't', 0, 'i', 'fStride2') %>
82 <%= cmul('v2', 's2', 't2') %>
83
84 <%= load('s3', 'output', 'outputOffset', 'i + m3', 'outputStride') %>
85 <%= load('t3', 't', 0, 'i', 'fStride3') %>
86 <%= cmul('v3', 's3', 't3') %>
87
88 <%= cadd('i0', 's0', 'v2') %>
89 <%= csub('i1', 's0', 'v2') %>
90 <%= cadd('i2', 'v1', 'v3') %>
91 <%= csub('i3', 'v1', 'v3') %>
92
93 <%= cadd('r0', 'i0', 'i2') %>
94
95 if (state.inverse) {
96 var r1_r = i1_r - i3_i
97 var r1_i = i1_i + i3_r
98 } else {
99 var r1_r = i1_r + i3_i
100 var r1_i = i1_i - i3_r
101 }
102
103 <%= csub('r2', 'i0', 'i2') %>
104
105 if (state.inverse) {
106 var r3_r = i1_r + i3_i
107 var r3_i = i1_i - i3_r
108 } else {
109 var r3_r = i1_r - i3_i
110 var r3_i = i1_i + i3_r
111 }
112
113 <%= store('r0', 'output', 'outputOffset', 'i', 'outputStride') %>
114 <%= store('r1', 'output', 'outputOffset', 'i + m1', 'outputStride') %>
115 <%= store('r2', 'output', 'outputOffset', 'i + m2', 'outputStride') %>
116 <%= store('r3', 'output', 'outputOffset', 'i + m3', 'outputStride') %>
117 }
118 }
119
120 function butterfly(output, outputOffset, outputStride, fStride, state, m, p) {
121 var t = state.twiddle, n = state.n, scratch = new Float64Array(2 * p)
122
123 for (var u = 0; u < m; u++) {
124 for (var q1 = 0, k = u; q1 < p; q1++, k += m) {
125 <%= load('x0', 'output', 'outputOffset', 'k', 'outputStride') %>
126 <%= store('x0', 'scratch', 'q1') %>
127 }
128
129 for (var q1 = 0, k = u; q1 < p; q1++, k += m) {
130 var tOffset = 0
131
132 <%= load('x0', 'scratch', 0) %>
133 <%= store('x0', 'output', 'outputOffset', 'k', 'outputStride') %>
134
135 for (var q = 1; q < p; q++) {
136 tOffset = (tOffset + fStride * k) % n
137
138 <%= load('s0', 'output', 'outputOffset', 'k', 'outputStride') %>
139
140 <%= load('s1', 'scratch', 'q') %>
141 <%= load('t1', 't', 'tOffset') %>
142 <%= cmul('v1', 's1', 't1') %>
143
144 <%= cadd('r0', 's0', 'v1') %>
145 <%= store('r0', 'output', 'outputOffset', 'k', 'outputStride') %>
146 }
147 }
148 }
149 }
150
151 function work(output, outputOffset, outputStride, f, fOffset, fStride, inputStride, factors, state) {
152 var p = factors.shift()
153 var m = factors.shift()
154
155 if (m == 1) {
156 for (var i = 0; i < p * m; i++) {
157 <%= load('x0', 'f', 'fOffset', 'i', 'fStride * inputStride') %>
158 <%= store('x0', 'output', 'outputOffset', 'i', 'outputStride') %>
159 }
160 } else {
161 for (var i = 0; i < p; i++) {
162 work(output, outputOffset + outputStride * i * m, outputStride, f, fOffset + i * fStride * inputStride, fStride * p, inputStride, factors.slice(), state)
163 }
164 }
165
166 switch (p) {
167 case 2: butterfly2(output, outputOffset, outputStride, fStride, state, m); break
168 case 3: butterfly3(output, outputOffset, outputStride, fStride, state, m); break
169 case 4: butterfly4(output, outputOffset, outputStride, fStride, state, m); break
170 default: butterfly(output, outputOffset, outputStride, fStride, state, m, p); break
171 }
172 }
173
174 var complex = function (n, inverse) {
175 if (arguments.length < 2) {
176 throw new RangeError("You didn't pass enough arguments, passed `" + arguments.length + "'")
177 }
178
179 var n = ~~n, inverse = !!inverse
180
181 if (n < 1) {
182 throw new RangeError("n is outside range, should be positive integer, was `" + n + "'")
183 }
184
185 var state = {
186 n: n,
187 inverse: inverse,
188
189 factors: [],
190 twiddle: new Float64Array(2 * n),
191 scratch: new Float64Array(2 * n)
192 }
193
194 var t = state.twiddle, theta = 2 * Math.PI / n
195
196 for (var i = 0; i < n; i++) {
197 if (inverse) {
198 var phase = theta * i
199 } else {
200 var phase = -theta * i
201 }
202
203 <%= real('t', 'i') %> = Math.cos(phase)
204 <%= imag('t', 'i') %> = Math.sin(phase)
205 }
206
207 var p = 4, v = Math.floor(Math.sqrt(n))
208
209 while (n > 1) {
210 while (n % p) {
211 switch (p) {
212 case 4: p = 2; break
213 case 2: p = 3; break
214 default: p += 2; break
215 }
216
217 if (p > v) {
218 p = n
219 }
220 }
221
222 n /= p
223
224 state.factors.push(p)
225 state.factors.push(n)
226 }
227
228 this.state = state
229 }
230
231 complex.prototype.simple = function (output, input, t) {
232 this.process(output, 0, 1, input, 0, 1, t)
233 }
234
235 complex.prototype.process = function(output, outputOffset, outputStride, input, inputOffset, inputStride, t) {
236 var outputStride = ~~outputStride, inputStride = ~~inputStride
237
238 var type = t == 'real' ? t : 'complex'
239
240 if (outputStride < 1) {
241 throw new RangeError("outputStride is outside range, should be positive integer, was `" + outputStride + "'")
242 }
243
244 if (inputStride < 1) {
245 throw new RangeError("inputStride is outside range, should be positive integer, was `" + inputStride + "'")
246 }
247
248 if (type == 'real') {
249 for (var i = 0; i < this.state.n; i++) {
250 var x0_r = input[inputOffset + inputStride * i]
251 var x0_i = 0.0
252
253 <%= store('x0', 'this.state.scratch', 'i') %>
254 }
255
256 work(output, outputOffset, outputStride, this.state.scratch, 0, 1, 1, this.state.factors.slice(), this.state)
257 } else {
258 if (input == output) {
259 work(this.state.scratch, 0, 1, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state)
260
261 for (var i = 0; i < this.state.n; i++) {
262 <%= load('x0', 'this.state.scratch', 'i') %>
263
264 <%= store('x0', 'output', 'outputOffset', 'i', 'outputStride') %>
265 }
266 } else {
267 work(output, outputOffset, outputStride, input, inputOffset, 1, inputStride, this.state.factors.slice(), this.state)
268 }
269 }
270 }
271
272 namespace.complex = complex
273 }(FFT)