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