| 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 |