|
| 1 | +/* eslint-disable no-unused-vars */ |
| 2 | +/* eslint-disable brace-style */ |
| 3 | +/* eslint-disable no-array-constructor */ |
| 4 | +/* eslint-disable require-jsdoc */ |
| 5 | +/* eslint-disable max-len */ |
| 6 | +/* |
| 7 | + * Free FFT and convolution (compiled from TypeScript) |
| 8 | + * |
| 9 | + * Copyright (c) 2022 Project Nayuki. (MIT License) |
| 10 | + * https://www.nayuki.io/page/free-small-fft-in-multiple-languages |
| 11 | + * |
| 12 | + * Permission is hereby granted, free of charge, to any person obtaining a copy of |
| 13 | + * this software and associated documentation files (the "Software"), to deal in |
| 14 | + * the Software without restriction, including without limitation the rights to |
| 15 | + * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of |
| 16 | + * the Software, and to permit persons to whom the Software is furnished to do so, |
| 17 | + * subject to the following conditions: |
| 18 | + * - The above copyright notice and this permission notice shall be included in |
| 19 | + * all copies or substantial portions of the Software. |
| 20 | + * - The Software is provided "as is", without warranty of any kind, express or |
| 21 | + * implied, including but not limited to the warranties of merchantability, |
| 22 | + * fitness for a particular purpose and noninfringement. In no event shall the |
| 23 | + * authors or copyright holders be liable for any claim, damages or other |
| 24 | + * liability, whether in an action of contract, tort or otherwise, arising from, |
| 25 | + * out of or in connection with the Software or the use or other dealings in the |
| 26 | + * Software. |
| 27 | + */ |
| 28 | +export const calculateFFT = (real, imag) => { |
| 29 | + const realexp = [...real]; |
| 30 | + const imagexp = [...imag]; |
| 31 | + const inputArray = [realexp, imagexp]; |
| 32 | + |
| 33 | + /* |
| 34 | + * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. |
| 35 | + * The vector can have any length. This is a wrapper function. |
| 36 | + */ |
| 37 | + function transform(real, imag) { |
| 38 | + const output = new Array(); |
| 39 | + const n = real.length; |
| 40 | + if (n != imag.length) throw new RangeError('Mismatched lengths'); |
| 41 | + if (n == 0) return; |
| 42 | + else if ((n & (n - 1)) == 0) { |
| 43 | + // Is power of 2 |
| 44 | + transformRadix2(real, imag); |
| 45 | + |
| 46 | + output.push([real, imag]); |
| 47 | + } |
| 48 | + // More complicated algorithm for arbitrary sizes |
| 49 | + else { |
| 50 | + transformBluestein(real, imag); |
| 51 | + output.push(real, imag); |
| 52 | + } |
| 53 | + |
| 54 | + // // Fake epsilon function |
| 55 | + // for (let i = 0; i < output[0][0].length; i++) { |
| 56 | + // if (Math.abs(output[0][0][i]) < Math.epsilon) { |
| 57 | + // output[0][0][i] = 0; |
| 58 | + // } |
| 59 | + // if (Math.abs(output[0][1][i]) < Math.epsilon) { |
| 60 | + // output[0][1][i] = 0; |
| 61 | + // } |
| 62 | + // } |
| 63 | + return output; |
| 64 | + } |
| 65 | + |
| 66 | + /* |
| 67 | + * Computes the inverse discrete Fourier transform (IDFT) of the given complex vector, storing the result back into the vector. |
| 68 | + * The vector can have any length. This is a wrapper function. This transform does not perform scaling, so the inverse is not a true inverse. |
| 69 | + */ |
| 70 | + function inverseTransform(real, imag) { |
| 71 | + transform(imag, real); |
| 72 | + } |
| 73 | + /* |
| 74 | + * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. |
| 75 | + * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm. |
| 76 | + */ |
| 77 | + function transformRadix2(real, imag) { |
| 78 | + // Length variables |
| 79 | + const n = real.length; |
| 80 | + if (n != imag.length) throw new RangeError('Mismatched lengths'); |
| 81 | + if (n == 1) { |
| 82 | + // Trivial transform |
| 83 | + return; |
| 84 | + } |
| 85 | + let levels = -1; |
| 86 | + for (let i = 0; i < 32; i++) { |
| 87 | + if (1 << i == n) levels = i; // Equal to log2(n) |
| 88 | + } |
| 89 | + if (levels == -1) throw new RangeError('Length is not a power of 2'); |
| 90 | + // Trigonometric tables |
| 91 | + const cosTable = new Array(n / 2); |
| 92 | + const sinTable = new Array(n / 2); |
| 93 | + |
| 94 | + for (let i = 0; i < n / 2; i++) { |
| 95 | + cosTable[i] = Math.cos((2 * Math.PI * i) / n); |
| 96 | + sinTable[i] = Math.sin((2 * Math.PI * i) / n); |
| 97 | + } |
| 98 | + // Bit-reversed addressing permutation |
| 99 | + for (let i = 0; i < n; i++) { |
| 100 | + const j = reverseBits(i, levels); |
| 101 | + |
| 102 | + if (j > i) { |
| 103 | + let temp = real[i]; |
| 104 | + real[i] = real[j]; |
| 105 | + real[j] = temp; |
| 106 | + temp = imag[i]; |
| 107 | + imag[i] = imag[j]; |
| 108 | + imag[j] = temp; |
| 109 | + } |
| 110 | + } |
| 111 | + // Cooley-Tukey decimation-in-time radix-2 FFT |
| 112 | + for (let size = 2; size <= n; size *= 2) { |
| 113 | + const halfsize = size / 2; |
| 114 | + const tablestep = n / size; |
| 115 | + for (let i = 0; i < n; i += size) { |
| 116 | + for (let j = i, k = 0; j < i + halfsize; j++, k += tablestep) { |
| 117 | + const l = j + halfsize; |
| 118 | + const tpre = real[l] * cosTable[k] + imag[l] * sinTable[k]; |
| 119 | + const tpim = -real[l] * sinTable[k] + imag[l] * cosTable[k]; |
| 120 | + real[l] = real[j] - tpre; |
| 121 | + imag[l] = imag[j] - tpim; |
| 122 | + real[j] += tpre; |
| 123 | + imag[j] += tpim; |
| 124 | + } |
| 125 | + } |
| 126 | + } |
| 127 | + // Returns the integer whose value is the reverse of the lowest 'width' bits of the integer 'val'. |
| 128 | + function reverseBits(val, width) { |
| 129 | + let result = 0; |
| 130 | + for (let i = 0; i < width; i++) { |
| 131 | + result = (result << 1) | (val & 1); |
| 132 | + val >>>= 1; |
| 133 | + } |
| 134 | + |
| 135 | + return result; |
| 136 | + } |
| 137 | + } |
| 138 | + /* |
| 139 | + * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. |
| 140 | + * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function. |
| 141 | + * Uses Bluestein's chirp z-transform algorithm. |
| 142 | + */ |
| 143 | + function transformBluestein(real, imag) { |
| 144 | + // Find a power-of-2 convolution length m such that m >= n * 2 + 1 |
| 145 | + const n = real.length; |
| 146 | + if (n != imag.length) throw new RangeError('Mismatched lengths'); |
| 147 | + let m = 1; |
| 148 | + while (m < n * 2 + 1) m *= 2; |
| 149 | + // Trigonometric tables |
| 150 | + const cosTable = new Array(n); |
| 151 | + const sinTable = new Array(n); |
| 152 | + for (let i = 0; i < n; i++) { |
| 153 | + const j = (i * i) % (n * 2); // This is more accurate than j = i * i |
| 154 | + cosTable[i] = Math.cos((Math.PI * j) / n); |
| 155 | + sinTable[i] = Math.sin((Math.PI * j) / n); |
| 156 | + } |
| 157 | + // Temporary vectors and preprocessing |
| 158 | + const areal = newArrayOfZeros(m); |
| 159 | + const aimag = newArrayOfZeros(m); |
| 160 | + for (let i = 0; i < n; i++) { |
| 161 | + areal[i] = real[i] * cosTable[i] + imag[i] * sinTable[i]; |
| 162 | + aimag[i] = -real[i] * sinTable[i] + imag[i] * cosTable[i]; |
| 163 | + } |
| 164 | + const breal = newArrayOfZeros(m); |
| 165 | + const bimag = newArrayOfZeros(m); |
| 166 | + breal[0] = cosTable[0]; |
| 167 | + bimag[0] = sinTable[0]; |
| 168 | + for (let i = 1; i < n; i++) { |
| 169 | + breal[i] = breal[m - i] = cosTable[i]; |
| 170 | + bimag[i] = bimag[m - i] = sinTable[i]; |
| 171 | + } |
| 172 | + // Convolution |
| 173 | + const creal = new Array(m); |
| 174 | + const cimag = new Array(m); |
| 175 | + convolveComplex(areal, aimag, breal, bimag, creal, cimag); |
| 176 | + // Postprocessing |
| 177 | + for (let i = 0; i < n; i++) { |
| 178 | + real[i] = creal[i] * cosTable[i] + cimag[i] * sinTable[i]; |
| 179 | + imag[i] = -creal[i] * sinTable[i] + cimag[i] * cosTable[i]; |
| 180 | + } |
| 181 | + } |
| 182 | + /* |
| 183 | + * Computes the circular convolution of the given real vectors. Each vector's length must be the same. |
| 184 | + */ |
| 185 | + function convolveReal(xvec, yvec, outvec) { |
| 186 | + const n = xvec.length; |
| 187 | + if (n != yvec.length || n != outvec.length) throw new RangeError('Mismatched lengths'); |
| 188 | + convolveComplex(xvec, newArrayOfZeros(n), yvec, newArrayOfZeros(n), outvec, newArrayOfZeros(n)); |
| 189 | + } |
| 190 | + /* |
| 191 | + * Computes the circular convolution of the given complex vectors. Each vector's length must be the same. |
| 192 | + */ |
| 193 | + function convolveComplex(xreal, ximag, yreal, yimag, outreal, outimag) { |
| 194 | + const n = xreal.length; |
| 195 | + if (n != ximag.length || n != yreal.length || n != yimag.length || n != outreal.length || n != outimag.length) throw new RangeError('Missmatched lengths'); |
| 196 | + xreal = xreal.slice(); |
| 197 | + ximag = ximag.slice(); |
| 198 | + yreal = yreal.slice(); |
| 199 | + yimag = yimag.slice(); |
| 200 | + transform(xreal, ximag); |
| 201 | + transform(yreal, yimag); |
| 202 | + for (let i = 0; i < n; i++) { |
| 203 | + const temp = xreal[i] * yreal[i] - ximag[i] * yimag[i]; |
| 204 | + ximag[i] = ximag[i] * yreal[i] + xreal[i] * yimag[i]; |
| 205 | + xreal[i] = temp; |
| 206 | + } |
| 207 | + inverseTransform(xreal, ximag); |
| 208 | + for (let i = 0; i < n; i++) { |
| 209 | + // Scaling (because this FFT implementation omits it) |
| 210 | + outreal[i] = xreal[i] / n; |
| 211 | + outimag[i] = ximag[i] / n; |
| 212 | + } |
| 213 | + } |
| 214 | + function newArrayOfZeros(n) { |
| 215 | + const result = []; |
| 216 | + for (let i = 0; i < n; i++) result.push(0); |
| 217 | + return result; |
| 218 | + } |
| 219 | + |
| 220 | + return transform(inputArray[0], inputArray[1]); |
| 221 | +}; |
0 commit comments