forked from slarmoo/slarmoosboxtesting
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathFFT.ts
362 lines (343 loc) · 14.7 KB
/
FFT.ts
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
// Copyright (c) 2012-2022 John Nesky and contributing authors, distributed under the MIT license, see accompanying the LICENSE.md file.
// interface shared by number[], Float32Array, and other typed arrays in JavaScript.
interface NumberArray {
length: number;
[index: number]: number;
}
// A basic FFT operation scales the overall magnitude of elements by the
// square root of the length of the array, √N. Performing a forward FFT and
// then an inverse FFT results in the original array, but multiplied by N.
// This helper function can be used to compensate for that.
export function scaleElementsByFactor(array: NumberArray, factor: number): void {
for (let i: number = 0; i < array.length; i++) {
array[i] *= factor;
}
}
function isPowerOf2(n: number): boolean {
return !!n && !(n & (n - 1));
}
function countBits(n: number): number {
if (!isPowerOf2(n)) throw new Error("FFT array length must be a power of 2.");
return Math.round(Math.log(n) / Math.log(2));
}
// Rearranges the elements of the array, swapping the element at an index
// with an element at an index that is the bitwise reverse of the first
// index in base 2. Useful for computing the FFT.
function reverseIndexBits(array: NumberArray, fullArrayLength: number): void {
const bitCount: number = countBits(fullArrayLength);
if (bitCount > 16) throw new Error("FFT array length must not be greater than 2^16.");
const finalShift: number = 16 - bitCount;
for (let i: number = 0; i < fullArrayLength; i++) {
// Dear Javascript: Please support bit order reversal intrinsics. Thanks! :D
let j: number;
j = ((i & 0xaaaa) >> 1) | ((i & 0x5555) << 1);
j = ((j & 0xcccc) >> 2) | ((j & 0x3333) << 2);
j = ((j & 0xf0f0) >> 4) | ((j & 0x0f0f) << 4);
j = ((j >> 8) | ((j & 0xff) << 8)) >> finalShift;
if (j > i) {
let temp: number = array[i];
array[i] = array[j];
array[j] = temp;
}
}
}
// Provided for educational purposes. Easier to read than
// fastFourierTransform(), but computes the same result.
// Takes two parallel arrays representing the real and imaginary elements,
// respectively, and returns an array containing two new arrays, which
// contain the complex result of the transform.
export function discreteFourierTransform(realArray: NumberArray, imagArray: NumberArray): number[][] {
const fullArrayLength: number = realArray.length;
if (fullArrayLength != imagArray.length) throw new Error("FFT arrays must be the same length.");
const realOut: number[] = [];
const imagOut: number[] = [];
for (let i: number = 0; i < fullArrayLength; i++) {
realOut[i] = 0.0;
imagOut[i] = 0.0;
for (let j: number = 0; j < fullArrayLength; j++) {
const radians: number = -6.2831853 * j * i / fullArrayLength;
const c: number = Math.cos(radians);
const s: number = Math.sin(radians);
realOut[i] += realArray[j] * c - imagArray[j] * s;
imagOut[i] += realArray[j] * s + imagArray[j] * c;
}
}
return [realOut, imagOut];
}
// Performs a Fourier transform in O(N log(N)) operations. Overwrites the
// input real and imaginary arrays. Can be used for both forward and inverse
// transforms: swap the order of the arguments for the inverse.
export function fastFourierTransform(realArray: NumberArray, imagArray: NumberArray): void {
const fullArrayLength: number = realArray.length;
if (!isPowerOf2(fullArrayLength)) throw new Error("FFT array length must be a power of 2.");
if (fullArrayLength < 4) throw new Error("FFT array length must be at least 4.");
if (fullArrayLength != imagArray.length) throw new Error("FFT arrays must be the same length.");
reverseIndexBits(realArray, fullArrayLength);
reverseIndexBits(imagArray, fullArrayLength);
// First two passes, with strides of 2 and 4, can be combined and optimized.
for (let startIndex: number = 0; startIndex < fullArrayLength; startIndex += 4) {
const startIndex1: number = startIndex + 1;
const startIndex2: number = startIndex + 2;
const startIndex3: number = startIndex + 3;
const real0: number = realArray[startIndex ];
const real1: number = realArray[startIndex1];
const real2: number = realArray[startIndex2];
const real3: number = realArray[startIndex3];
const imag0: number = imagArray[startIndex ];
const imag1: number = imagArray[startIndex1];
const imag2: number = imagArray[startIndex2];
const imag3: number = imagArray[startIndex3];
const realTemp0: number = real0 + real1;
const realTemp1: number = real0 - real1;
const realTemp2: number = real2 + real3;
const realTemp3: number = real2 - real3;
const imagTemp0: number = imag0 + imag1;
const imagTemp1: number = imag0 - imag1;
const imagTemp2: number = imag2 + imag3;
const imagTemp3: number = imag2 - imag3;
realArray[startIndex ] = realTemp0 + realTemp2;
realArray[startIndex1] = realTemp1 + imagTemp3;
realArray[startIndex2] = realTemp0 - realTemp2;
realArray[startIndex3] = realTemp1 - imagTemp3;
imagArray[startIndex ] = imagTemp0 + imagTemp2;
imagArray[startIndex1] = imagTemp1 - realTemp3;
imagArray[startIndex2] = imagTemp0 - imagTemp2;
imagArray[startIndex3] = imagTemp1 + realTemp3;
}
for (let stride: number = 8; stride <= fullArrayLength; stride += stride) {
const halfLength: number = stride >>> 1;
const radiansIncrement: number = Math.PI * 2.0 / stride;
const cosIncrement: number = Math.cos(radiansIncrement);
const sinIncrement: number = Math.sin(radiansIncrement);
const oscillatorMultiplier: number = 2.0 * cosIncrement;
for (let startIndex: number = 0; startIndex < fullArrayLength; startIndex += stride) {
let c: number = 1.0;
let s: number = 0.0;
let cPrev: number = cosIncrement;
let sPrev: number = sinIncrement;
const secondHalf: number = startIndex + halfLength;
for (let i: number = startIndex; i < secondHalf; i++) {
const j: number = i + halfLength;
const real0: number = realArray[i];
const imag0: number = imagArray[i];
const real1: number = realArray[j] * c - imagArray[j] * s;
const imag1: number = realArray[j] * s + imagArray[j] * c;
realArray[i] = real0 + real1;
imagArray[i] = imag0 + imag1;
realArray[j] = real0 - real1;
imagArray[j] = imag0 - imag1;
const cTemp: number = oscillatorMultiplier * c - cPrev;
const sTemp: number = oscillatorMultiplier * s - sPrev;
cPrev = c;
sPrev = s;
c = cTemp;
s = sTemp;
}
}
}
}
// Computes the Fourier transform from an array of real-valued time-domain
// samples. The output is specially formatted for space efficieny: elements
// 0 through N/2 represent cosine wave amplitudes in ascending frequency,
// and elements N/2+1 through N-1 represent sine wave amplitudes in
// descending frequency. Overwrites the input array.
export function forwardRealFourierTransform(array: NumberArray): void {
const fullArrayLength: number = array.length;
const totalPasses: number = countBits(fullArrayLength);
if (fullArrayLength < 4) throw new Error("FFT array length must be at least 4.");
reverseIndexBits(array, fullArrayLength);
// First and second pass.
for (let index: number = 0; index < fullArrayLength; index += 4) {
const index1: number = index + 1;
const index2: number = index + 2;
const index3: number = index + 3;
const real0: number = array[index ];
const real1: number = array[index1];
const real2: number = array[index2];
const real3: number = array[index3];
// no imaginary elements yet since the input is fully real.
const tempA: number = real0 + real1;
const tempB: number = real2 + real3;
array[index ] = tempA + tempB;
array[index1] = real0 - real1;
array[index2] = tempA - tempB;
array[index3] = real2 - real3;
}
// Third pass.
const sqrt2over2: number = Math.sqrt(2.0) / 2.0;
for (let index: number = 0; index < fullArrayLength; index += 8) {
const index1: number = index + 1;
const index3: number = index + 3;
const index4: number = index + 4;
const index5: number = index + 5;
const index7: number = index + 7;
const real0: number = array[index ];
const real1: number = array[index1];
const imag3: number = array[index3];
const real4: number = array[index4];
const real5: number = array[index5];
const imag7: number = array[index7];
const tempA: number = (real5 - imag7) * sqrt2over2;
const tempB: number = (real5 + imag7) * sqrt2over2;
array[index ] = real0 + real4;
array[index1] = real1 + tempA;
array[index3] = real1 - tempA;
array[index4] = real0 - real4;
array[index5] = tempB - imag3;
array[index7] = tempB + imag3;
}
// Handle remaining passes.
for (let pass: number = 3; pass < totalPasses; pass++) {
const subStride: number = 1 << pass;
const midSubStride: number = subStride >> 1;
const stride: number = subStride << 1;
const radiansIncrement: number = Math.PI * 2.0 / stride;
const cosIncrement: number = Math.cos(radiansIncrement);
const sinIncrement: number = Math.sin(radiansIncrement);
const oscillatorMultiplier: number = 2.0 * cosIncrement;
for (let startIndex: number = 0; startIndex < fullArrayLength; startIndex += stride) {
const startIndexA: number = startIndex;
const startIndexB: number = startIndexA + subStride;
const stopIndex: number = startIndexB + subStride;
const realStartA: number = array[startIndexA];
const realStartB: number = array[startIndexB];
array[startIndexA] = realStartA + realStartB;
array[startIndexB] = realStartA - realStartB;
let c: number = cosIncrement;
let s: number = -sinIncrement;
let cPrev: number = 1.0;
let sPrev: number = 0.0;
for (let index: number = 1; index < midSubStride; index++) {
const indexA0: number = startIndexA + index;
const indexA1: number = startIndexB - index;
const indexB0: number = startIndexB + index;
const indexB1: number = stopIndex - index;
const real0: number = array[indexA0];
const imag0: number = array[indexA1];
const real1: number = array[indexB0];
const imag1: number = array[indexB1];
const tempA: number = real1 * c + imag1 * s;
const tempB: number = real1 * s - imag1 * c;
array[indexA0] = real0 + tempA;
array[indexA1] = real0 - tempA;
array[indexB0] =-imag0 - tempB;
array[indexB1] = imag0 - tempB;
const cTemp: number = oscillatorMultiplier * c - cPrev;
const sTemp: number = oscillatorMultiplier * s - sPrev;
cPrev = c;
sPrev = s;
c = cTemp;
s = sTemp;
}
}
}
}
// Computes the inverse Fourier transform from a specially formatted array of
// scalar values. Elements 0 through N/2 are expected to be the real values of
// the corresponding complex elements, representing cosine wave amplitudes in
// ascending frequency, and elements N/2+1 through N-1 correspond to the
// imaginary values, representing sine wave amplitudes in descending frequency.
// Generates real-valued time-domain samples. Overwrites the input array.
export function inverseRealFourierTransform(array: NumberArray, fullArrayLength: number): void {
const totalPasses: number = countBits(fullArrayLength);
if (fullArrayLength < 4) throw new Error("FFT array length must be at least 4.");
// Perform all but the last few passes in reverse.
for (let pass: number = totalPasses - 1; pass >= 2; pass--) {
const subStride: number = 1 << pass;
const midSubStride: number = subStride >> 1;
const stride: number = subStride << 1;
const radiansIncrement: number = Math.PI * 2.0 / stride;
const cosIncrement: number = Math.cos(radiansIncrement);
const sinIncrement: number = Math.sin(radiansIncrement);
const oscillatorMultiplier: number = 2.0 * cosIncrement;
for (let startIndex: number = 0; startIndex < fullArrayLength; startIndex += stride) {
const startIndexA: number = startIndex;
const midIndexA: number = startIndexA + midSubStride;
const startIndexB: number = startIndexA + subStride;
const midIndexB: number = startIndexB + midSubStride;
const stopIndex: number = startIndexB + subStride;
const realStartA: number = array[startIndexA];
const imagStartB: number = array[startIndexB];
array[startIndexA] = realStartA + imagStartB;
array[midIndexA] *= 2;
array[startIndexB] = realStartA - imagStartB;
array[midIndexB] *= 2;
let c: number = cosIncrement;
let s: number = -sinIncrement;
let cPrev: number = 1.0;
let sPrev: number = 0.0;
for (let index: number = 1; index < midSubStride; index++) {
const indexA0: number = startIndexA + index;
const indexA1: number = startIndexB - index;
const indexB0: number = startIndexB + index;
const indexB1: number = stopIndex - index;
const real0: number = array[indexA0];
const real1: number = array[indexA1];
const imag0: number = array[indexB0];
const imag1: number = array[indexB1];
const tempA: number = real0 - real1;
const tempB: number = imag0 + imag1;
array[indexA0] = real0 + real1;
array[indexA1] = imag1 - imag0;
array[indexB0] = tempA * c - tempB * s;
array[indexB1] = tempB * c + tempA * s;
const cTemp: number = oscillatorMultiplier * c - cPrev;
const sTemp: number = oscillatorMultiplier * s - sPrev;
cPrev = c;
sPrev = s;
c = cTemp;
s = sTemp;
}
}
}
/*
// Commented out this block (and compensated with an extra pass above)
// because it's slower in my testing so far.
// Pass with stride 8.
const sqrt2over2: number = Math.sqrt(2.0) / 2.0;
for (let index: number = 0; index < fullArrayLength; index += 8) {
const index1: number = index + 1;
const index2: number = index + 2;
const index3: number = index + 3;
const index4: number = index + 4;
const index5: number = index + 5;
const index6: number = index + 6;
const index7: number = index + 7;
const real0: number = array[index ];
const real1: number = array[index1];
const real2: number = array[index2];
const real3: number = array[index3];
const imag4: number = array[index4];
const imag5: number = array[index5];
const imag6: number = array[index6];
const imag7: number = array[index7];
const tempA: number = real1 - real3;
const tempB: number = imag5 + imag7;
array[index ] = real0 + imag4;
array[index1] = real1 + real3;
array[index2] = real2 * 2;
array[index3] = imag7 - imag5;
array[index4] = real0 - imag4;
array[index5] = (tempA + tempB) * sqrt2over2;
array[index6] = imag6 * 2;
array[index7] = (tempB - tempA) * sqrt2over2;
}
*/
// The final passes with strides 4 and 2, combined into one loop.
for (let index: number = 0; index < fullArrayLength; index += 4) {
const index1: number = index + 1;
const index2: number = index + 2;
const index3: number = index + 3;
const real0: number = array[index ];
const real1: number = array[index1] * 2;
const imag2: number = array[index2];
const imag3: number = array[index3] * 2;
const tempA: number = real0 + imag2;
const tempB: number = real0 - imag2;
array[index ] = tempA + real1;
array[index1] = tempA - real1;
array[index2] = tempB + imag3;
array[index3] = tempB - imag3;
}
reverseIndexBits(array, fullArrayLength);
}