30:T1d44,fft_in_place = { /* * Free FFT and convolution (JavaScript) * * Copyright (c) 2017 Project Nayuki. (MIT License) * https://www.nayuki.io/page/free-small-fft-in-multiple-languages * * Permission is hereby granted, free of charge, to any person obtaining a copy of * this software and associated documentation files (the "Software"), to deal in * the Software without restriction, including without limitation the rights to * use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of * the Software, and to permit persons to whom the Software is furnished to do so, * subject to the following conditions: * - The above copyright notice and this permission notice shall be included in * all copies or substantial portions of the Software. * - The Software is provided "as is", without warranty of any kind, express or * implied, including but not limited to the warranties of merchantability, * fitness for a particular purpose and noninfringement. In no event shall the * authors or copyright holders be liable for any claim, damages or other * liability, whether in an action of contract, tort or otherwise, arising from, * out of or in connection with the Software or the use or other dealings in the * Software. */ /* * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. * The vector can have any length. This is a wrapper function. */ function transform(real, imag) { var n = real.length; if (n != imag.length) throw "Mismatched lengths"; if (n == 0) return; else if ((n & (n - 1)) == 0) // Is power of 2 transformRadix2(real, imag); else // More complicated algorithm for arbitrary sizes transformBluestein(real, imag); } /* * Computes the inverse discrete Fourier transform (IDFT) of the given complex vector, storing the result back into the vector. * 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. */ function inverseTransform(real, imag) { transform(imag, real); } /* * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. * The vector's length must be a power of 2. Uses the Cooley-Tukey decimation-in-time radix-2 algorithm. */ function transformRadix2(real, imag) { // Length variables var n = real.length; if (n != imag.length) throw "Mismatched lengths"; if (n == 1) // Trivial transform return; var levels = -1; for (var i = 0; i < 32; i++) { if (1 << i == n) levels = i; // Equal to log2(n) } if (levels == -1) throw "Length is not a power of 2"; // Trigonometric tables var cosTable = new Array(n / 2); var sinTable = new Array(n / 2); for (var i = 0; i < n / 2; i++) { cosTable[i] = Math.cos(2 * Math.PI * i / n); sinTable[i] = Math.sin(2 * Math.PI * i / n); } // Bit-reversed addressing permutation for (var i = 0; i < n; i++) { var j = reverseBits(i, levels); if (j > i) { var temp = real[i]; real[i] = real[j]; real[j] = temp; temp = imag[i]; imag[i] = imag[j]; imag[j] = temp; } } // Cooley-Tukey decimation-in-time radix-2 FFT for (var size = 2; size <= n; size *= 2) { var halfsize = size / 2; var tablestep = n / size; for (var i = 0; i < n; i += size) { for (var j = i, k = 0; j < i + halfsize; j++, k += tablestep) { var l = j + halfsize; var tpre = real[l] * cosTable[k] + imag[l] * sinTable[k]; var tpim = -real[l] * sinTable[k] + imag[l] * cosTable[k]; real[l] = real[j] - tpre; imag[l] = imag[j] - tpim; real[j] += tpre; imag[j] += tpim; } } } // Returns the integer whose value is the reverse of the lowest 'bits' bits of the integer 'x'. function reverseBits(x, bits) { var y = 0; for (var i = 0; i < bits; i++) { y = (y << 1) | (x & 1); x >>>= 1; } return y; } } /* * Computes the discrete Fourier transform (DFT) of the given complex vector, storing the result back into the vector. * The vector can have any length. This requires the convolution function, which in turn requires the radix-2 FFT function. * Uses Bluestein's chirp z-transform algorithm. */ function transformBluestein(real, imag) { // Find a power-of-2 convolution length m such that m >= n * 2 + 1 var n = real.length; if (n != imag.length) throw "Mismatched lengths"; var m = 1; while (m < n * 2 + 1) m *= 2; // Trignometric tables var cosTable = new Array(n); var sinTable = new Array(n); for (var i = 0; i < n; i++) { var j = i * i % (n * 2); // This is more accurate than j = i * i cosTable[i] = Math.cos(Math.PI * j / n); sinTable[i] = Math.sin(Math.PI * j / n); } // Temporary vectors and preprocessing var areal = newArrayOfZeros(m); var aimag = newArrayOfZeros(m); for (var i = 0; i < n; i++) { areal[i] = real[i] * cosTable[i] + imag[i] * sinTable[i]; aimag[i] = -real[i] * sinTable[i] + imag[i] * cosTable[i]; } var breal = newArrayOfZeros(m); var bimag = newArrayOfZeros(m); breal[0] = cosTable[0]; bimag[0] = sinTable[0]; for (var i = 1; i < n; i++) { breal[i] = breal[m - i] = cosTable[i]; bimag[i] = bimag[m - i] = sinTable[i]; } // Convolution var creal = new Array(m); var cimag = new Array(m); convolveComplex(areal, aimag, breal, bimag, creal, cimag); // Postprocessing for (var i = 0; i < n; i++) { real[i] = creal[i] * cosTable[i] + cimag[i] * sinTable[i]; imag[i] = -creal[i] * sinTable[i] + cimag[i] * cosTable[i]; } } /* * Computes the circular convolution of the given real vectors. Each vector's length must be the same. */ function convolveReal(x, y, out) { var n = x.length; if (n != y.length || n != out.length) throw "Mismatched lengths"; convolveComplex(x, newArrayOfZeros(n), y, newArrayOfZeros(n), out, newArrayOfZeros(n)); } /* * Computes the circular convolution of the given complex vectors. Each vector's length must be the same. */ function convolveComplex(xreal, ximag, yreal, yimag, outreal, outimag) { var n = xreal.length; if (n != ximag.length || n != yreal.length || n != yimag.length || n != outreal.length || n != outimag.length) throw "Mismatched lengths"; xreal = xreal.slice(); ximag = ximag.slice(); yreal = yreal.slice(); yimag = yimag.slice(); transform(xreal, ximag); transform(yreal, yimag); for (var i = 0; i < n; i++) { var temp = xreal[i] * yreal[i] - ximag[i] * yimag[i]; ximag[i] = ximag[i] * yreal[i] + xreal[i] * yimag[i]; xreal[i] = temp; } inverseTransform(xreal, ximag); for (var i = 0; i < n; i++) { // Scaling (because this FFT implementation omits it) outreal[i] = xreal[i] / n; outimag[i] = ximag[i] / n; } } function newArrayOfZeros(n) { var result = []; for (var i = 0; i < n; i++) result.push(0); return result; } return transform; }31:[] 2e:{"title":"Function Approximation Using Chebyshev Polynomials","cells":[{"id":0,"value":"Title = md`# Function Approximation Using Chebyshev Polynomials`","pinned":false,"hidden":false,"mode":"ojs"},{"id":450,"value":"{\n const height = 400;\n const svg = d3.select(DOM.svg(width, height));\n \n const heightscale = d3.scaleLinear().domain([-1, 1]).range([height-5, 5]);\n const widthscale = d3.scaleLinear().domain([-1, 1]).range([5, width-5]);\n \n const k = 40;\n const data = d3.range(k, 0, -1).map(i => {\n const v = []; v.length = k; v.fill(0); v[i] = 1; return v;});\n \n const d3line = d3.line().x(d => widthscale(d[0])).y(d => heightscale(d[1]));\n const npts = 400;\n const pts = chebpts(npts);\n \n const curve = svg.append(\"g\").attr(\"fill\", \"none\")\n .selectAll(\"path\")\n .data(data)\n .enter().append(\"path\")\n .attr(\"d\", d => d3line(pts.map(p => [p, evaluate(d, p)])))\n .attr(\"stroke\", (d, i) => d3.interpolatePurples(Math.pow(i/k, 3)));\n \n return svg.node();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":225,"value":"md`A tiny subset of the tools in [Chebfun](http://www.chebfun.org). Just barely enough here to be generally useful so far. For the moment mostly a place for me to park functions for use elsewhere. Eventually would be nice to put up some explanatory interactive doodles.\n\n\\`\\`\\`js\nimport {\n evaluate, sum, cumsum, diff,\n chebpts, vals2coeffs, coeffs2vals\n} from '@jrus/cheb'\n\\`\\`\\`\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":389,"value":"md`\n\nTo understand these, compare to implementations ([1](https://github.com/chebfun/chebfun/tree/development/%40chebtech), [2](https://github.com/chebfun/chebfun/tree/development/%40chebtech2)) in Chebfun, read the [Chebfun guide](http://www.chebfun.org/docs/guide/guide04.html), or look at Mason & Handscomb (2002) ([amzn](https://amzn.com/0849303559), [crc press](https://crcpress.com/p/book/9780849303555)).\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2,"value":"// Clenshaw's method.\nevaluate = function evaluate(coeffs, x) {\n const x2 = 2 * x;\n let d = coeffs.length - 1, b2 = 0, b1 = (d % 2) ? coeffs[d--] : 0;\n for (let i = d; i >= 2; i -= 2) {\n b2 = coeffs[i] + x2 * b1 - b2;\n b1 = coeffs[i-1] + x2 * b2 - b1;\n }\n return coeffs[0] + x * b1 - b2;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":105,"value":"// definite integral from -1 to 1\nsum = function sum(coeffs) {\n let s = 0, i = coeffs.length;\n i = i + i % 2;\n while ((i -= 2) >= 0) s += 2 * coeffs[i] / (1 - i*i);\n return s;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":9,"value":"// indefinite integral\ncumsum = function cumsum(coeffs) {\n let c = coeffs.slice(), i = c.push(0, 0) - 1, k = 0, newcoeffs = [];\n c[0] *= 2;\n while (--i) k = (newcoeffs[i] = (c[i-1] - c[i+1]) / (2*i)) - k;\n newcoeffs[0] = (1 - 2 * coeffs.length % 2) * k;\n return newcoeffs;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":95,"value":"// derivative\ndiff = function diff(coeffs) {\n let newcoeffs = [], i = coeffs.length;\n newcoeffs[i--] = newcoeffs[i] = 0;\n while (i--) newcoeffs[i] = newcoeffs[i+2] + 2*(i+1)*coeffs[i+1];\n newcoeffs[0] *= .5;\n newcoeffs.splice(-2);\n return newcoeffs;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":115,"value":"// Chebyshev nodes of the 2nd kind in [-1, 1]\nchebpts = function chebpts (n) {\n let pts = [], scale = Math.PI * .5 / (n - 1);\n for (let i = n - 1; i >= 0; i--) {\n pts[i] = Math.sin((1 - n + 2*i) * scale);\n }\n return pts;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":392,"value":"md `\n\nThe \\`vals2coeffs\\` and \\`coeffs2vals\\` functions are basically a [discrete cosine transform](https://en.wikipedia.org/wiki/Discrete_cosine_transform) and its inverse. This can be implemented (with some extra computation) using the discrete fourier transform (FFT) by first mirroring the data, then applying the FFT, then truncating the result. Getting the details right in JavaScript makes this look more complicated than it is, because the syntax is not very conducive to array math.\n\nWhenever possible try to use a number of points which is 1 more than a power of 2 (e.g. pad your coefficients out with zeros to the next larger such size), because the FFT algorithm used for general lengths is several times slower.\n\nIf we had the real-valued DCT implemented directly we could make power-of-2 sizes at least 4x faster. We are also leaving at least another order of magnitude in performance on the table vs. a carefully vectorized and cache-optimized native code implementation. Such is browser life.\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":298,"value":"coeffs2vals = function coeffs2vals (coeffs) {\n // for now, only deal with real-valued coefficients.\n coeffs = coeffs.slice();\n const n = coeffs.length - 1;\n if (n <= 0) return coeffs;\n \n for (let i = 1; i < n; i++) {\n coeffs[2*n - i] = coeffs[i] *= .5;\n }\n let imag = []; imag.length = 2*n; imag.fill(0);\n fft_in_place(coeffs, imag);\n \n let values = [];\n for (let i = n; i >=0; i--) {\n values[i] = coeffs[n-i];\n }\n return values;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":339,"value":"vals2coeffs = function vals2coeffs (vals) {\n const n = vals.length - 1;\n if (n <= 0) return vals.slice();\n\n let coeffs = [];\n for (let i = n - 1; i >= 1; i--) {\n coeffs[n + i] = coeffs[n - i] = vals[i];\n }\n coeffs[0] = vals[n]; coeffs[n] = vals[0];\n \n let imag = []; imag.length = 2*n; imag.fill(0);\n fft_in_place(imag, coeffs);\n \n const scale = 1/n;\n for (let i = 1; i <= n; i++) {\n coeffs[i] *= scale;\n }\n coeffs[0] *= scale * .5; coeffs[n] *= scale * .5;\n coeffs.splice(n+1);\n return coeffs;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":596,"value":"// Either truncate the array or pad it with zeros to match the desired length.\nfixlength = function fixlength(array, n) {\n const k = array.length;\n if (k < n) { array.length = n; array.fill(0, k); }\n else { array.splice(n); }\n return array\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":378,"value":"md`\nBasic FFT code from [Project Nayuki](https://www.nayuki.io/page/free-small-fft-in-multiple-languages) under MIT license (eventually perhaps replace this with a purpose-made and optimized DCT implementation):\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":256,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":451,"value":"d3 = require(\"https://d3js.org/d3.v5.min.js\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W31","stdlib":"1","titling":"auto"}