OLD | NEW |
1 // FFT routines. Obtained from https://github.com/rtoy/js-hacks/blob/master/fft
/fft-r2-nn.js. | 1 // FFT routines. Obtained from |
| 2 // https://github.com/rtoy/js-hacks/blob/master/fft/fft-r2-nn.js. |
2 // | 3 // |
3 // Create an FFT object of order |order|. This will | 4 // Create an FFT object of order |order|. This will |
4 // compute forward and inverse FFTs of length 2^|order|. | 5 // compute forward and inverse FFTs of length 2^|order|. |
5 function FFT(order) { | 6 function FFT(order) { |
6 if (order <= 1) { | 7 if (order <= 1) { |
7 throw new this.FFTException(order); | 8 throw new this.FFTException(order); |
8 } | 9 } |
9 | 10 |
10 this.order = order; | 11 this.order = order; |
11 this.N = 1 << order; | 12 this.N = 1 << order; |
12 this.halfN = 1 << (order - 1); | 13 this.halfN = 1 << (order - 1); |
13 | 14 |
14 // Internal variables needed for computing each stage of the FFT. | 15 // Internal variables needed for computing each stage of the FFT. |
15 this.pairsInGroup = 0; | 16 this.pairsInGroup = 0; |
16 this.numberOfGroups = 0; | 17 this.numberOfGroups = 0; |
17 this.distance = 0; | 18 this.distance = 0; |
18 this.notSwitchInput = 0; | 19 this.notSwitchInput = 0; |
19 | 20 |
20 // Work arrays | 21 // Work arrays |
21 this.aReal = new Float32Array(this.N); | 22 this.aReal = new Float32Array(this.N); |
22 this.aImag = new Float32Array(this.N); | 23 this.aImag = new Float32Array(this.N); |
23 this.debug = 0; | 24 this.debug = 0; |
24 | 25 |
25 // Twiddle tables for the FFT. | 26 // Twiddle tables for the FFT. |
26 this.twiddleCos = new Float32Array(this.N); | 27 this.twiddleCos = new Float32Array(this.N); |
27 this.twiddleSin = new Float32Array(this.N); | 28 this.twiddleSin = new Float32Array(this.N); |
28 | 29 |
29 var omega = -2 * Math.PI / this.N; | 30 let omega = -2 * Math.PI / this.N; |
30 for (var k = 0; k < this.N; ++k) { | 31 for (let k = 0; k < this.N; ++k) { |
31 this.twiddleCos[k] = Math.fround(Math.cos(omega * k)); | 32 this.twiddleCos[k] = Math.fround(Math.cos(omega * k)); |
32 this.twiddleSin[k] = Math.fround(Math.sin(omega * k)); | 33 this.twiddleSin[k] = Math.fround(Math.sin(omega * k)); |
33 } | 34 } |
34 } | 35 } |
35 | 36 |
36 FFT.prototype.FFTException = function (order) { | 37 FFT.prototype.FFTException = |
37 this.value = order; | 38 function(order) { |
38 this.message = "Order must be greater than 1: "; | 39 this.value = order; |
39 this.toString = function () { | 40 this.message = 'Order must be greater than 1: '; |
40 return this.message + this.value; | 41 this.toString = function() { |
41 }; | 42 return this.message + this.value; |
42 } | 43 }; |
43 | 44 } |
44 // Core routine that does one stage of the FFT, implementing all of | 45 |
45 // the butterflies for that stage. | 46 // Core routine that does one stage of the FFT, implementing all of |
46 FFT.prototype.FFTRadix2Core = function (aReal, aImag, bReal, bImag) { | 47 // the butterflies for that stage. |
47 var index = 0; | 48 FFT.prototype.FFTRadix2Core = |
48 for (var k = 0; k < this.numberOfGroups; ++k) { | 49 function(aReal, aImag, bReal, bImag) { |
49 var jfirst = 2 * k * this.pairsInGroup; | 50 let index = 0; |
50 var jlast = jfirst + this.pairsInGroup - 1; | 51 for (let k = 0; k < this.numberOfGroups; ++k) { |
51 var jtwiddle = k * this.pairsInGroup; | 52 let jfirst = 2 * k * this.pairsInGroup; |
52 var wr = this.twiddleCos[jtwiddle]; | 53 let jlast = jfirst + this.pairsInGroup - 1; |
53 var wi = this.twiddleSin[jtwiddle]; | 54 let jtwiddle = k * this.pairsInGroup; |
54 | 55 let wr = this.twiddleCos[jtwiddle]; |
55 for (var j = jfirst; j <= jlast; ++j) { | 56 let wi = this.twiddleSin[jtwiddle]; |
56 // temp = w * a[j + distance] | 57 |
57 var idx = j + this.distance; | 58 for (let j = jfirst; j <= jlast; ++j) { |
58 var tr = wr * aReal[idx] - wi * aImag[idx]; | 59 // temp = w * a[j + distance] |
59 var ti = wr * aImag[idx] + wi * aReal[idx]; | 60 let idx = j + this.distance; |
60 | 61 let tr = wr * aReal[idx] - wi * aImag[idx]; |
61 bReal[index] = aReal[j] + tr; | 62 let ti = wr * aImag[idx] + wi * aReal[idx]; |
62 bImag[index] = aImag[j] + ti; | 63 |
63 bReal[index + this.halfN] = aReal[j] - tr; | 64 bReal[index] = aReal[j] + tr; |
64 bImag[index + this.halfN] = aImag[j] - ti; | 65 bImag[index] = aImag[j] + ti; |
65 ++index; | 66 bReal[index + this.halfN] = aReal[j] - tr; |
66 } | 67 bImag[index + this.halfN] = aImag[j] - ti; |
67 } | 68 ++index; |
68 } | 69 } |
69 | 70 } |
70 // Forward out-of-place complex FFT. | 71 } |
71 // | 72 |
72 // Computes the forward FFT, b, of a complex vector x: | 73 // Forward out-of-place complex FFT. |
73 // | 74 // |
74 // b[k] = sum(x[n] * W^(k*n), n, 0, N - 1), k = 0, 1,..., N-1 | 75 // Computes the forward FFT, b, of a complex vector x: |
75 // | 76 // |
76 // where | 77 // b[k] = sum(x[n] * W^(k*n), n, 0, N - 1), k = 0, 1,..., N-1 |
77 // N = length of x, which must be a power of 2 and | 78 // |
78 // W = exp(-2*i*pi/N) | 79 // where |
79 // | 80 // N = length of x, which must be a power of 2 and |
80 // x = |xr| + i*|xi| | 81 // W = exp(-2*i*pi/N) |
81 // b = |bReal| + i*|bImag| | 82 // |
82 FFT.prototype.fft = function (xr, xi, bReal, bImag) { | 83 // x = |xr| + i*|xi| |
83 this.pairsInGroup = this.halfN; | 84 // b = |bReal| + i*|bImag| |
84 this.numberOfGroups = 1; | 85 FFT.prototype.fft = |
85 this.distance = this.halfN; | 86 function(xr, xi, bReal, bImag) { |
86 this.notSwitchInput = true; | 87 this.pairsInGroup = this.halfN; |
87 | 88 this.numberOfGroups = 1; |
88 // Arrange it so that the last iteration puts the desired output | 89 this.distance = this.halfN; |
89 // in bReal/bImag. | 90 this.notSwitchInput = true; |
90 if (this.order & 1 == 1) { | 91 |
91 this.FFTRadix2Core(xr, xi, bReal, bImag); | 92 // Arrange it so that the last iteration puts the desired output |
92 this.notSwitchInput = !this.notSwitchInput; | 93 // in bReal/bImag. |
| 94 if (this.order & 1 == 1) { |
| 95 this.FFTRadix2Core(xr, xi, bReal, bImag); |
| 96 this.notSwitchInput = !this.notSwitchInput; |
| 97 } else { |
| 98 this.FFTRadix2Core(xr, xi, this.aReal, this.aImag); |
| 99 } |
| 100 |
| 101 this.pairsInGroup >>= 1; |
| 102 this.numberOfGroups <<= 1; |
| 103 this.distance >>= 1; |
| 104 |
| 105 while (this.numberOfGroups < this.N) { |
| 106 if (this.notSwitchInput) { |
| 107 this.FFTRadix2Core(this.aReal, this.aImag, bReal, bImag); |
93 } else { | 108 } else { |
94 this.FFTRadix2Core(xr, xi, this.aReal, this.aImag); | 109 this.FFTRadix2Core(bReal, bImag, this.aReal, this.aImag); |
95 } | 110 } |
96 | 111 |
| 112 this.notSwitchInput = !this.notSwitchInput; |
97 this.pairsInGroup >>= 1; | 113 this.pairsInGroup >>= 1; |
98 this.numberOfGroups <<= 1; | 114 this.numberOfGroups <<= 1; |
99 this.distance >>= 1; | 115 this.distance >>= 1; |
100 | 116 } |
101 while (this.numberOfGroups < this.N) { | 117 } |
102 if (this.notSwitchInput) { | 118 |
103 this.FFTRadix2Core(this.aReal, this.aImag, bReal, bImag); | 119 // Core routine that does one stage of the FFT, implementing all of |
104 } else { | 120 // the butterflies for that stage. This is identical to |
105 this.FFTRadix2Core(bReal, bImag, this.aReal, this.aImag); | 121 // FFTRadix2Core, except the twiddle factor, w, is the conjugate. |
106 } | 122 FFT.prototype.iFFTRadix2Core = |
107 | 123 function(aReal, aImag, bReal, bImag) { |
108 this.notSwitchInput = !this.notSwitchInput; | 124 let index = 0; |
109 this.pairsInGroup >>= 1; | 125 for (let k = 0; k < this.numberOfGroups; ++k) { |
110 this.numberOfGroups <<= 1; | 126 let jfirst = 2 * k * this.pairsInGroup; |
111 this.distance >>= 1; | 127 let jlast = jfirst + this.pairsInGroup - 1; |
112 } | 128 let jtwiddle = k * this.pairsInGroup; |
113 } | 129 let wr = this.twiddleCos[jtwiddle]; |
114 | 130 let wi = -this.twiddleSin[jtwiddle]; |
115 // Core routine that does one stage of the FFT, implementing all of | 131 |
116 // the butterflies for that stage. This is identical to | 132 for (let j = jfirst; j <= jlast; ++j) { |
117 // FFTRadix2Core, except the twiddle factor, w, is the conjugate. | 133 // temp = w * a[j + distance] |
118 FFT.prototype.iFFTRadix2Core = function (aReal, aImag, bReal, bImag) { | 134 let idx = j + this.distance; |
119 var index = 0; | 135 let tr = wr * aReal[idx] - wi * aImag[idx]; |
120 for (var k = 0; k < this.numberOfGroups; ++k) { | 136 let ti = wr * aImag[idx] + wi * aReal[idx]; |
121 var jfirst = 2 * k * this.pairsInGroup; | 137 |
122 var jlast = jfirst + this.pairsInGroup - 1; | 138 bReal[index] = aReal[j] + tr; |
123 var jtwiddle = k * this.pairsInGroup; | 139 bImag[index] = aImag[j] + ti; |
124 var wr = this.twiddleCos[jtwiddle]; | 140 bReal[index + this.halfN] = aReal[j] - tr; |
125 var wi = -this.twiddleSin[jtwiddle]; | 141 bImag[index + this.halfN] = aImag[j] - ti; |
126 | 142 ++index; |
127 for (var j = jfirst; j <= jlast; ++j) { | 143 } |
128 // temp = w * a[j + distance] | 144 } |
129 var idx = j + this.distance; | 145 } |
130 var tr = wr * aReal[idx] - wi * aImag[idx]; | 146 |
131 var ti = wr * aImag[idx] + wi * aReal[idx]; | 147 // Inverse out-of-place complex FFT. |
132 | 148 // |
133 bReal[index] = aReal[j] + tr; | 149 // Computes the inverse FFT, b, of a complex vector x: |
134 bImag[index] = aImag[j] + ti; | 150 // |
135 bReal[index + this.halfN] = aReal[j] - tr; | 151 // b[k] = sum(x[n] * W^(-k*n), n, 0, N - 1), k = 0, 1,..., N-1 |
136 bImag[index + this.halfN] = aImag[j] - ti; | 152 // |
137 ++index; | 153 // where |
138 } | 154 // N = length of x, which must be a power of 2 and |
139 } | 155 // W = exp(-2*i*pi/N) |
140 } | 156 // |
141 | 157 // x = |xr| + i*|xi| |
142 // Inverse out-of-place complex FFT. | 158 // b = |bReal| + i*|bImag| |
143 // | 159 // |
144 // Computes the inverse FFT, b, of a complex vector x: | 160 // Note that ifft(fft(x)) = N * x. To get x, call ifftScale to |
145 // | 161 // scale the output of the ifft appropriately. |
146 // b[k] = sum(x[n] * W^(-k*n), n, 0, N - 1), k = 0, 1,..., N-1 | 162 FFT.prototype.ifft = |
147 // | 163 function(xr, xi, bReal, bImag) { |
148 // where | 164 this.pairsInGroup = this.halfN; |
149 // N = length of x, which must be a power of 2 and | 165 this.numberOfGroups = 1; |
150 // W = exp(-2*i*pi/N) | 166 this.distance = this.halfN; |
151 // | 167 this.notSwitchInput = true; |
152 // x = |xr| + i*|xi| | 168 |
153 // b = |bReal| + i*|bImag| | 169 // Arrange it so that the last iteration puts the desired output |
154 // | 170 // in bReal/bImag. |
155 // Note that ifft(fft(x)) = N * x. To get x, call ifftScale to scale | 171 if (this.order & 1 == 1) { |
156 // the output of the ifft appropriately. | 172 this.iFFTRadix2Core(xr, xi, bReal, bImag); |
157 FFT.prototype.ifft = function (xr, xi, bReal, bImag) { | 173 this.notSwitchInput = !this.notSwitchInput; |
158 this.pairsInGroup = this.halfN; | 174 } else { |
159 this.numberOfGroups = 1; | 175 this.iFFTRadix2Core(xr, xi, this.aReal, this.aImag); |
160 this.distance = this.halfN; | 176 } |
161 this.notSwitchInput = true; | 177 |
162 | 178 this.pairsInGroup >>= 1; |
163 // Arrange it so that the last iteration puts the desired output | 179 this.numberOfGroups <<= 1; |
164 // in bReal/bImag. | 180 this.distance >>= 1; |
165 if (this.order & 1 == 1) { | 181 |
166 this.iFFTRadix2Core(xr, xi, bReal, bImag); | 182 while (this.numberOfGroups < this.N) { |
167 this.notSwitchInput = !this.notSwitchInput; | 183 if (this.notSwitchInput) { |
| 184 this.iFFTRadix2Core(this.aReal, this.aImag, bReal, bImag); |
168 } else { | 185 } else { |
169 this.iFFTRadix2Core(xr, xi, this.aReal, this.aImag); | 186 this.iFFTRadix2Core(bReal, bImag, this.aReal, this.aImag); |
170 } | 187 } |
171 | 188 |
| 189 this.notSwitchInput = !this.notSwitchInput; |
172 this.pairsInGroup >>= 1; | 190 this.pairsInGroup >>= 1; |
173 this.numberOfGroups <<= 1; | 191 this.numberOfGroups <<= 1; |
174 this.distance >>= 1; | 192 this.distance >>= 1; |
175 | 193 } |
176 while (this.numberOfGroups < this.N) { | 194 } |
177 | 195 |
178 if (this.notSwitchInput) { | 196 // |
179 this.iFFTRadix2Core(this.aReal, this.aImag, bReal, bImag); | 197 // Scales the IFFT by 1/N, done in place. |
180 } else { | 198 FFT.prototype.ifftScale = |
181 this.iFFTRadix2Core(bReal, bImag, this.aReal, this.aImag); | 199 function(xr, xi) { |
182 } | 200 let factor = 1 / this.N; |
183 | 201 for (let k = 0; k < this.N; ++k) { |
184 this.notSwitchInput = !this.notSwitchInput; | 202 xr[k] *= factor; |
185 this.pairsInGroup >>= 1; | 203 xi[k] *= factor; |
186 this.numberOfGroups <<= 1; | 204 } |
187 this.distance >>= 1; | 205 } |
188 } | 206 |
189 } | 207 // First stage for the RFFT. Basically the same as |
190 | 208 // FFTRadix2Core, but we assume aImag is 0, and adjust |
191 // | 209 // the code accordingly. |
192 // Scales the IFFT by 1/N, done in place. | 210 FFT.prototype.RFFTRadix2CoreStage1 = |
193 FFT.prototype.ifftScale = function (xr, xi) { | 211 function(aReal, bReal, bImag) { |
194 var factor = 1 / this.N; | 212 let index = 0; |
195 for (var k = 0; k < this.N; ++k) { | 213 for (let k = 0; k < this.numberOfGroups; ++k) { |
196 xr[k] *= factor; | 214 let jfirst = 2 * k * this.pairsInGroup; |
197 xi[k] *= factor; | 215 let jlast = jfirst + this.pairsInGroup - 1; |
198 } | 216 let jtwiddle = k * this.pairsInGroup; |
199 } | 217 let wr = this.twiddleCos[jtwiddle]; |
200 | 218 let wi = this.twiddleSin[jtwiddle]; |
201 // First stage for the RFFT. Basically the same as FFTRadix2Core, but | 219 |
202 // we assume aImag is 0, and adjust the code accordingly. | 220 for (let j = jfirst; j <= jlast; ++j) { |
203 FFT.prototype.RFFTRadix2CoreStage1 = function (aReal, bReal, bImag) { | 221 // temp = w * a[j + distance] |
204 var index = 0; | 222 let idx = j + this.distance; |
205 for (var k = 0; k < this.numberOfGroups; ++k) { | 223 let tr = wr * aReal[idx]; |
206 var jfirst = 2 * k * this.pairsInGroup; | 224 let ti = wi * aReal[idx]; |
207 var jlast = jfirst + this.pairsInGroup - 1; | 225 |
208 var jtwiddle = k * this.pairsInGroup; | 226 bReal[index] = aReal[j] + tr; |
209 var wr = this.twiddleCos[jtwiddle]; | 227 bImag[index] = ti; |
210 var wi = this.twiddleSin[jtwiddle]; | 228 bReal[index + this.halfN] = aReal[j] - tr; |
211 | 229 bImag[index + this.halfN] = -ti; |
212 for (var j = jfirst; j <= jlast; ++j) { | 230 ++index; |
213 // temp = w * a[j + distance] | 231 } |
214 var idx = j + this.distance; | 232 } |
215 var tr = wr * aReal[idx]; | 233 } |
216 var ti = wi * aReal[idx]; | 234 |
217 | 235 // Forward Real FFT. Like fft, but the signal is |
218 bReal[index] = aReal[j] + tr; | 236 // assumed to be real, so the imaginary part is not |
219 bImag[index] = ti; | 237 // supplied. The output, however, is still the same |
220 bReal[index + this.halfN] = aReal[j] - tr; | 238 // and is returned in two arrays. (This could be |
221 bImag[index + this.halfN] = -ti; | 239 // optimized to use less storage, both internally |
222 ++index; | 240 // and for the output, but we don't do that.) |
223 } | 241 FFT.prototype.rfft = function(xr, bReal, bImag) { |
224 } | 242 this.pairsInGroup = this.halfN; |
225 } | 243 this.numberOfGroups = 1; |
226 | 244 this.distance = this.halfN; |
227 // Forward Real FFT. Like fft, but the signal is assumed to be real, | 245 this.notSwitchInput = true; |
228 // so the imaginary part is not supplied. The output, however, is | 246 |
229 // still the same and is returned in two arrays. (This could be | 247 // Arrange it so that the last iteration puts the desired output |
230 // optimized to use less storage, both internally and for the output, | 248 // in bReal/bImag. |
231 // but we don't do that.) | 249 if (this.order & 1 == 1) { |
232 FFT.prototype.rfft = function (xr, bReal, bImag) { | 250 this.RFFTRadix2CoreStage1(xr, bReal, bImag); |
233 this.pairsInGroup = this.halfN; | 251 this.notSwitchInput = !this.notSwitchInput; |
234 this.numberOfGroups = 1; | 252 } else { |
235 this.distance = this.halfN; | 253 this.RFFTRadix2CoreStage1(xr, this.aReal, this.aImag); |
236 this.notSwitchInput = true; | 254 } |
237 | 255 |
238 // Arrange it so that the last iteration puts the desired output | 256 this.pairsInGroup >>= 1; |
239 // in bReal/bImag. | 257 this.numberOfGroups <<= 1; |
240 if (this.order & 1 == 1) { | 258 this.distance >>= 1; |
241 this.RFFTRadix2CoreStage1(xr, bReal, bImag); | 259 |
242 this.notSwitchInput = !this.notSwitchInput; | 260 while (this.numberOfGroups < this.N) { |
| 261 if (this.notSwitchInput) { |
| 262 this.FFTRadix2Core(this.aReal, this.aImag, bReal, bImag); |
243 } else { | 263 } else { |
244 this.RFFTRadix2CoreStage1(xr, this.aReal, this.aImag); | 264 this.FFTRadix2Core(bReal, bImag, this.aReal, this.aImag); |
245 } | 265 } |
246 | 266 |
| 267 this.notSwitchInput = !this.notSwitchInput; |
247 this.pairsInGroup >>= 1; | 268 this.pairsInGroup >>= 1; |
248 this.numberOfGroups <<= 1; | 269 this.numberOfGroups <<= 1; |
249 this.distance >>= 1; | 270 this.distance >>= 1; |
250 | 271 } |
251 while (this.numberOfGroups < this.N) { | 272 } |
252 if (this.notSwitchInput) { | |
253 this.FFTRadix2Core(this.aReal, this.aImag, bReal, bImag); | |
254 } else { | |
255 this.FFTRadix2Core(bReal, bImag, this.aReal, this.aImag); | |
256 } | |
257 | |
258 this.notSwitchInput = !this.notSwitchInput; | |
259 this.pairsInGroup >>= 1; | |
260 this.numberOfGroups <<= 1; | |
261 this.distance >>= 1; | |
262 } | |
263 } | |
OLD | NEW |