From d2a3c99894e200a6d8ed96d6430f7d6b0a9e552e Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Fri, 8 May 2026 21:27:49 +0530 Subject: [PATCH 1/5] feat: copy older implementation, add radf3 --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../fft/base/fftpack/rfftf/lib/index.js | 49 +++ .../fft/base/fftpack/rfftf/lib/main.js | 124 ++++++ .../fft/base/fftpack/rfftf/lib/radf2.js | 366 ++++++++++++++++ .../fft/base/fftpack/rfftf/lib/radf3.js | 397 ++++++++++++++++++ 4 files changed, 936 insertions(+) create mode 100644 lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/index.js create mode 100644 lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/main.js create mode 100644 lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js create mode 100644 lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/index.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/index.js new file mode 100644 index 000000000000..d4e7f52daa0a --- /dev/null +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/index.js @@ -0,0 +1,49 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +*/ + +'use strict'; + +/** +* Compute the forward real-valued fast Fourier transform (FFT). +* +* @module @stdlib/fft/base/fftpack/rfftf +* +* @example +* var Float64Array = require( '@stdlib/array/float64' ); +* var rffti = require( '@stdlib/fft/base/fftpack/rffti' ); +* var rfftf = require( '@stdlib/fft/base/fftpack/rfftf' ); +* +* var N = 4; +* +* var workspace = new Float64Array( ( 2*N ) + 34 ); +* rffti( N, workspace, 1, 0 ); +* +* var r = new Float64Array( [ 1.0, 2.0, 3.0, 4.0 ] ); +* +* rfftf( N, r, 1, 0, workspace, 1, 0 ); +* // r => [ 10.0, -2.0, 2.0, -2.0 ] +*/ + +// MODULES // + +var main = require( './main.js' ); + + +// EXPORTS // + +module.exports = main; diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/main.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/main.js new file mode 100644 index 000000000000..0e54338933a1 --- /dev/null +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/main.js @@ -0,0 +1,124 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +* +* +* ## Notice +* +* The original C code and copyright notice are from the [PFFFT library]{@link https://github.com/marton78/pffft/blob/master/fftpack.c}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (c) 2004 the University Corporation for Atmospheric +* Research ("UCAR"). All rights reserved. Developed by NCAR's +* Computational and Information Systems Laboratory, UCAR, +* www.cisl.ucar.edu. +* +* Redistribution and use of the Software in source and binary forms, +* with or without modification, is permitted provided that the +* following conditions are met: +* +* - Neither the names of NCAR's Computational and Information Systems +* Laboratory, the University Corporation for Atmospheric Research, +* nor the names of its sponsors or contributors may be used to +* endorse or promote products derived from this Software without +* specific prior written permission. +* +* - Redistributions of source code must retain the above copyright +* notices, this list of conditions, and the disclaimer below. +* +* - Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions, and the disclaimer below in the +* documentation and/or other materials provided with the +* distribution. +* +* THIS 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 CONTRIBUTORS OR COPYRIGHT +* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL, +* EXEMPLARY, OR CONSEQUENTIAL 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 WITH THE +* SOFTWARE. +* ``` +*/ + +'use strict'; + +// MODULES // + +var isNonNegativeInteger = require( '@stdlib/assert/is-nonnegative-integer' ); +var isFloat64Array = require( '@stdlib/assert/is-float64array' ); +var isInteger = require( '@stdlib/assert/is-integer' ); +var rfftf1 = require( './rfftf1.js' ); + + +// MAIN // + +/** +* Computes the forward real-valued fast Fourier transform (FFT). +* +* @param {NonNegativeInteger} N - length of the sequence to transform +* @param {Float64Array} r - input array containing the sequence to transform +* @param {integer} strideR - stride length for `r` +* @param {NonNegativeInteger} offsetR - starting index for `r` +* @param {Float64Array} workspace - workspace array containing pre-computed values +* @param {integer} strideW - stride length for `workspace` +* @param {NonNegativeInteger} offsetW - starting index for `workspace` +* @returns {void} +* +* @example +* var Float64Array = require( '@stdlib/array/float64' ); +* var rffti = require( '@stdlib/fft/base/fftpack/rffti' ); +* +* var N = 4; +* +* var workspace = new Float64Array( ( 2*N ) + 34 ); +* rffti( N, workspace, 1, 0 ); +* +* var r = new Float64Array( [ 1.0, 2.0, 3.0, 4.0 ] ); +* +* rfftf( N, r, 1, 0, workspace, 1, 0 ); +* // r => [ 10.0, -2.0, 2.0, -2.0 ] +*/ +function rfftf( N, r, strideR, offsetR, workspace, strideW, offsetW ) { + var offsetT; + var offsetF; + + if ( !isNonNegativeInteger( N ) || !isFloat64Array( r ) || + !isInteger( strideR ) || !isNonNegativeInteger( offsetR ) || + !isFloat64Array( workspace ) || !isInteger( strideW ) || + !isNonNegativeInteger( offsetW ) || + workspace.length < offsetW + ( ( ( 2*N ) + 34 ) * strideW ) ) { + return; + } + + // When a sub-sequence is a single data point, the FFT is the identity, so no transformation necessary... + if ( N === 1 ) { + return; + } + + // Resolve the starting indices for storing twiddle factors and factorization results: + offsetT = offsetW + (N * strideW); // index offset for twiddle factors + offsetF = offsetT + (N * strideW); // index offset for factors describing the sub-transforms + + rfftf1( N, r, offsetR, workspace, offsetW, workspace, offsetT, workspace, offsetF ); // eslint-disable-line max-len +} + + +// EXPORTS // + +module.exports = rfftf; diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js new file mode 100644 index 000000000000..0d0af15bbf79 --- /dev/null +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js @@ -0,0 +1,366 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +* +* +* ## Notice +* +* The original C code and copyright notice are from the [PFFFT library]{@link https://github.com/marton78/pffft/blob/0b4ee12c4ba45a4a8e567550c16d96d1679f50ce/src/fftpack.c}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (c) 2004 the University Corporation for Atmospheric +* Research ("UCAR"). All rights reserved. Developed by NCAR's +* Computational and Information Systems Laboratory, UCAR, +* www.cisl.ucar.edu. +* +* Redistribution and use of the Software in source and binary forms, +* with or without modification, is permitted provided that the +* following conditions are met: +* +* - Neither the names of NCAR's Computational and Information Systems +* Laboratory, the University Corporation for Atmospheric Research, +* nor the names of its sponsors or contributors may be used to +* endorse or promote products derived from this Software without +* specific prior written permission. +* +* - Redistributions of source code must retain the above copyright +* notices, this list of conditions, and the disclaimer below. +* +* - Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions, and the disclaimer below in the +* documentation and/or other materials provided with the +* distribution. +* +* THIS 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 CONTRIBUTORS OR COPYRIGHT +* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL, +* EXEMPLARY, OR CONSEQUENTIAL 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 WITH THE +* SOFTWARE. +* ``` +*/ + +/* eslint-disable max-len */ + +'use strict'; + +// FUNCTIONS // + +/** +* Resolves an index into the input array. +* +* ## Notes +* +* In a forward real FFT, the previous stage writes its results as two "rows" per sub-sequence. +* +* Thus, when reading from an input array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having two "rows" (`even + odd` and `even - odd`, respectively) and where each "row" is arranged as `M*L` contiguous elements corresponding to interleaved real and imaginary components. +* +* Accordingly, the following is a logical view of an input array (zero-based indexing) which contains `L = 3` transforms and in which each sub-sequence has length `M = 4`: +* +* ```text +* │ k = 0 k = 1 k = 2 +* │ ──────────────────────────────────────────────────────────────────────────→ k +* j = 0 (even+odd) │ cc(0,0,0) ... cc(3,0,0) cc(0,1,0) ... cc(3,1,0) cc(0,2,0) ... cc(3,2,0) +* │ +* j = 1 (even-odd) │ cc(0,0,1) ... cc(3,0,1) cc(0,1,1) ... cc(3,1,1) cc(0,2,1) ... cc(3,2,1) +* └───────────────────────────────────────────────────────────────────────────→ i +* ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 M-1 0 M-1 0 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short sub-sequence corresponding to either the `even + odd` or `even - odd` row. +* - `j` selects between the `even + odd` and `even - odd` row. +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | cc(0,0,0)...cc(3,0,0) ... cc(0,2,0)...cc(3,2,0) | cc(0,0,1)...cc(3,0,1) ... cc(0,2,1)...cc(3,2,1) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 LM LM-1 (L+1)M (L+1)M-1 (2L-1)M 2LM-1 +* ``` +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} j - input row +* @param {NonNegativeInteger} L - number of sub-sequences +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the input array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = iptr( 0, 0, 0, L, M, stride, offset ); +* // returns 0 +* +* idx = iptr( 1, 0, 0, L, M, stride, offset ); +* // returns 1 +* +* idx = iptr( M-1, 0, 0, L, M, stride, offset ); +* // returns 3 +* +* idx = iptr( 0, 1, 0, L, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = iptr( M-1, L-1, 1, L, M, stride, offset ); +* // returns 23 +*/ +function iptr( i, k, j, L, M, stride, offset ) { + var n = i + ( ( k+(j*L) ) * M ); + return ( n*stride ) + offset; +} + +/** +* Resolves an index into the output array. +* +* ## Notes +* +* When writing to an output array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having two "columns" corresponding to the even and odd half-spectrum of a folded complex vector (with real and imaginary parts of each half-spectrum interleaved along each sub-sequence) and where each "column" has `M` elements. +* +* Accordingly, the following is a logical view of an output array (zero-based indexing) which contains `L = 3` transforms and in which each "column" sub-sequence has length `M = 4`: +* +* ```text +* j = 0 ("even" column) j = 1 ("odd" column) +* k = 0 ─┬───────────────────────────────────────┬───────────────────────────────────────┐ +* │ out(0,0,0) out(1,0,0) ... out(3,0,0) │ out(0,1,0) out(1,1,0) ... out(3,1,0) │ +* └───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 1 ─┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,1) out(1,0,1) ... out(3,0,1) │ out(0,1,1) out(1,1,1) ... out(3,1,1) │ +* └───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 2 ─┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,2) out(1,0,2) ... out(3,0,2) │ out(0,1,2) out(1,1,2) ... out(3,1,2) │ +* └───────────────────────────────────────┴───────────────────────────────────────┘ +* ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 1 M-1 0 1 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short "column" sub-sequence. +* - `j` selects between the even (0) and odd (1) column. +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | out(0,0,0)...out(3,0,0) | out(0,1,0)...out(3,1,0) | out(0,0,1)...out(3,0,1) | ... | out(0,1,2)...out(3,1,2) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 M 2M-1 2M 3M-1 (2L-1)M 2LM-1 +* ``` +* +* As may be observed, when resolving an index in the output array, the `j` and `k` dimensions are swapped relative index resolution in the input array. This stems from `radf2` being only one stage in a multi-stage driver which alternates between using `cc` and `out` as workspace buffers. After each stage, the next stage reads what the previous stage wrote. +* +* Each stage expects a transpose, and, in order to avoid explicit transposition between the stages, we swap the last two logical dimensions while still maintaining cache locality within the inner loop logical dimension, as indexed by `i`. +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} j - index specifying which of the two complex halves we are in (either `0` or `1`) +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = optr( 0, 0, 0, M, stride, offset ); +* // returns 0 +* +* idx = optr( 1, 0, 0, M, stride, offset ); +* // returns 1 +* +* idx = optr( M-1, 0, 0, M, stride, offset ); +* // returns 3 +* +* idx = optr( 0, 1, 0, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = optr( M-1, 1, L-1, M, stride, offset ); +* // returns 23 +*/ +function optr( i, j, k, M, stride, offset ) { + var n = i + ( ( j+(k*2) ) * M ); + return ( n*stride ) + offset; +} + + +// MAIN // + +/** +* Performs one radix-2 stage within a forward Fourier transform for a real-valued sequence. +* +* @private +* @param {NonNegativeInteger} M - number of elements in each sub-sequence to be transformed +* @param {NonNegativeInteger} L - number of sub-sequences to be transformed +* @param {Float64Array} cc - input array containing the sub-sequences to be transformed +* @param {integer} sc - stride length for `cc` +* @param {NonNegativeInteger} oc - index offset for `cc` +* @param {Float64Array} out - output array containing transformed sub-sequences +* @param {integer} so - stride length for `out` +* @param {NonNegativeInteger} oo - index offset for `out` +* @param {Float64Array} twiddles - array containing twiddle factors +* @param {integer} st - stride length for `twiddles` +* @param {NonNegativeInteger} ot - index offset for `twiddles` +* @returns {void} +*/ +function radf2( M, L, cc, sc, oc, out, so, oo, twiddles, st, ot ) { // eslint-disable-line max-params + var tr2; + var ti2; + var ip1; + var ip2; + var it1; + var it2; + var io; + var im; + var i; + var k; + + /* + * First, perform the core butterfly for each sub-sequence being transformed. + * + * In the following loop, we handle harmonic `n = 0` for every transform. As described for `iptr`, the input array is interpreted as a three-dimensional array, containing two "rows" per sub-sequence. + * + * row0 = even + odd (j = 0) + * row1 = even - odd (j = 1) + * + * For a radix-2 forward FFT of a real-valued signal `x` and `n = 0`, + * + * x[0] = row0 + row1 = 2⋅even + * x[M] = row0 - row1 = 2⋅odd + * + * Because `W_2^0 = 1`, no twiddle multiplication is necessary. + */ + for ( k = 0; k < L; k++ ) { + ip1 = iptr( 0, k, 0, L, M, sc, oc ); // real part row 0 + ip2 = iptr( 0, k, 1, L, M, sc, oc ); // real part row 1 + + io = optr( 0, 0, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + cc[ ip2 ]; // even + odd + + io = optr( M-1, 1, k, M, so, oo ); + out[ io ] = cc[ ip1 ] - cc[ ip2 ]; // even - odd + } + // When the number of elements in a sub-sequence is less than `2`, there is nothing more to do, as the above butterfly produced the full result... + if ( M < 2 ) { + return; + } + /* + * Next, apply the general case where we need to loop through the non-trivial harmonics. + * + * For each harmonic `n = 1, ..., M/2-1`, we need to + * + * - recover even/odd spectra from the two input rows. + * - multiply the odd part by the twiddle factor `W_n = cos(θ) - j sin(θ)`. + * - form the following + * + * x[2n] = E_n + W_n⋅O_n => column 0 + * x[2n+1] = E_n - W_n⋅O_n => column 1 + * + * The mirror index `im = M - i` selects the conjugate-symmetric partner, thus allowing the routine to read each symmetry pair only once. + */ + if ( M >= 3 ) { + // Loop over each sub-sequence to be transformed... + for ( k = 0; k < L; k++ ) { + // Loop over the elements in each sub-sequence... + for ( i = 2; i < M; i += 2 ) { + im = M - i; // "mirror" index + + // Resolve twiddle factor indices: + it1 = ( (i-2)*st ) + ot; // cos(θ) + it2 = ( (i-1)*st ) + ot; // sin(θ) + + // Load the `even-odd` row... + ip1 = iptr( i-1, k, 1, L, M, sc, oc ); // c = Re(row1) + ip2 = iptr( i, k, 1, L, M, sc, oc ); // d = Im(row1) + + // eslint-disable-next-line stdlib/capitalized-comments, stdlib/empty-line-before-comment + // tmp = W_n ⋅ (c + j⋅d) + tr2 = ( twiddles[ it1 ] * cc[ ip1 ] ) + ( twiddles[ it2 ] * cc[ ip2 ] ); // Re(tmp) + ti2 = ( twiddles[ it1 ] * cc[ ip2 ] ) - ( twiddles[ it2 ] * cc[ ip1 ] ); // Im(tmp) + + // Load the `even+odd` row... + ip1 = iptr( i, k, 0, L, M, sc, oc ); // b = Im(row0) + io = optr( i, 0, k, M, so, oo ); // Im(x[2n]) + out[ io ] = cc[ ip1 ] + ti2; + + io = optr( im, 1, k, M, so, oo ); // Im(x[2n+1]) + out[ io ] = ti2 - cc[ ip1 ]; + + ip1 = iptr( i-1, k, 0, L, M, sc, oc ); // a = Re(row0) + io = optr( i-1, 0, k, M, so, oo ); // Re(x[2n]) + out[ io ] = cc[ ip1 ] + tr2; + + io = optr( im-1, 1, k, M, so, oo ); // Re(x[2n+1]) + out[ io ] = cc[ ip1 ] - tr2; + } + } + // When `M` is odd, there is no Nyquist pair to process, so we do not need to perform any further transformations... + if ( M%2 === 1 ) { + return; + } + } + /* + * Lastly, handle the Nyquist frequency where `i = M-1` (i.e., the last element of each sub-sequence). + * + * When `M` is even, the Nyquist index is `i = M/2`. In this stage, we've stored that element at the end of each sub-sequence (i.e., `i = M-1` in the packed layout). + * + * At this point, the cosine term is -1 and the sine term is 0, so the twiddle multiplication collapses to simple addition/subtraction. + * + * In this case, + * + * W_n⋅O_n = ±Re(O_n) + * + * where + * + * row0 = Re(E_n) + * row1 = -Re(O_n) + */ + for ( k = 0; k < L; k++ ) { + ip1 = iptr( M-1, k, 1, L, M, sc, oc ); // -Re(O_n) + io = optr( 0, 1, k, M, so, oo ); + out[ io ] = -cc[ ip1 ]; // -(-Re(O_n)) = Re(O_n) + + ip1 = iptr( M-1, k, 0, L, M, sc, oc ); // Re(E_n) + io = optr( M-1, 0, k, M, so, oo ); + out[ io ] = cc[ ip1 ]; + } +} + + +// EXPORTS // + +module.exports = radf2; diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js new file mode 100644 index 000000000000..2fced4b2c683 --- /dev/null +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js @@ -0,0 +1,397 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +* +* +* ## Notice +* +* The original C code and copyright notice are from the [PFFFT library]{@link https://github.com/marton78/pffft/blob/0b4ee12c4ba45a4a8e567550c16d96d1679f50ce/src/fftpack.c}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (c) 2004 the University Corporation for Atmospheric +* Research ("UCAR"). All rights reserved. Developed by NCAR's +* Computational and Information Systems Laboratory, UCAR, +* www.cisl.ucar.edu. +* +* Redistribution and use of the Software in source and binary forms, +* with or without modification, is permitted provided that the +* following conditions are met: +* +* - Neither the names of NCAR's Computational and Information Systems +* Laboratory, the University Corporation for Atmospheric Research, +* nor the names of its sponsors or contributors may be used to +* endorse or promote products derived from this Software without +* specific prior written permission. +* +* - Redistributions of source code must retain the above copyright +* notices, this list of conditions, and the disclaimer below. +* +* - Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions, and the disclaimer below in the +* documentation and/or other materials provided with the +* distribution. +* +* THIS 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 CONTRIBUTORS OR COPYRIGHT +* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL, +* EXEMPLARY, OR CONSEQUENTIAL 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 WITH THE +* SOFTWARE. +* ``` +*/ + +/* eslint-disable max-len */ + +'use strict'; + +// VARIABLES // + +var TAUR = -0.5; // cos( 2π/3 ) +var TAUI = 0.866025403784439; // sin( 2π/3 ) + + +// FUNCTIONS // + +/** +* Resolves an index into the input array. +* +* ## Notes +* +* In a forward real FFT, the previous stage writes its results as three "rows" per sub-sequence. +* +* Thus, when reading from an input array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having three "rows" (components 0, 1, and 2) and where each "row" is arranged as `M*L` contiguous elements corresponding to interleaved real and imaginary components. +* +* Accordingly, the following is a logical view of an input array (zero-based indexing) which contains `L = 3` transforms and in which each sub-sequence has length `M = 4`: +* +* ```text +* │ k = 0 k = 1 k = 2 +* │ ──────────────────────────────────────────────────────────────────────────→ k +* j = 0 │ cc(0,0,0) ... cc(3,0,0) cc(0,1,0) ... cc(3,1,0) cc(0,2,0) ... cc(3,2,0) +* │ +* j = 1 │ cc(0,0,1) ... cc(3,0,1) cc(0,1,1) ... cc(3,1,1) cc(0,2,1) ... cc(3,2,1) +* │ +* j = 2 │ cc(0,0,2) ... cc(3,0,2) cc(0,1,2) ... cc(3,1,2) cc(0,2,2) ... cc(3,2,2) +* └───────────────────────────────────────────────────────────────────────────→ i +* ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 M-1 0 M-1 0 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short sub-sequence corresponding to one of the three component rows. +* - `j` selects between the three component rows (0, 1, or 2). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | cc(0,0,0)...cc(3,0,0) ... cc(0,2,0)...cc(3,2,0) | cc(0,0,1)...cc(3,0,1) ... cc(0,2,1)...cc(3,2,1) | cc(0,0,2)...cc(3,0,2) ... cc(0,2,2)...cc(3,2,2) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 LM LM-1 (L+1)M (L+1)M-1 (2L-1)M 2LM-1 2LM (2L+1)M-1 (3L-1)M 3LM-1 +* ``` +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} j - input row +* @param {NonNegativeInteger} L - number of sub-sequences +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the input array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = iptr( 0, 0, 0, L, M, stride, offset ); +* // returns 0 +* +* idx = iptr( 1, 0, 0, L, M, stride, offset ); +* // returns 1 +* +* idx = iptr( M-1, 0, 0, L, M, stride, offset ); +* // returns 3 +* +* idx = iptr( 0, 1, 0, L, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = iptr( M-1, L-1, 1, L, M, stride, offset ); +* // returns 23 +*/ +function iptr( i, k, j, L, M, stride, offset ) { + var n = i + ( ( k+(j*L) ) * M ); + return ( n*stride ) + offset; +} + +/** +* Resolves an index into the output array. +* +* ## Notes +* +* When writing to an output array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having three "columns" corresponding to the three components of a radix-3 stage (with real and imaginary parts of each component interleaved along each sub-sequence) and where each "column" has `M` elements. +* +* Accordingly, the following is a logical view of an output array (zero-based indexing) which contains `L = 3` transforms and in which each "column" sub-sequence has length `M = 4`: +* +* ```text +* j = 0 (component 0) j = 1 (component 1) j = 2 (component 2) +* k = 0 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┐ +* │ out(0,0,0) out(1,0,0) ... out(3,0,0) │ out(0,1,0) out(1,1,0) ... out(3,1,0) │ out(0,2,0) out(1,2,0) ... out(3,2,0) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 1 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,1) out(1,0,1) ... out(3,0,1) │ out(0,1,1) out(1,1,1) ... out(3,1,1) │ out(0,2,1) out(1,2,1) ... out(3,2,1) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 2 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,2) out(1,0,2) ... out(3,0,2) │ out(0,1,2) out(1,1,2) ... out(3,1,2) │ out(0,2,2) out(1,2,2) ... out(3,2,2) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┘ +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 1 M-1 0 1 M-1 0 1 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short "column" sub-sequence. +* - `j` selects between one of the three components (0, 1, or 2). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | out(0,0,0)...out(3,0,0) | out(0,1,0)...out(3,1,0) | out(0,2,0)...out(3,2,0) | out(0,0,1)...out(3,0,1) | ... | out(0,2,2)...out(3,2,2) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 M 2M-1 2M 3M-1 3M 4M-1 (3L-1)M 3LM-1 +* ``` +* +* As may be observed, when resolving an index in the output array, the `j` and `k` dimensions are swapped relative index resolution in the input array. This stems from `radf3` being only one stage in a multi-stage driver which alternates between using `cc` and `out` as workspace buffers. After each stage, the next stage reads what the previous stage wrote. +* +* Each stage expects a transpose, and, in order to avoid explicit transposition between the stages, we swap the last two logical dimensions while still maintaining cache locality within the inner loop logical dimension, as indexed by `i`. +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} j - index specifying which of the three complex components we are in (0, 1, or 2) +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = optr( 0, 0, 0, M, stride, offset ); +* // returns 0 +* +* idx = optr( 1, 0, 0, M, stride, offset ); +* // returns 1 +* +* idx = optr( M-1, 0, 0, M, stride, offset ); +* // returns 3 +* +* idx = optr( 0, 1, 0, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = optr( M-1, 2, L-1, M, stride, offset ); +* // returns 35 +*/ +function optr( i, j, k, M, stride, offset ) { + var n = i + ( ( j+(k*3) ) * M ); + return ( n*stride ) + offset; +} + + +// MAIN // + +/** +* Performs one radix-3 stage within a forward Fourier transform for a real-valued sequence. +* +* @private +* @param {NonNegativeInteger} M - number of elements in each sub-sequence to be transformed +* @param {NonNegativeInteger} L - number of sub-sequences to be transformed +* @param {Float64Array} cc - input array containing the sub-sequences to be transformed +* @param {integer} sc - stride length for `cc` +* @param {NonNegativeInteger} oc - index offset for `cc` +* @param {Float64Array} out - output array containing transformed sub-sequences +* @param {integer} so - stride length for `out` +* @param {NonNegativeInteger} oo - index offset for `out` +* @param {Float64Array} twiddles1 - first array containing twiddle factors +* @param {integer} st1 - stride length for `twiddles1` +* @param {NonNegativeInteger} ot1 - index offset for `twiddles1` +* @param {Float64Array} twiddles2 - second array containing twiddle factors +* @param {integer} st2 - stride length for `twiddles2` +* @param {NonNegativeInteger} ot2 - index offset for `twiddles2` +* @returns {void} +*/ +function radf3( M, L, cc, sc, oc, out, so, oo, twiddles1, st1, ot1, twiddles2, st2, ot2 ) { // eslint-disable-line max-params + var it11; + var it12; + var it21; + var it22; + var tr2; + var ti2; + var tr3; + var ti3; + var cr2; + var ci2; + var dr2; + var di2; + var dr3; + var di3; + var ip1; + var ip2; + var ip3; + var io; + var im; + var i; + var k; + + /* + * First, perform the core butterfly for each sub-sequence being transformed. + * + * In the following loop, we handle harmonic `n = 0` for every transform. As described for `iptr`, the input array is interpreted as a three-dimensional array, containing three "rows" per sub-sequence. + * + * row0 = input component 0 (j = 0) + * row1 = input component 1 (j = 1) + * row2 = input component 2 (j = 2) + * + * For a radix-3 forward FFT of a real-valued signal `x` and `n = 0`, + * + * x[0] = row0 + row1 + row2 + * x[2M] = taui * ( row2-row1 ) + * x[M] = row0 + ( taur * ( row1+row2 ) ) + * + * where `taur = cos(2π/3)` and `taui = sin(2π/3)`. + * + * Because `W_3^0 = 1`, no twiddle multiplication is necessary beyond these constant factors (`taur` and `taui`). + */ + for ( k = 0; k < L; k++ ) { + ip1 = iptr( 0, k, 0, L, M, sc, oc ); // real part row 0 + ip2 = iptr( 0, k, 1, L, M, sc, oc ); // real part row 1 + ip3 = iptr( 0, k, 2, L, M, sc, oc ); // real part row 2 + + cr2 = cc[ ip2 ] + cc[ ip3 ]; + + io = optr( 0, 0, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + cr2; // row0 + row1 + row2 + + io = optr( 0, 2, k, M, so, oo ); + out[ io ] = TAUI * ( cc[ ip3 ] - cc[ ip2 ] ); // taui * ( row2-row1 ) + + io = optr( M-1, 1, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + ( TAUR * cr2 ); // row0 + ( taur * ( row1+row2 ) ) + } + // When the number of elements in a sub-sequence is less than `2`, there is nothing more to do, as the above butterfly produced the full result... + if ( M < 2 ) { + return; + } + /* + * Next, apply the general case where we need to loop through the non-trivial harmonics. + * + * For each harmonic `n = 1, ..., M/2-1`, we need to + * + * - recover three spectra from the three input rows. + * - apply radix-3 twiddle factors to rows 1 and 2, then combine with row 0 to form three output columns. + * - form the following + * + * x[3n] = X₀ + (W₁⋅X₁ + W₂⋅X₂) => column 0 + * x[3n+1] = X₀ + taur⋅(W₁⋅X₁ + W₂⋅X₂) + taui⋅i⋅(W₁⋅X₁ - W₂⋅X₂) => column 2 + * x[3n+2] = X₀ + taur⋅(W₁⋅X₁ + W₂⋅X₂) - taui⋅i⋅(W₁⋅X₁ - W₂⋅X₂) => column 1 + * + * The mirror index `im = M - i` selects the conjugate-symmetric partner, thus allowing the routine to read each symmetry pair only once. + */ + // Loop over each sub-sequence to be transformed... + for ( k = 0; k < L; k++ ) { + // Loop over the elements in each sub-sequence... + for ( i = 2; i < M; i += 2 ) { + im = M - i; // "mirror" index + + // Resolve twiddle factor indices for component 1: + it11 = ( (i-2)*st1 ) + ot1; // cos(θ) for component 1 + it12 = ( (i-1)*st1 ) + ot1; // sin(θ) for component 1 + + // Resolve twiddle factor indices for component 2: + it21 = ( (i-2)*st2 ) + ot2; // cos(θ) for component 2 + it22 = ( (i-1)*st2 ) + ot2; // sin(θ) for component 2 + + // Load component 1 data... + ip1 = iptr( i-1, k, 1, L, M, sc, oc ); // real part component 1 + ip2 = iptr( i, k, 1, L, M, sc, oc ); // imag part component 1 + + // Apply twiddle factor for component 1 + dr2 = ( twiddles1[ it11 ] * cc[ ip1 ] ) + ( twiddles1[ it12 ] * cc[ ip2 ] ); + di2 = ( twiddles1[ it11 ] * cc[ ip2 ] ) - ( twiddles1[ it12 ] * cc[ ip1 ] ); + + // Load component 2 data... + ip1 = iptr( i-1, k, 2, L, M, sc, oc ); // real part component 2 + ip2 = iptr( i, k, 2, L, M, sc, oc ); // imag part component 2 + + // Apply twiddle factor for component 2 + dr3 = ( twiddles2[ it21 ] * cc[ ip1 ] ) + ( twiddles2[ it22 ] * cc[ ip2 ] ); + di3 = ( twiddles2[ it21 ] * cc[ ip2 ] ) - ( twiddles2[ it22 ] * cc[ ip1 ] ); + + // Combine the twiddled components + cr2 = dr2 + dr3; // sum of real parts + ci2 = di2 + di3; // sum of imag parts + + // Load component 0 data... + ip1 = iptr( i-1, k, 0, L, M, sc, oc ); // real part component 0 + io = optr( i-1, 0, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + cr2; + + ip2 = iptr( i, k, 0, L, M, sc, oc ); // imag part component 0 + io = optr( i, 0, k, M, so, oo ); + out[ io ] = cc[ ip2 ] + ci2; + + // Intermediate results for components 1 and 2 + tr2 = cc[ ip1 ] + ( TAUR * cr2 ); + ti2 = cc[ ip2 ] + ( TAUR * ci2 ); + tr3 = TAUI * ( di2 - di3 ); + ti3 = TAUI * ( dr3 - dr2 ); + + // Output component 1 + io = optr( i-1, 2, k, M, so, oo ); + out[ io ] = tr2 + tr3; + + io = optr( im-1, 1, k, M, so, oo ); + out[ io ] = tr2 - tr3; + + // Output component 2 + io = optr( i, 2, k, M, so, oo ); + out[ io ] = ti2 + ti3; + + io = optr( im, 1, k, M, so, oo ); + out[ io ] = ti3 - ti2; + } + } +} + + +// EXPORTS // + +module.exports = radf3; From dd7b9bcd3d95f62bf94112e07d0bc1739521ef9f Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Mon, 11 May 2026 00:53:35 +0530 Subject: [PATCH 2/5] feat: add radf4 --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../fft/base/fftpack/rfftf/lib/radf2.js | 14 +- .../fft/base/fftpack/rfftf/lib/radf3.js | 11 +- .../fft/base/fftpack/rfftf/lib/radf4.js | 492 ++++++++++++++++++ 3 files changed, 505 insertions(+), 12 deletions(-) create mode 100644 lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js index 0d0af15bbf79..4b41a74e2ae3 100644 --- a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js @@ -304,7 +304,7 @@ function radf2( M, L, cc, sc, oc, out, so, oo, twiddles, st, ot ) { // eslint-di it2 = ( (i-1)*st ) + ot; // sin(θ) // Load the `even-odd` row... - ip1 = iptr( i-1, k, 1, L, M, sc, oc ); // c = Re(row1) + ip1 = iptr( i-1, k, 1, L, M, sc, oc ); // c = Re(row1) ip2 = iptr( i, k, 1, L, M, sc, oc ); // d = Im(row1) // eslint-disable-next-line stdlib/capitalized-comments, stdlib/empty-line-before-comment @@ -314,17 +314,17 @@ function radf2( M, L, cc, sc, oc, out, so, oo, twiddles, st, ot ) { // eslint-di // Load the `even+odd` row... ip1 = iptr( i, k, 0, L, M, sc, oc ); // b = Im(row0) - io = optr( i, 0, k, M, so, oo ); // Im(x[2n]) + io = optr( i, 0, k, M, so, oo ); // Im(x[2n]) out[ io ] = cc[ ip1 ] + ti2; - io = optr( im, 1, k, M, so, oo ); // Im(x[2n+1]) + io = optr( im, 1, k, M, so, oo ); // Im(x[2n+1]) out[ io ] = ti2 - cc[ ip1 ]; - ip1 = iptr( i-1, k, 0, L, M, sc, oc ); // a = Re(row0) - io = optr( i-1, 0, k, M, so, oo ); // Re(x[2n]) + ip1 = iptr( i-1, k, 0, L, M, sc, oc ); // a = Re(row0) + io = optr( i-1, 0, k, M, so, oo ); // Re(x[2n]) out[ io ] = cc[ ip1 ] + tr2; - io = optr( im-1, 1, k, M, so, oo ); // Re(x[2n+1]) + io = optr( im-1, 1, k, M, so, oo ); // Re(x[2n+1]) out[ io ] = cc[ ip1 ] - tr2; } } @@ -352,7 +352,7 @@ function radf2( M, L, cc, sc, oc, out, so, oo, twiddles, st, ot ) { // eslint-di for ( k = 0; k < L; k++ ) { ip1 = iptr( M-1, k, 1, L, M, sc, oc ); // -Re(O_n) io = optr( 0, 1, k, M, so, oo ); - out[ io ] = -cc[ ip1 ]; // -(-Re(O_n)) = Re(O_n) + out[ io ] = -cc[ ip1 ]; // -(-Re(O_n)) = Re(O_n) ip1 = iptr( M-1, k, 0, L, M, sc, oc ); // Re(E_n) io = optr( M-1, 0, k, M, so, oo ); diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js index 2fced4b2c683..3aae60305bf8 100644 --- a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js @@ -282,12 +282,13 @@ function radf3( M, L, cc, sc, oc, out, so, oo, twiddles1, st1, ot1, twiddles2, s * * For a radix-3 forward FFT of a real-valued signal `x` and `n = 0`, * + * taur = cos(2π/3) + * taui = sin(2π/3) + * * x[0] = row0 + row1 + row2 * x[2M] = taui * ( row2-row1 ) * x[M] = row0 + ( taur * ( row1+row2 ) ) * - * where `taur = cos(2π/3)` and `taui = sin(2π/3)`. - * * Because `W_3^0 = 1`, no twiddle multiplication is necessary beyond these constant factors (`taur` and `taui`). */ for ( k = 0; k < L; k++ ) { @@ -341,7 +342,7 @@ function radf3( M, L, cc, sc, oc, out, so, oo, twiddles1, st1, ot1, twiddles2, s // Load component 1 data... ip1 = iptr( i-1, k, 1, L, M, sc, oc ); // real part component 1 - ip2 = iptr( i, k, 1, L, M, sc, oc ); // imag part component 1 + ip2 = iptr( i, k, 1, L, M, sc, oc ); // imag part component 1 // Apply twiddle factor for component 1 dr2 = ( twiddles1[ it11 ] * cc[ ip1 ] ) + ( twiddles1[ it12 ] * cc[ ip2 ] ); @@ -349,7 +350,7 @@ function radf3( M, L, cc, sc, oc, out, so, oo, twiddles1, st1, ot1, twiddles2, s // Load component 2 data... ip1 = iptr( i-1, k, 2, L, M, sc, oc ); // real part component 2 - ip2 = iptr( i, k, 2, L, M, sc, oc ); // imag part component 2 + ip2 = iptr( i, k, 2, L, M, sc, oc ); // imag part component 2 // Apply twiddle factor for component 2 dr3 = ( twiddles2[ it21 ] * cc[ ip1 ] ) + ( twiddles2[ it22 ] * cc[ ip2 ] ); @@ -364,7 +365,7 @@ function radf3( M, L, cc, sc, oc, out, so, oo, twiddles1, st1, ot1, twiddles2, s io = optr( i-1, 0, k, M, so, oo ); out[ io ] = cc[ ip1 ] + cr2; - ip2 = iptr( i, k, 0, L, M, sc, oc ); // imag part component 0 + ip2 = iptr( i, k, 0, L, M, sc, oc ); // imag part component 0 io = optr( i, 0, k, M, so, oo ); out[ io ] = cc[ ip2 ] + ci2; diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js new file mode 100644 index 000000000000..c3e52effe63f --- /dev/null +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js @@ -0,0 +1,492 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +* +* +* ## Notice +* +* The original C code and copyright notice are from the [PFFFT library]{@link https://github.com/marton78/pffft/blob/0b4ee12c4ba45a4a8e567550c16d96d1679f50ce/src/fftpack.c}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (c) 2004 the University Corporation for Atmospheric +* Research ("UCAR"). All rights reserved. Developed by NCAR's +* Computational and Information Systems Laboratory, UCAR, +* www.cisl.ucar.edu. +* +* Redistribution and use of the Software in source and binary forms, +* with or without modification, is permitted provided that the +* following conditions are met: +* +* - Neither the names of NCAR's Computational and Information Systems +* Laboratory, the University Corporation for Atmospheric Research, +* nor the names of its sponsors or contributors may be used to +* endorse or promote products derived from this Software without +* specific prior written permission. +* +* - Redistributions of source code must retain the above copyright +* notices, this list of conditions, and the disclaimer below. +* +* - Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions, and the disclaimer below in the +* documentation and/or other materials provided with the +* distribution. +* +* THIS 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 CONTRIBUTORS OR COPYRIGHT +* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL, +* EXEMPLARY, OR CONSEQUENTIAL 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 WITH THE +* SOFTWARE. +* ``` +*/ + +/* eslint-disable max-len, max-statements */ + +'use strict'; + +// MODULES // + +var SQRT_HALF = require( '@stdlib/constants/float64/sqrt-half' ); + + +// FUNCTIONS // + +/** +* Resolves an index into the input array. +* +* ## Notes +* +* In a forward real FFT, the previous stage writes its results as four "rows" per sub-sequence. +* +* Thus, when reading from an input array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having four "rows" (components 0, 1, 2, and 3) and where each "row" is arranged as `M*L` contiguous elements corresponding to interleaved real and imaginary components. +* +* Accordingly, the following is a logical view of an input array (zero-based indexing) which contains `L = 3` transforms and in which each sub-sequence has length `M = 4`: +* +* ```text +* │ k = 0 k = 1 k = 2 +* │ ──────────────────────────────────────────────────────────────────────────→ k +* j = 0 │ cc(0,0,0) ... cc(3,0,0) cc(0,1,0) ... cc(3,1,0) cc(0,2,0) ... cc(3,2,0) +* │ +* j = 1 │ cc(0,0,1) ... cc(3,0,1) cc(0,1,1) ... cc(3,1,1) cc(0,2,1) ... cc(3,2,1) +* │ +* j = 2 │ cc(0,0,2) ... cc(3,0,2) cc(0,1,2) ... cc(3,1,2) cc(0,2,2) ... cc(3,2,2) +* │ +* j = 3 │ cc(0,0,3) ... cc(3,0,3) cc(0,1,3) ... cc(3,1,3) cc(0,2,3) ... cc(3,2,3) +* └───────────────────────────────────────────────────────────────────────────→ i +* ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 M-1 0 M-1 0 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short sub-sequence corresponding to one of the four component rows. +* - `j` selects between the four component rows (0, 1, 2, or 3). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | cc(0,0,0)...cc(3,0,0) ... cc(0,2,0)...cc(3,2,0) | cc(0,0,1)...cc(3,0,1) ... cc(0,2,1)...cc(3,2,1) | ... | cc(0,0,3)...cc(3,0,3) ... cc(0,2,3)...cc(3,2,3) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 LM LM-1 (L+1)M (L+1)M-1 (2L-1)M 2LM-1 3LM (3L+1)M-1 (4L-1)M 4LM-1 +* ``` +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} j - input row +* @param {NonNegativeInteger} L - number of sub-sequences +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the input array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = iptr( 0, 0, 0, L, M, stride, offset ); +* // returns 0 +* +* idx = iptr( 1, 0, 0, L, M, stride, offset ); +* // returns 1 +* +* idx = iptr( M-1, 0, 0, L, M, stride, offset ); +* // returns 3 +* +* idx = iptr( 0, 1, 0, L, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = iptr( M-1, L-1, 1, L, M, stride, offset ); +* // returns 23 +*/ +function iptr( i, k, j, L, M, stride, offset ) { + var n = i + ( ( k+(j*L) ) * M ); + return ( n*stride ) + offset; +} + +/** +* Resolves an index into the output array. +* +* ## Notes +* +* When writing to an output array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having four "columns" corresponding to the four components of a radix-4 stage (with real and imaginary parts of each component interleaved along each sub-sequence) and where each "column" has `M` elements. +* +* Accordingly, the following is a logical view of an output array (zero-based indexing) which contains `L = 3` transforms and in which each "column" sub-sequence has length `M = 4`: +* +* ```text +* j = 0 (component 0) j = 1 (component 1) j = 2 (component 2) j = 3 (component 3) +* k = 0 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┐ +* │ out(0,0,0) out(1,0,0) ... out(3,0,0) │ out(0,1,0) out(1,1,0) ... out(3,1,0) │ out(0,2,0) out(1,2,0) ... out(3,2,0) │ out(0,3,0) out(1,3,0) ... out(3,3,0) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 1 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,1) out(1,0,1) ... out(3,0,1) │ out(0,1,1) out(1,1,1) ... out(3,1,1) │ out(0,2,1) out(1,2,1) ... out(3,2,1) │ out(0,3,1) out(1,3,1) ... out(3,3,1) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 2 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,2) out(1,0,2) ... out(3,0,2) │ out(0,1,2) out(1,1,2) ... out(3,1,2) │ out(0,2,2) out(1,2,2) ... out(3,2,2) │ out(0,3,2) out(1,3,2) ... out(3,3,2) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┘ +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 1 M-1 0 1 M-1 0 1 M-1 0 1 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short "column" sub-sequence. +* - `j` selects which of the four components (0, 1, 2, or 3). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | out(0,0,0)...out(3,0,0) | out(0,1,0)...out(3,1,0) | out(0,2,0)...out(3,2,0) | out(0,3,0)...out(3,3,0) | out(0,0,1)...out(3,0,1) | ... | out(0,3,2)...out(3,3,2) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 M 2M-1 2M 3M-1 3M 4M-1 4M 5M-1 (4L-1)M 4LM-1 +* ``` +* +* As may be observed, when resolving an index in the output array, the `j` and `k` dimensions are swapped relative index resolution in the input array. This stems from `radf4` being only one stage in a multi-stage driver which alternates between using `cc` and `out` as workspace buffers. After each stage, the next stage reads what the previous stage wrote. +* +* Each stage expects a transpose, and, in order to avoid explicit transposition between the stages, we swap the last two logical dimensions while still maintaining cache locality within the inner loop logical dimension, as indexed by `i`. +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} j - index specifying which of the four complex components we are in (0, 1, 2, or 3) +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = optr( 0, 0, 0, M, stride, offset ); +* // returns 0 +* +* idx = optr( 1, 0, 0, M, stride, offset ); +* // returns 1 +* +* idx = optr( M-1, 0, 0, M, stride, offset ); +* // returns 3 +* +* idx = optr( 0, 1, 0, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = optr( M-1, 3, L-1, M, stride, offset ); +* // returns 47 +*/ +function optr( i, j, k, M, stride, offset ) { + var n = i + ( ( j+(k*4) ) * M ); + return ( n*stride ) + offset; +} + + +// MAIN // + +/** +* Performs one radix-4 stage within a forward Fourier transform for a real-valued sequence. +* +* @private +* @param {NonNegativeInteger} M - number of elements in each sub-sequence to be transformed +* @param {NonNegativeInteger} L - number of sub-sequences to be transformed +* @param {Float64Array} cc - input array containing the sub-sequences to be transformed +* @param {integer} sc - stride length for `cc` +* @param {NonNegativeInteger} oc - index offset for `cc` +* @param {Float64Array} out - output array containing transformed sub-sequences +* @param {integer} so - stride length for `out` +* @param {NonNegativeInteger} oo - index offset for `out` +* @param {Float64Array} twiddles1 - first array containing twiddle factors +* @param {integer} st1 - stride length for `twiddles1` +* @param {NonNegativeInteger} ot1 - index offset for `twiddles1` +* @param {Float64Array} twiddles2 - second array containing twiddle factors +* @param {integer} st2 - stride length for `twiddles2` +* @param {NonNegativeInteger} ot2 - index offset for `twiddles2` +* @param {Float64Array} twiddles3 - third array containing twiddle factors +* @param {integer} st3 - stride length for `twiddles3` +* @param {NonNegativeInteger} ot3 - index offset for `twiddles3` +* @returns {void} +*/ +function radf4( M, L, cc, sc, oc, out, so, oo, twiddles1, st1, ot1, twiddles2, st2, ot2, twiddles3, st3, ot3 ) { // eslint-disable-line max-params + var it11; + var it12; + var it21; + var it22; + var it31; + var it32; + var tr1; + var tr2; + var tr3; + var tr4; + var ti1; + var ti2; + var ti3; + var ti4; + var cr2; + var cr3; + var cr4; + var ci2; + var ci3; + var ci4; + var ip1; + var ip2; + var ip3; + var ip4; + var io; + var im; + var i; + var k; + + /* + * First, perform the core butterfly for each sub-sequence being transformed. + * + * In the following loop, we handle harmonic `n = 0` for every transform. As described for `iptr`, the input array is interpreted as a three-dimensional array, containing four "rows" per sub-sequence. + * + * row0 = input component 0 (j = 0) + * row1 = input component 1 (j = 1) + * row2 = input component 2 (j = 2) + * row3 = input component 3 (j = 3) + * + * For a radix-4 forward FFT of a real-valued signal `x` and `n = 0`, + * + * tr1 = row1 + row3 + * tr2 = row0 + row2 + * + * x[0] = tr1 + tr2 + * x[2M-1] = row0 - row2 + * x[2M] = row3 - row1 + * x[4M-1] = tr2 - tr1 + * + * Because `W_4^0 = 1`, no twiddle multiplication is necessary. + */ + for ( k = 0; k < L; k++ ) { + ip1 = iptr( 0, k, 0, L, M, sc, oc ); // real part row 0 + ip2 = iptr( 0, k, 1, L, M, sc, oc ); // real part row 1 + ip3 = iptr( 0, k, 2, L, M, sc, oc ); // real part row 2 + ip4 = iptr( 0, k, 3, L, M, sc, oc ); // real part row 3 + + tr1 = cc[ ip2 ] + cc[ ip4 ]; // row1 + row3 + tr2 = cc[ ip1 ] + cc[ ip3 ]; // row0 + row2 + + io = optr( 0, 0, k, M, so, oo ); + out[ io ] = tr1 + tr2; + + io = optr( M-1, 3, k, M, so, oo ); + out[ io ] = tr2 - tr1; + + io = optr( M-1, 1, k, M, so, oo ); + out[ io ] = cc[ ip1 ] - cc[ ip3 ]; + + io = optr( 0, 2, k, M, so, oo ); + out[ io ] = cc[ ip4 ] - cc[ ip2 ]; + } + // When the number of elements in a sub-sequence is equal to `1`, there is nothing more to do, as the above butterfly produced the full result... + if ( M < 2 ) { + return; + } + /* + * Next, apply the general case where we need to loop through the non-trivial harmonics. + * + * For each harmonic `n = 1, ..., M/2-1`, we need to + * + * - recover four spectra from the four input rows. + * - apply radix-4 twiddle factors to rows 1, 2, and 3, then combine with row 0 to form four output columns. + * - form the following + * + * x[4n] = X₀ + W₂·X₂ + (W₁·X₁ + W₃·X₃) => column 0 + * x[4n+1] = X₀ + W₂·X₂ - (W₁·X₁ + W₃·X₃) => column 3 + * x[4n+2] = X₀ - W₂·X₂ + i·(W₁·X₁ - W₃·X₃) => column 2 + * x[4n+3] = X₀ - W₂·X₂ - i·(W₁·X₁ - W₃·X₃) => column 1 + * + * The mirror index `im = M - i` selects the conjugate-symmetric partner, thus allowing the routine to read each symmetry pair only once. + */ + if ( M >= 3 ) { + // Loop over each sub-sequence to be transformed... + for ( k = 0; k < L; k++ ) { + // Loop over the elements in each sub-sequence... + for ( i = 2; i < M; i += 2 ) { + im = M - i; // "mirror" index + + // Resolve twiddle factor indices for component 1: + it11 = ( (i-2)*st1 ) + ot1; // cos(θ) for component 1 + it12 = ( (i-1)*st1 ) + ot1; // sin(θ) for component 1 + + // Resolve twiddle factor indices for component 2: + it21 = ( (i-2)*st2 ) + ot2; // cos(θ) for component 2 + it22 = ( (i-1)*st2 ) + ot2; // sin(θ) for component 2 + + // Resolve twiddle factor indices for component 3: + it31 = ( (i-2)*st3 ) + ot3; // cos(θ) for component 3 + it32 = ( (i-1)*st3 ) + ot3; // sin(θ) for component 3 + + // Load component 1 data and apply twiddle... + ip1 = iptr( i-1, k, 1, L, M, sc, oc ); // real part component 1 + ip2 = iptr( i, k, 1, L, M, sc, oc ); // imag part component 1 + + // Apply the twiddle factors for component 1: + cr2 = ( twiddles1[ it11 ] * cc[ ip1 ] ) + ( twiddles1[ it12 ] * cc[ ip2 ] ); + ci2 = ( twiddles1[ it11 ] * cc[ ip2 ] ) - ( twiddles1[ it12 ] * cc[ ip1 ] ); + + // Load component 2 data and apply twiddle... + ip1 = iptr( i-1, k, 2, L, M, sc, oc ); // real part component 2 + ip2 = iptr( i, k, 2, L, M, sc, oc ); // imag part component 2 + + // Apply the twiddle factors for component 2: + cr3 = ( twiddles2[ it21 ] * cc[ ip1 ] ) + ( twiddles2[ it22 ] * cc[ ip2 ] ); + ci3 = ( twiddles2[ it21 ] * cc[ ip2 ] ) - ( twiddles2[ it22 ] * cc[ ip1 ] ); + + // Load component 3 data and apply twiddle... + ip1 = iptr( i-1, k, 3, L, M, sc, oc ); // real part component 3 + ip2 = iptr( i, k, 3, L, M, sc, oc ); // imag part component 3 + + // Apply the twiddle factors for component 3: + cr4 = ( twiddles3[ it31 ] * cc[ ip1 ] ) + ( twiddles3[ it32 ] * cc[ ip2 ] ); + ci4 = ( twiddles3[ it31 ] * cc[ ip2 ] ) - ( twiddles3[ it32 ] * cc[ ip1 ] ); + + // Radix-4 butterfly intermediate computations: + tr1 = cr2 + cr4; + tr4 = cr4 - cr2; + ti1 = ci2 + ci4; + ti4 = ci2 - ci4; + + // Load component 0 data and compute intermediates... + ip1 = iptr( i-1, k, 0, L, M, sc, oc ); // real part component 0 + tr2 = cc[ ip1 ] + cr3; + tr3 = cc[ ip1 ] - cr3; + + ip2 = iptr( i, k, 0, L, M, sc, oc ); // imag part component 0 + ti2 = cc[ ip2 ] + ci3; + ti3 = cc[ ip2 ] - ci3; + + // Output column 0 (real part) + io = optr( i-1, 0, k, M, so, oo ); + out[ io ] = tr1 + tr2; + + // Output column 3 (real part) + io = optr( im-1, 3, k, M, so, oo ); + out[ io ] = tr2 - tr1; + + // Output column 0 (imag part) + io = optr( i, 0, k, M, so, oo ); + out[ io ] = ti1 + ti2; + + // Output column 3 (imag part) + io = optr( im, 3, k, M, so, oo ); + out[ io ] = ti1 - ti2; + + // Output column 2 (real part) + io = optr( i-1, 2, k, M, so, oo ); + out[ io ] = ti4 + tr3; + + // Output column 1 (real part) + io = optr( im-1, 1, k, M, so, oo ); + out[ io ] = tr3 - ti4; + + // Output column 2 (imag part) + io = optr( i, 2, k, M, so, oo ); + out[ io ] = tr4 + ti3; + + // Output column 1 (imag part) + io = optr( im, 1, k, M, so, oo ); + out[ io ] = tr4 - ti3; + } + } + + // When `M` is odd, there is no Nyquist pair to process, so we do not need to perform any further transformations... + if ( M%2 === 1 ) { + return; + } + } + + /* + * Lastly, handle the Nyquist frequency where `i = M-1` (i.e., the last element of each sub-sequence). + * + * When `M` is even, the Nyquist index is `i = M/2`. In this stage, we've stored that element at the end of each sub-sequence (i.e., `i = M-1` in the packed layout). + * + * At this point, only real components exist for rows 1 and 3 (imaginary parts are zero due to symmetry), and the twiddle multiplication simplifies to involve √(1/2) factors. + * + * In this case, + * + * ti1 = -√(1/2) ⋅ (row1 + row3) + * tr1 = √(1/2) ⋅ (row1 - row3) + * + * where + * + * row0 = Re(X₀) + * row1 = Re(X₁) (imaginary part is zero) + * row2 = Re(X₂) + * row3 = Re(X₃) (imaginary part is zero) + */ + for ( k = 0; k < L; k++ ) { + ip1 = iptr( M-1, k, 1, L, M, sc, oc ); // Re(X₁) + ip2 = iptr( M-1, k, 3, L, M, sc, oc ); // Re(X₃) + + ti1 = -SQRT_HALF * ( cc[ ip1 ] + cc[ ip2 ] ); + tr1 = SQRT_HALF * ( cc[ ip1 ] - cc[ ip2 ] ); + + ip3 = iptr( M-1, k, 0, L, M, sc, oc ); // Re(X₀) + ip4 = iptr( M-1, k, 2, L, M, sc, oc ); // Re(X₂) + + io = optr( M-1, 0, k, M, so, oo ); + out[ io ] = cc[ ip3 ] + tr1; + + io = optr( M-1, 2, k, M, so, oo ); + out[ io ] = cc[ ip3 ] - tr1; + + io = optr( 0, 1, k, M, so, oo ); + out[ io ] = ti1 - cc[ ip4 ]; + + io = optr( 0, 3, k, M, so, oo ); + out[ io ] = cc[ ip4 ] + ti1; + } +} + + +// EXPORTS // + +module.exports = radf4; From b1762929a192685d9db7ad41ffc328ecea7a1357 Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Tue, 12 May 2026 00:45:19 +0530 Subject: [PATCH 3/5] feat: add radf5 --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../fft/base/fftpack/rfftf/lib/radf2.js | 2 +- .../fft/base/fftpack/rfftf/lib/radf3.js | 2 +- .../fft/base/fftpack/rfftf/lib/radf4.js | 2 +- .../fft/base/fftpack/rfftf/lib/radf5.js | 517 ++++++++++++++++++ 4 files changed, 520 insertions(+), 3 deletions(-) create mode 100644 lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf5.js diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js index 4b41a74e2ae3..d7ac34b18d32 100644 --- a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf2.js @@ -105,7 +105,7 @@ * @param {NonNegativeInteger} L - number of sub-sequences * @param {NonNegativeInteger} M - sub-sequence length * @param {integer} stride - stride length of the input array -* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the input array * @returns {NonNegativeInteger} computed index * * @example diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js index 3aae60305bf8..54783ad72d3b 100644 --- a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js @@ -113,7 +113,7 @@ var TAUI = 0.866025403784439; // sin( 2π/3 ) * @param {NonNegativeInteger} L - number of sub-sequences * @param {NonNegativeInteger} M - sub-sequence length * @param {integer} stride - stride length of the input array -* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the input array * @returns {NonNegativeInteger} computed index * * @example diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js index c3e52effe63f..d657a5cf8005 100644 --- a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf4.js @@ -114,7 +114,7 @@ var SQRT_HALF = require( '@stdlib/constants/float64/sqrt-half' ); * @param {NonNegativeInteger} L - number of sub-sequences * @param {NonNegativeInteger} M - sub-sequence length * @param {integer} stride - stride length of the input array -* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the input array * @returns {NonNegativeInteger} computed index * * @example diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf5.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf5.js new file mode 100644 index 000000000000..6ad5dea1f911 --- /dev/null +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf5.js @@ -0,0 +1,517 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +* +* +* ## Notice +* +* The original C code and copyright notice are from the [PFFFT library]{@link https://github.com/marton78/pffft/blob/0b4ee12c4ba45a4a8e567550c16d96d1679f50ce/src/fftpack.c}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (c) 2004 the University Corporation for Atmospheric +* Research ("UCAR"). All rights reserved. Developed by NCAR's +* Computational and Information Systems Laboratory, UCAR, +* www.cisl.ucar.edu. +* +* Redistribution and use of the Software in source and binary forms, +* with or without modification, is permitted provided that the +* following conditions are met: +* +* - Neither the names of NCAR's Computational and Information Systems +* Laboratory, the University Corporation for Atmospheric Research, +* nor the names of its sponsors or contributors may be used to +* endorse or promote products derived from this Software without +* specific prior written permission. +* +* - Redistributions of source code must retain the above copyright +* notices, this list of conditions, and the disclaimer below. +* +* - Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions, and the disclaimer below in the +* documentation and/or other materials provided with the +* distribution. +* +* THIS 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 CONTRIBUTORS OR COPYRIGHT +* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL, +* EXEMPLARY, OR CONSEQUENTIAL 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 WITH THE +* SOFTWARE. +* ``` +*/ + +/* eslint-disable max-len, max-statements */ + +'use strict'; + +// MODULES // + +var cospi = require( '@stdlib/math/base/special/cospi' ); +var sinpi = require( '@stdlib/math/base/special/sinpi' ); + + +// VARIABLES // + +// cos( 2π/5 ) = 0.30901699437494734 +var TR11 = cospi( 2.0 / 5.0 ); + +// sin( 2π/5 ) = 0.9510565162951536 +var TI11 = sinpi( 2.0 / 5.0 ); + +// cos( 4π/5 ) = -0.8090169943749475 +var TR12 = cospi( 4.0 / 5.0 ); + +// sin( 4π/5 ) = 0.587785252292473 +var TI12 = sinpi( 4.0 / 5.0 ); + + +// FUNCTIONS // + +/** +* Resolves an index into the input array. +* +* ## Notes +* +* In a forward real FFT, the previous stage writes its results as five "rows" per sub-sequence. +* +* Thus, when reading from an input array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having five "rows" (components 0, 1, 2, 3, and 4) and where each "row" is arranged as `M*L` contiguous elements corresponding to interleaved real and imaginary components. +* +* Accordingly, the following is a logical view of an input array (zero-based indexing) which contains `L = 3` transforms and in which each sub-sequence has length `M = 4`: +* +* ```text +* │ k = 0 k = 1 k = 2 +* │ ──────────────────────────────────────────────────────────────────────────→ k +* j = 0 │ cc(0,0,0) ... cc(3,0,0) cc(0,1,0) ... cc(3,1,0) cc(0,2,0) ... cc(3,2,0) +* │ +* j = 1 │ cc(0,0,1) ... cc(3,0,1) cc(0,1,1) ... cc(3,1,1) cc(0,2,1) ... cc(3,2,1) +* │ +* j = 2 │ cc(0,0,2) ... cc(3,0,2) cc(0,1,2) ... cc(3,1,2) cc(0,2,2) ... cc(3,2,2) +* │ +* j = 3 │ cc(0,0,3) ... cc(3,0,3) cc(0,1,3) ... cc(3,1,3) cc(0,2,3) ... cc(3,2,3) +* │ +* j = 4 │ cc(0,0,4) ... cc(3,0,4) cc(0,1,4) ... cc(3,1,4) cc(0,2,4) ... cc(3,2,4) +* └───────────────────────────────────────────────────────────────────────────→ i +* ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 M-1 0 M-1 0 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short sub-sequence corresponding to one of the five component rows. +* - `j` selects between the five component rows (0, 1, 2, 3, or 4). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | cc(0,0,0)...cc(3,0,0) ... cc(0,2,0)...cc(3,2,0) | cc(0,0,1)...cc(3,0,1) ... cc(0,2,1)...cc(3,2,1) | ... | cc(0,0,4)...cc(3,0,4) ... cc(0,2,4)...cc(3,2,4) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 LM LM-1 (L+1)M (L+1)M-1 (2L-1)M 2LM-1 4LM (4L+1)M-1 (5L-1)M 5LM-1 +* ``` +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} j - input row +* @param {NonNegativeInteger} L - number of sub-sequences +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the input array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the input array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = iptr( 0, 0, 0, L, M, stride, offset ); +* // returns 0 +* +* idx = iptr( 1, 0, 0, L, M, stride, offset ); +* // returns 1 +* +* idx = iptr( M-1, 0, 0, L, M, stride, offset ); +* // returns 3 +* +* idx = iptr( 0, 1, 0, L, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = iptr( M-1, L-1, 1, L, M, stride, offset ); +* // returns 23 +*/ +function iptr( i, k, j, L, M, stride, offset ) { + var n = i + ( ( k+(j*L) ) * M ); + return ( n*stride ) + offset; +} + +/** +* Resolves an index into the output array. +* +* ## Notes +* +* When writing to an output array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having five "columns" corresponding to the five components of a radix-5 stage (with real and imaginary parts of each component interleaved along each sub-sequence) and where each "column" has `M` elements. +* +* Accordingly, the following is a logical view of an output array (zero-based indexing) which contains `L = 3` transforms and in which each "column" sub-sequence has length `M = 4`: +* +* ```text +* j = 0 (component 0) j = 1 (component 1) j = 2 (component 2) j = 3 (component 3) j = 4 (component 4) +* k = 0 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┐ +* │ out(0,0,0) out(1,0,0) ... out(3,0,0) │ out(0,1,0) out(1,1,0) ... out(3,1,0) │ out(0,2,0) out(1,2,0) ... out(3,2,0) │ out(0,3,0) out(1,3,0) ... out(3,3,0) │ out(0,4,0) out(1,4,0) ... out(3,4,0) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 1 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,1) out(1,0,1) ... out(3,0,1) │ out(0,1,1) out(1,1,1) ... out(3,1,1) │ out(0,2,1) out(1,2,1) ... out(3,2,1) │ out(0,3,1) out(1,3,1) ... out(3,3,1) │ out(0,4,1) out(1,4,1) ... out(3,4,1) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┤ +* k = 2 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┤ +* │ out(0,0,2) out(1,0,2) ... out(3,0,2) │ out(0,1,2) out(1,1,2) ... out(3,1,2) │ out(0,2,2) out(1,2,2) ... out(3,2,2) │ out(0,3,2) out(1,3,2) ... out(3,3,2) │ out(0,4,2) out(1,4,2) ... out(3,4,2) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┘ +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 1 M-1 0 1 M-1 0 1 M-1 0 1 M-1 0 1 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short "column" sub-sequence. +* - `j` selects which of the five components (0, 1, 2, 3, or 4). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | out(0,0,0)...out(3,0,0) | out(0,1,0)...out(3,1,0) | out(0,2,0)...out(3,2,0) | out(0,3,0)...out(3,3,0) | out(0,4,0)...out(3,4,0) | out(0,0,1)...out(3,0,1) | ... | out(0,4,2)...out(3,4,2) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 M 2M-1 2M 3M-1 3M 4M-1 4M 5M-1 5M 6M-1 (5L-1)M 5LM-1 +* ``` +* +* As may be observed, when resolving an index in the output array, the `j` and `k` dimensions are swapped relative index resolution in the input array. This stems from `radf5` being only one stage in a multi-stage driver which alternates between using `cc` and `out` as workspace buffers. After each stage, the next stage reads what the previous stage wrote. +* +* Each stage expects a transpose, and, in order to avoid explicit transposition between the stages, we swap the last two logical dimensions while still maintaining cache locality within the inner loop logical dimension, as indexed by `i`. +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} j - index specifying which of the five complex components we are in (0, 1, 2, 3, or 4) +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = optr( 0, 0, 0, M, stride, offset ); +* // returns 0 +* +* idx = optr( 1, 0, 0, M, stride, offset ); +* // returns 1 +* +* idx = optr( M-1, 0, 0, M, stride, offset ); +* // returns 3 +* +* idx = optr( 0, 1, 0, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = optr( M-1, 4, L-1, M, stride, offset ); +* // returns 59 +*/ +function optr( i, j, k, M, stride, offset ) { + var n = i + ( ( j+(k*5) ) * M ); + return ( n*stride ) + offset; +} + + +// MAIN // + +/** +* Performs one radix-5 stage within a forward Fourier transform for a real-valued sequence. +* +* @private +* @param {NonNegativeInteger} M - number of elements in each sub-sequence to be transformed +* @param {NonNegativeInteger} L - number of sub-sequences to be transformed +* @param {Float64Array} cc - input array containing the sub-sequences to be transformed +* @param {integer} sc - stride length for `cc` +* @param {NonNegativeInteger} oc - index offset for `cc` +* @param {Float64Array} out - output array containing transformed sub-sequences +* @param {integer} so - stride length for `out` +* @param {NonNegativeInteger} oo - index offset for `out` +* @param {Float64Array} twiddles1 - first array containing twiddle factors +* @param {integer} st1 - stride length for `twiddles1` +* @param {NonNegativeInteger} ot1 - index offset for `twiddles1` +* @param {Float64Array} twiddles2 - second array containing twiddle factors +* @param {integer} st2 - stride length for `twiddles2` +* @param {NonNegativeInteger} ot2 - index offset for `twiddles2` +* @param {Float64Array} twiddles3 - third array containing twiddle factors +* @param {integer} st3 - stride length for `twiddles3` +* @param {NonNegativeInteger} ot3 - index offset for `twiddles3` +* @param {Float64Array} twiddles4 - fourth array containing twiddle factors +* @param {integer} st4 - stride length for `twiddles4` +* @param {NonNegativeInteger} ot4 - index offset for `twiddles4` +* @returns {void} +*/ +function radf5( M, L, cc, sc, oc, out, so, oo, twiddles1, st1, ot1, twiddles2, st2, ot2, twiddles3, st3, ot3, twiddles4, st4, ot4 ) { // eslint-disable-line max-params + var it11; + var it12; + var it21; + var it22; + var it31; + var it32; + var it41; + var it42; + var ci2; + var ci3; + var ci4; + var ci5; + var cr2; + var cr3; + var cr4; + var cr5; + var dr2; + var dr3; + var dr4; + var dr5; + var di2; + var di3; + var di4; + var di5; + var ti2; + var ti3; + var ti4; + var ti5; + var tr2; + var tr3; + var tr4; + var tr5; + var ip1; + var ip2; + var ip3; + var ip4; + var ip5; + var io; + var im; + var i; + var k; + + /* + * First, perform the core butterfly for each sub-sequence being transformed. + * + * In the following loop, we handle harmonic `n = 0` for every transform. As described for `iptr`, the input array is interpreted as a three-dimensional array, containing five "rows" per sub-sequence. + * + * row0 = input component 0 (j = 0) + * row1 = input component 1 (j = 1) + * row2 = input component 2 (j = 2) + * row3 = input component 3 (j = 3) + * row4 = input component 4 (j = 4) + * + * For a radix-5 forward FFT of a real-valued signal `x` and `n = 0`, + * + * cr2 = row4 + row1 + * ci5 = row4 - row1 + * cr3 = row3 + row2 + * ci4 = row3 - row2 + * + * x[0] = row0 + cr2 + cr3 + * x[M] = row0 + ( TR11*cr2 ) + ( TR12*cr3 ) + * x[2M] = ( TI11*ci5 ) + ( TI12*ci4 ) + * x[3M] = row0 + ( TR12*cr2 ) + ( TR11*cr3 ) + * x[4M] = ( TI12*ci5 ) - ( TI11*ci4 ) + * + * Because `W_5^0 = 1`, no frequency-dependent twiddle multiplication is necessary beyond these constant factors. + */ + for ( k = 0; k < L; k++ ) { + ip1 = iptr( 0, k, 0, L, M, sc, oc ); // row 0 + ip2 = iptr( 0, k, 1, L, M, sc, oc ); // row 1 + ip3 = iptr( 0, k, 2, L, M, sc, oc ); // row 2 + ip4 = iptr( 0, k, 3, L, M, sc, oc ); // row 3 + ip5 = iptr( 0, k, 4, L, M, sc, oc ); // row 4 + + cr2 = cc[ ip5 ] + cc[ ip2 ]; // row4 + row1 + ci5 = cc[ ip5 ] - cc[ ip2 ]; // row4 - row1 + cr3 = cc[ ip4 ] + cc[ ip3 ]; // row3 + row2 + ci4 = cc[ ip4 ] - cc[ ip3 ]; // row3 - row2 + + io = optr( 0, 0, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + cr2 + cr3; + + io = optr( M-1, 1, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + ( TR11 * cr2 ) + ( TR12 * cr3 ); + + io = optr( 0, 2, k, M, so, oo ); + out[ io ] = ( TI11 * ci5 ) + ( TI12 * ci4 ); + + io = optr( M-1, 3, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + ( TR12 * cr2 ) + ( TR11 * cr3 ); + + io = optr( 0, 4, k, M, so, oo ); + out[ io ] = ( TI12 * ci5 ) - ( TI11 * ci4 ); + } + // When the number of elements in a sub-sequence is less than `2`, there is nothing more to do, as the above butterfly produced the full result... + if ( M < 2 ) { + return; + } + /* + * Next, apply the general case where we need to loop through the non-trivial harmonics. + * + * For each harmonic `n = 1, ..., M/2-1`, we need to + * + * - recover five spectra from the five input rows. + * - apply radix-5 twiddle factors to rows 1, 2, 3, and 4, then combine with row 0 to form five output columns. + * - form the following + * + * x[5n] = X₀ + (W₁⋅X₁ + W₄⋅X₄) + (W₂⋅X₂ + W₃⋅X₃) => column 0 + * x[5n+1] = X₀ + TR11⋅(W₁⋅X₁ + W₄⋅X₄) + TR12⋅(W₂⋅X₂ + W₃⋅X₃) - i⋅(TI11⋅(W₁⋅X₁ - W₄⋅X₄) + TI12⋅(W₂⋅X₂ - W₃⋅X₃)) => column 1 + * x[5n+2] = X₀ + TR11⋅(W₁⋅X₁ + W₄⋅X₄) + TR12⋅(W₂⋅X₂ + W₃⋅X₃) + i⋅(TI11⋅(W₁⋅X₁ - W₄⋅X₄) + TI12⋅(W₂⋅X₂ - W₃⋅X₃)) => column 2 + * x[5n+3] = X₀ + TR12⋅(W₁⋅X₁ + W₄⋅X₄) + TR11⋅(W₂⋅X₂ + W₃⋅X₃) - i⋅(TI12⋅(W₁⋅X₁ - W₄⋅X₄) - TI11⋅(W₂⋅X₂ - W₃⋅X₃)) => column 3 + * x[5n+4] = X₀ + TR12⋅(W₁⋅X₁ + W₄⋅X₄) + TR11⋅(W₂⋅X₂ + W₃⋅X₃) + i⋅(TI12⋅(W₁⋅X₁ - W₄⋅X₄) - TI11⋅(W₂⋅X₂ - W₃⋅X₃)) => column 4 + * + * The mirror index `im = M - i` selects the conjugate-symmetric partner, thus allowing the routine to read each symmetry pair only once. + */ + // Loop over each sub-sequence to be transformed... + for ( k = 0; k < L; k++ ) { + // Loop over the elements in each sub-sequence... + for ( i = 2; i < M; i += 2 ) { + im = M - i; // "mirror" index + + // Resolve twiddle factor indices for component 1: + it11 = ( (i-2)*st1 ) + ot1; // cos(θ) for component 1 + it12 = ( (i-1)*st1 ) + ot1; // sin(θ) for component 1 + + // Resolve twiddle factor indices for component 2: + it21 = ( (i-2)*st2 ) + ot2; // cos(θ) for component 2 + it22 = ( (i-1)*st2 ) + ot2; // sin(θ) for component 2 + + // Resolve twiddle factor indices for component 3: + it31 = ( (i-2)*st3 ) + ot3; // cos(θ) for component 3 + it32 = ( (i-1)*st3 ) + ot3; // sin(θ) for component 3 + + // Resolve twiddle factor indices for component 4: + it41 = ( (i-2)*st4 ) + ot4; // cos(θ) for component 4 + it42 = ( (i-1)*st4 ) + ot4; // sin(θ) for component 4 + + // Load component 1 data and apply twiddle... + ip1 = iptr( i-1, k, 1, L, M, sc, oc ); // real part component 1 + ip2 = iptr( i, k, 1, L, M, sc, oc ); // imag part component 1 + + // Apply the twiddle factors for component 1: + dr2 = ( twiddles1[ it11 ] * cc[ ip1 ] ) + ( twiddles1[ it12 ] * cc[ ip2 ] ); // Re(dr2) + di2 = ( twiddles1[ it11 ] * cc[ ip2 ] ) - ( twiddles1[ it12 ] * cc[ ip1 ] ); // Im(di2) + + // Load component 2 data and apply twiddle... + ip1 = iptr( i-1, k, 2, L, M, sc, oc ); // real part component 2 + ip2 = iptr( i, k, 2, L, M, sc, oc ); // imag part component 2 + + // Apply the twiddle factors for component 2: + dr3 = ( twiddles2[ it21 ] * cc[ ip1 ] ) + ( twiddles2[ it22 ] * cc[ ip2 ] ); // Re(dr3) + di3 = ( twiddles2[ it21 ] * cc[ ip2 ] ) - ( twiddles2[ it22 ] * cc[ ip1 ] ); // Im(di3) + + // Load component 3 data and apply twiddle... + ip1 = iptr( i-1, k, 3, L, M, sc, oc ); // real part component 3 + ip2 = iptr( i, k, 3, L, M, sc, oc ); // imag part component 3 + + // Apply the twiddle factors for component 3: + dr4 = ( twiddles3[ it31 ] * cc[ ip1 ] ) + ( twiddles3[ it32 ] * cc[ ip2 ] ); // Re(dr4) + di4 = ( twiddles3[ it31 ] * cc[ ip2 ] ) - ( twiddles3[ it32 ] * cc[ ip1 ] ); // Im(di4) + + // Load component 4 data and apply twiddle... + ip1 = iptr( i-1, k, 4, L, M, sc, oc ); // real part component 4 + ip2 = iptr( i, k, 4, L, M, sc, oc ); // imag part component 4 + + // Apply the twiddle factors for component 4: + dr5 = ( twiddles4[ it41 ] * cc[ ip1 ] ) + ( twiddles4[ it42 ] * cc[ ip2 ] ); // Re(dr5) + di5 = ( twiddles4[ it41 ] * cc[ ip2 ] ) - ( twiddles4[ it42 ] * cc[ ip1 ] ); // Im(di5) + + // Radix-5 butterfly intermediate computations (symmetric/antisymmetric combinations): + cr2 = dr2 + dr5; + ci5 = dr5 - dr2; + cr5 = di2 - di5; + ci2 = di5 + di2; + cr3 = dr3 + dr4; + ci4 = dr4 - dr3; + cr4 = di3 - di4; + ci3 = di3 + di4; + + // Output column 0 (real part) + ip1 = iptr( i-1, k, 0, L, M, sc, oc ); + io = optr( i-1, 0, k, M, so, oo ); + out[ io ] = cc[ ip1 ] + cr2 + cr3; + + // Output column 0 (imag part) + ip2 = iptr( i, k, 0, L, M, sc, oc ); + io = optr( i, 0, k, M, so, oo ); + out[ io ] = cc[ ip2 ] + ci2 + ci3; + + // Load component 0 data and compute intermediates... + tr2 = cc[ ip1 ] + ( TR11 * cr2 ) + ( TR12 * cr3 ); + ti2 = cc[ ip2 ] + ( TR11 * ci2 ) + ( TR12 * ci3 ); + tr3 = cc[ ip1 ] + ( TR12 * cr2 ) + ( TR11 * cr3 ); + ti3 = cc[ ip2 ] + ( TR12 * ci2 ) + ( TR11 * ci3 ); + + // Additional intermediates: + tr5 = ( TI11 * cr5 ) + ( TI12 * cr4 ); + ti5 = ( TI11 * ci5 ) + ( TI12 * ci4 ); + tr4 = ( TI12 * cr5 ) - ( TI11 * cr4 ); + ti4 = ( TI12 * ci5 ) - ( TI11 * ci4 ); + + // Output column 2 (real part) + io = optr( i-1, 2, k, M, so, oo ); + out[ io ] = tr2 + tr5; + + // Output column 1 (real part) + io = optr( im-1, 1, k, M, so, oo ); + out[ io ] = tr2 - tr5; + + // Output column 2 (imag part) + io = optr( i, 2, k, M, so, oo ); + out[ io ] = ti2 + ti5; + + // Output column 1 (imag part) + io = optr( im, 1, k, M, so, oo ); + out[ io ] = ti5 - ti2; + + // Output column 4 (real part) + io = optr( i-1, 4, k, M, so, oo ); + out[ io ] = tr3 + tr4; + + // Output column 3 (real part) + io = optr( im-1, 3, k, M, so, oo ); + out[ io ] = tr3 - tr4; + + // Output column 4 (imag part) + io = optr( i, 4, k, M, so, oo ); + out[ io ] = ti3 + ti4; + + // Output column 3 (imag part) + io = optr( im, 3, k, M, so, oo ); + out[ io ] = ti4 - ti3; + } + } +} + + +// EXPORTS // + +module.exports = radf5; From 5f3fc63cf0db06f23c35f665b887db06f0caac93 Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Tue, 12 May 2026 00:49:35 +0530 Subject: [PATCH 4/5] refactor: use sinpi and cospi --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../@stdlib/fft/base/fftpack/rfftf/lib/radf3.js | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js index 54783ad72d3b..914ca19853df 100644 --- a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radf3.js @@ -60,10 +60,19 @@ 'use strict'; +// MODULES // + +var cospi = require( '@stdlib/math/base/special/cospi' ); +var sinpi = require( '@stdlib/math/base/special/sinpi' ); + + // VARIABLES // -var TAUR = -0.5; // cos( 2π/3 ) -var TAUI = 0.866025403784439; // sin( 2π/3 ) +// cos( 2π/3 ) = -0.5 +var TAUR = cospi( 2.0 / 3.0 ); + +// sin( 2π/3 ) = 0.866025403784439 +var TAUI = sinpi( 2.0 / 3.0 ); // FUNCTIONS // From dd6338d1fdde9ebea792d0adccae7cb85b9750ed Mon Sep 17 00:00:00 2001 From: GUNJ JOSHI Date: Wed, 20 May 2026 01:55:20 +0530 Subject: [PATCH 5/5] feat: add radfg --- type: pre_commit_static_analysis_report description: Results of running static analysis checks when committing changes. report: - task: lint_filenames status: passed - task: lint_editorconfig status: passed - task: lint_markdown status: na - task: lint_package_json status: na - task: lint_repl_help status: na - task: lint_javascript_src status: passed - task: lint_javascript_cli status: na - task: lint_javascript_examples status: na - task: lint_javascript_tests status: na - task: lint_javascript_benchmarks status: na - task: lint_python status: na - task: lint_r status: na - task: lint_c_src status: na - task: lint_c_examples status: na - task: lint_c_benchmarks status: na - task: lint_c_tests_fixtures status: na - task: lint_shell status: na - task: lint_typescript_declarations status: passed - task: lint_typescript_tests status: na - task: lint_license_headers status: passed --- --- .../fft/base/fftpack/rfftf/lib/radfg.js | 761 ++++++++++++++++++ 1 file changed, 761 insertions(+) create mode 100644 lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radfg.js diff --git a/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radfg.js b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radfg.js new file mode 100644 index 000000000000..c6f7dfc925a8 --- /dev/null +++ b/lib/node_modules/@stdlib/fft/base/fftpack/rfftf/lib/radfg.js @@ -0,0 +1,761 @@ +/** +* @license Apache-2.0 +* +* Copyright (c) 2026 The Stdlib Authors. +* +* Licensed under the Apache License, Version 2.0 (the "License"); +* you may not use this file except in compliance with the License. +* You may obtain a copy of the License at +* +* http://www.apache.org/licenses/LICENSE-2.0 +* +* Unless required by applicable law or agreed to in writing, software +* distributed under the License is distributed on an "AS IS" BASIS, +* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +* See the License for the specific language governing permissions and +* limitations under the License. +* +* +* ## Notice +* +* The original C code and copyright notice are from the [PFFFT library]{@link https://github.com/marton78/pffft/blob/0b4ee12c4ba45a4a8e567550c16d96d1679f50ce/src/fftpack.c}. The implementation follows the original, but has been modified for JavaScript. +* +* ```text +* Copyright (c) 2004 the University Corporation for Atmospheric +* Research ("UCAR"). All rights reserved. Developed by NCAR's +* Computational and Information Systems Laboratory, UCAR, +* www.cisl.ucar.edu. +* +* Redistribution and use of the Software in source and binary forms, +* with or without modification, is permitted provided that the +* following conditions are met: +* +* - Neither the names of NCAR's Computational and Information Systems +* Laboratory, the University Corporation for Atmospheric Research, +* nor the names of its sponsors or contributors may be used to +* endorse or promote products derived from this Software without +* specific prior written permission. +* +* - Redistributions of source code must retain the above copyright +* notices, this list of conditions, and the disclaimer below. +* +* - Redistributions in binary form must reproduce the above copyright +* notice, this list of conditions, and the disclaimer below in the +* documentation and/or other materials provided with the +* distribution. +* +* THIS 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 CONTRIBUTORS OR COPYRIGHT +* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL, +* EXEMPLARY, OR CONSEQUENTIAL 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 WITH THE +* SOFTWARE. +* ``` +*/ + +/* eslint-disable max-len, max-statements, max-lines-per-function */ + +'use strict'; + +// MODULES // + +var cospi = require( '@stdlib/math/base/special/cospi' ); +var sinpi = require( '@stdlib/math/base/special/sinpi' ); +var floor = require( '@stdlib/math/base/special/floor' ); + + +// FUNCTIONS // + +/** +* Resolves an index into the input array. +* +* ## Notes +* +* In a forward real FFT, the previous stage writes its results as `R` "rows" per sub-sequence. +* +* Thus, when reading from an input array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having `R` "rows" (components 0, 1, 2, ..., R-1) and where each "row" is arranged as `M*L` contiguous elements corresponding to interleaved real and imaginary components. +* +* Accordingly, the following is a logical view of an input array (zero-based indexing) which contains `L = 3` transforms and in which each sub-sequence has length `M = 4`: +* +* ```text +* │ k = 0 k = 1 k = 2 +* │ ─────────────────────────────────────────────────────────────────────────────────────→ k +* j = 0 │ cc(0,0,0) ... cc(3,0,0) cc(0,1,0) ... cc(3,1,0) cc(0,2,0) ... cc(3,2,0) +* │ +* j = 1 │ cc(0,0,1) ... cc(3,0,1) cc(0,1,1) ... cc(3,1,1) cc(0,2,1) ... cc(3,2,1) +* │ +* j = 2 │ cc(0,0,2) ... cc(3,0,2) cc(0,1,2) ... cc(3,1,2) cc(0,2,2) ... cc(3,2,2) +* │ +* j = 3 │ cc(0,0,3) ... cc(3,0,3) cc(0,1,3) ... cc(3,1,3) cc(0,2,3) ... cc(3,2,3) +* │ +* ... │ ... ... ... +* │ +* j=R-1 │ cc(0,0,R-1) ... cc(3,0,R-1) cc(0,1,R-1) ... cc(3,1,R-1) cc(0,2,R-1) ... cc(3,2,R-1) +* └──────────────────────────────────────────────────────────────────────────────────────→ i +* ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 M-1 0 M-1 0 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short sub-sequence corresponding to one of the R component rows. +* - `j` selects between the R component rows (0, 1, 2, ..., R-1). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | cc(0,0,0)...cc(3,0,0) ... cc(0,2,0)...cc(3,2,0) | cc(0,0,1)...cc(3,0,1) ... cc(0,2,1)...cc(3,2,1) | ... | cc(0,0,R-1)...cc(3,0,R-1) ... cc(0,2,R-1)...cc(3,2,R-1) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 LM LM-1 (L+1)M (L+1)M-1 (2L-1)M 2LM-1 (R-1)LM ((R-1)L+1)M-1 (RL-1)M RLM-1 +* ``` +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} j - input row +* @param {NonNegativeInteger} L - number of sub-sequences +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the input array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the input array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = iptr( 0, 0, 0, L, M, stride, offset ); +* // returns 0 +* +* idx = iptr( 1, 0, 0, L, M, stride, offset ); +* // returns 1 +* +* idx = iptr( M-1, 0, 0, L, M, stride, offset ); +* // returns 3 +* +* idx = iptr( 0, 1, 0, L, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = iptr( M-1, L-1, 6, L, M, stride, offset ); +* // returns 83 +*/ +function iptr( i, k, j, L, M, stride, offset ) { + var n = i + ( ( k+(j*L) ) * M ); + return ( n*stride ) + offset; +} + +/** +* Resolves an index into the output array. +* +* ## Notes +* +* When writing to an output array stored in linear memory, we can reinterpret the array as a three-dimensional logical view containing `L` independent sub-sequences having `R` "columns" corresponding to the `R` components of a radix-R stage (with real and imaginary parts of each component interleaved along each sub-sequence) and where each "column" has `M` elements. +* +* Accordingly, the following is a logical view of an output array (zero-based indexing) which contains `L = 3` transforms and in which each "column" sub-sequence has length `M = 4` for an arbitrary radix `R`: +* +* ```text +* j = 0 (component 0) j = 1 (component 1) j = 2 (component 2) ... j = R-1 (component R-1) +* k = 0 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────────────┐ +* │ out(0,0,0) out(1,0,0) ... out(3,0,0) │ out(0,1,0) out(1,1,0) ... out(3,1,0) │ out(0,2,0) out(1,2,0) ... out(3,2,0) │ ... │ out(0,R-1,0) out(1,R-1,0) ... out(3,R-1,0) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────────────┤ +* k = 1 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────────────┤ +* │ out(0,0,1) out(1,0,1) ... out(3,0,1) │ out(0,1,1) out(1,1,1) ... out(3,1,1) │ out(0,2,1) out(1,2,1) ... out(3,2,1) │ ... │ out(0,R-1,1) out(1,R-1,1) ... out(3,R-1,1) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────────────┤ +* k = 2 ─┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────┬───────────────────────────────────────────────┤ +* │ out(0,0,2) out(1,0,2) ... out(3,0,2) │ out(0,1,2) out(1,1,2) ... out(3,1,2) │ out(0,2,2) out(1,2,2) ... out(3,2,2) │ ... │ out(0,R-1,2) out(1,R-1,2) ... out(3,R-1,2) │ +* └───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────┴───────────────────────────────────────────────┘ +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* i = 0 1 M-1 0 1 M-1 0 1 M-1 0 1 M-1 +* ``` +* +* In the above, +* +* - `i` is the fastest varying index, which walks within one short "column" sub-sequence. +* - `j` selects which of the R components we are in (0, 1, 2, ..., R-1). +* - `k` specifies the index of one of the `L` independent transforms we are processing. +* +* In linear memory, the three-dimensional logical view is arranged as follows: +* +* ```text +* | out(0,0,0)...out(3,0,0) | out(0,1,0)...out(3,1,0) | out(0,2,0)...out(3,2,0) | ... | out(0,R-1,0)...out(3,R-1,0) | out(0,0,1)...out(3,0,1) | ... | out(0,R-1,2)...out(3,R-1,2) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 M-1 M 2M-1 2M 3M-1 (R-1)M RM-1 RM (R+1)M-1 (RL-1)M RLM-1 +* ``` +* +* As may be observed, when resolving an index in the output array, the `j` and `k` dimensions are swapped relative index resolution in the input array. This stems from `radfg` being only one stage in a multi-stage driver which alternates between using `cc` and `out` as workspace buffers. After each stage, the next stage reads what the previous stage wrote. +* +* Each stage expects a transpose, and, in order to avoid explicit transposition between the stages, we swap the last two logical dimensions while still maintaining cache locality within the inner loop logical dimension, as indexed by `i`. +* +* @private +* @param {NonNegativeInteger} i - index of an element within a sub-sequence +* @param {NonNegativeInteger} j - index specifying which of the R complex components we are in (0, 1, ..., R-1) +* @param {NonNegativeInteger} k - index of the sub-sequence being transformed +* @param {NonNegativeInteger} R - radix +* @param {NonNegativeInteger} M - sub-sequence length +* @param {integer} stride - stride length of the output array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* +* var idx = optr( 0, 0, 0, 7, M, stride, offset ); +* // returns 0 +* +* idx = optr( 1, 0, 0, 7, M, stride, offset ); +* // returns 1 +* +* idx = optr( M-1, 0, 0, 7, M, stride, offset ); +* // returns 3 +* +* idx = optr( 0, 1, 0, 7, M, stride, offset ); +* // returns 4 +* +* // ... +* +* idx = optr( M-1, 6, L-1, 7, M, stride, offset ); +* // returns 83 +*/ +function optr( i, j, k, R, M, stride, offset ) { + var n = i + ( ( j+(k*R) ) * M ); + return ( n*stride ) + offset; +} + +/** +* Resolves an index into a flattened workspace array. +* +* ## Notes +* +* During the general radix stage, flattened workspace arrays store `ML = M * L` elements for each of the `R` component columns. Both the flattened input and output workspace arrays share the same `[ML, R]` logical layout. +* +* Thus, when accessing a flattened workspace array stored in linear memory, we can reinterpret the array as a two-dimensional logical view having `R` component columns, where each column has `ML` elements. +* +* Accordingly, the following is a logical view of a flattened workspace array (zero-based indexing) having flattened dimension `ML` and arbitrary radix `R`: +* +* ```text +* │ j = 0 j = 1 j = 2 ... j = R-1 +* │ ────────────────────────────────────────────────────────────────────────→ j +* ik = 0 │ ws(0,0) ws(0,1) ws(0,2) ... ws(0,R-1) +* │ +* ik = 1 │ ws(1,0) ws(1,1) ws(1,2) ... ws(1,R-1) +* │ +* ik = 2 │ ws(2,0) ws(2,1) ws(2,2) ... ws(2,R-1) +* │ +* ik = 3 │ ws(3,0) ws(3,1) ws(3,2) ... ws(3,R-1) +* │ +* ... │ ... ... ... ... ... +* │ +* ik = ML-1 │ ws(ML-1,0) ws(ML-1,1) ws(ML-1,2) ... ws(ML-1,R-1) +* ``` +* +* In the above, +* +* - `ik` is the fastest varying index, which walks along the flattened sub-sequence dimension. +* - `j` selects between the `R` workspace component columns (0, 1, 2, ..., R-1). +* +* In linear memory, the two-dimensional logical view is arranged as follows: +* +* ```text +* | ws(0,0)...ws(ML-1,0) | ws(0,1)...ws(ML-1,1) | ws(0,2)...ws(ML-1,2) | ... | ws(0,R-1)...ws(ML-1,R-1) | +* ↑ ↑ ↑ ↑ ↑ ↑ ↑ ↑ +* 0 ML-1 ML 2ML-1 2ML 3ML-1 (R-1)ML RML-1 +* ``` +* +* Here, the original `M` and `L` dimensions have been collapsed into the single flattened `ik` dimension. Each workspace component `j` therefore occupies one contiguous block of `ML` elements in linear memory. +* +* @private +* @param {NonNegativeInteger} ik - index of an element within a flattened sub-sequence +* @param {NonNegativeInteger} j - index specifying which of the R workspace component columns we are in (0, 1, ..., R-1) +* @param {NonNegativeInteger} ML - flattened sub-sequence length +* @param {integer} stride - stride length of the flattened workspace array +* @param {NonNegativeInteger} offset - index specifying the first indexed element in the flattened workspace array +* @returns {NonNegativeInteger} computed index +* +* @example +* var stride = 1; +* var offset = 0; +* +* var M = 4; // sub-sequence length +* var L = 3; // number of sub-sequences +* var ML = M * L; // flattened dimension ( 4*3 = 12 ) +* +* var idx = fptr( 0, 0, ML, stride, offset ); +* // returns 0 +* +* idx = fptr( 1, 0, ML, stride, offset ); +* // returns 1 +* +* idx = fptr( ML-1, 0, ML, stride, offset ); +* // returns 11 +* +* idx = fptr( 0, 1, ML, stride, offset ); +* // returns 12 +* +* // ... +* +* idx = fptr( ML-1, 6, ML, stride, offset ); +* // returns 83 +*/ +function fptr( ik, j, ML, stride, offset ) { + var n = ik + ( j*ML ); + return ( n*stride ) + offset; +} + + +// MAIN // + +/** +* Performs one general radix stage within a forward Fourier transform for a real-valued sequence. +* +* @private +* @param {NonNegativeInteger} M - number of elements in each sub-sequence to be transformed +* @param {NonNegativeInteger} R - radix of the transform +* @param {NonNegativeInteger} L - number of sub-sequences to be transformed +* @param {NonNegativeInteger} ML - number of elements in each flattened sub-sequence (`M*L`) +* @param {Float64Array} cc - input array containing the sub-sequences to be transformed +* @param {integer} sc - stride length for `cc` +* @param {NonNegativeInteger} oc - index offset for `cc` +* @param {Float64Array} c1 - input workspace array containing intermediate sub-sequences to be transformed +* @param {integer} sc1 - stride length for `c1` +* @param {NonNegativeInteger} oc1 - index offset for `c1` +* @param {Float64Array} c2 - flattened input workspace array containing intermediate sub-sequences +* @param {integer} sc2 - stride length for `c2` +* @param {NonNegativeInteger} oc2 - index offset for `c2` +* @param {Float64Array} ch - output array containing transformed sequences +* @param {integer} sch - stride length for `ch` +* @param {NonNegativeInteger} och - index offset for `ch` +* @param {Float64Array} ch2 - flattened output workspace array containing intermediate sub-sequences +* @param {integer} sch2 - stride length for `ch2` +* @param {NonNegativeInteger} och2 - index offset for `ch2` +* @param {Float64Array} twiddles - array of twiddle factors +* @param {integer} stw - stride length for `twiddles` +* @param {NonNegativeInteger} otw - index offset for `twiddles` +* @returns {void} +*/ +function radfg( M, R, L, ML, cc, sc, oc, c1, sc1, oc1, c2, sc2, oc2, ch, sch, och, ch2, sch2, och2, twiddles, stw, otw ) { // eslint-disable-line max-params + var ar1h; + var ar2h; + var idij; + var ic2o; + var ic1o; + var ich; + var icc; + var Rph; + var Rp1; + var ih1; + var ih2; + var ih3; + var ih4; + var ic1; + var ic2; + var ic3; + var ic4; + var it1; + var it2; + var io1; + var io2; + var io3; + var io4; + var dc2; + var ai1; + var ai2; + var ar1; + var ar2; + var ds2; + var nbd; + var dcp; + var dsp; + var im; + var ik; + var is; + var jc; + var lc; + var j2; + var i; + var j; + var k; + var l; + + // Compute the basic rotation angle for the radix-R DFT matrix: + dcp = cospi( 2.0 / R ); // cos(2π/R) + dsp = sinpi( 2.0 / R ); // sin(2π/R) + + // Compute half the radix, rounded up, which is used as the exclusive upper bound for conjugate-symmetric component loops: + Rph = floor( ( R + 1 ) / 2 ); + Rp1 = R + 1; + + // Compute the number of non-DC complex harmonics in each sub-sequence: + nbd = floor( ( M - 1 ) / 2 ); + + /* + * First, initialize the work arrays for the general-radix butterfly. + * + * At this stage, each sub-sequence has already been split into `R` radix components. The `j = 0` component is carried in the flattened work arrays, while the remaining components `j = 1, ..., R-1` are carried in the three-dimensional work arrays. + * + * For harmonic `n = 0`, we only copy the DC terms. + * + * For harmonics `n = 1, ..., floor((M-1)/2)`, each nonzero radix component is multiplied by its stage twiddle factor and stored in `ch`. + */ + if ( M === 1 ) { + // Copy the DC column from the flattened output workspace back to the flattened input workspace... + for ( ik = 0; ik < ML; ik++ ) { + ic2o = fptr( ik, 0, ML, sc2, oc2 ); + ich = fptr( ik, 0, ML, sch2, och2 ); + c2[ ic2o ] = ch2[ ich ]; + } + } else { + // Copy the DC column into the flattened output workspace... + for ( ik = 0; ik < ML; ik++ ) { + ich = fptr( ik, 0, ML, sch2, och2 ); + ic2o = fptr( ik, 0, ML, sc2, oc2 ); + ch2[ ich ] = c2[ ic2o ]; + } + + // Copy the DC terms for the remaining radix components... + for ( j = 1; j < R; j++ ) { + for ( k = 0; k < L; k++ ) { + ih1 = optr( 0, k, j, L, M, sch, och ); + ic1 = iptr( 0, k, j, L, M, sc1, oc1 ); + ch[ ih1 ] = c1[ ic1 ]; + } + } + + // Apply twiddle factors to the non-DC harmonics of the remaining radix components... + if ( nbd <= L ) { + // Loop over harmonics before sub-sequences... + is = 0; + for ( j = 1; j < R; j++ ) { + idij = is; + for ( i = 2; i < M; i += 2 ) { + for ( k = 0; k < L; k++ ) { + // Resolve output indices in ch: + ih1 = optr( i-1, k, j, L, M, sch, och ); // real part + ih2 = optr( i, k, j, L, M, sch, och ); // imaginary part + + // Resolve indices in the input workspace: + ic1 = iptr( i-1, k, j, L, M, sc1, oc1 ); // real part + ic2 = iptr( i, k, j, L, M, sc1, oc1 ); // imaginary part + + // Resolve twiddle factor indices: + it1 = ( idij * stw ) + otw; // cos(θ) + it2 = ( (idij+1) * stw ) + otw; // sin(θ) + + // Apply the twiddle factor: + ch[ ih1 ] = ( twiddles[ it1 ] * c1[ ic1 ] ) + ( twiddles[ it2 ] * c1[ ic2 ] ); // Re(ch) + ch[ ih2 ] = ( twiddles[ it1 ] * c1[ ic2 ] ) - ( twiddles[ it2 ] * c1[ ic1 ] ); // Im(ch) + } + idij += 2; + } + is += M; + } + } else { + // Loop over sub-sequences before harmonics... + is = 0; + for ( j = 1; j < R; j++ ) { + for ( k = 0; k < L; k++ ) { + idij = is; + for ( i = 2; i < M; i += 2 ) { + // Resolve output indices in ch: + ih1 = optr( i-1, k, j, L, M, sch, och ); // real part + ih2 = optr( i, k, j, L, M, sch, och ); // imaginary part + + // Resolve indices in the input workspace: + ic1 = iptr( i-1, k, j, L, M, sc1, oc1 ); // real part + ic2 = iptr( i, k, j, L, M, sc1, oc1 ); // imaginary part + + // Resolve twiddle factor indices: + it1 = ( idij * stw ) + otw; // cos(θ) + it2 = ( (idij+1) * stw ) + otw; // sin(θ) + + // Apply the twiddle factor: + ch[ ih1 ] = ( twiddles[ it1 ] * c1[ ic1 ] ) + ( twiddles[ it2 ] * c1[ ic2 ] ); // Re(ch) + ch[ ih2 ] = ( twiddles[ it1 ] * c1[ ic2 ] ) - ( twiddles[ it2 ] * c1[ ic1 ] ); // Im(ch) + + idij += 2; + } + } + is += M; + } + } + + /* + * Next, combine mirrored radix components, storing their sums and differences in c1. + */ + if ( nbd >= L ) { + // Loop over sub-sequences before harmonics... + for ( j = 1; j < Rph; j++ ) { + jc = Rp1 - j - 1; // "mirror" index + for ( k = 0; k < L; k++ ) { + for ( i = 2; i < M; i += 2 ) { + // Resolve ch indices for component j: + ih1 = optr( i-1, k, j, L, M, sch, och ); // Re(ch[j]) + ih2 = optr( i, k, j, L, M, sch, och ); // Im(ch[j]) + + // Resolve ch indices for conjugate component jc: + ih3 = optr( i-1, k, jc, L, M, sch, och ); // Re(ch[jc]) + ih4 = optr( i, k, jc, L, M, sch, och ); // Im(ch[jc]) + + // Resolve input-workspace indices for component j: + ic1 = iptr( i-1, k, j, L, M, sc1, oc1 ); // Re(component j) + ic2 = iptr( i, k, j, L, M, sc1, oc1 ); // Im(component j) + + // Resolve input-workspace indices for component jc: + ic3 = iptr( i-1, k, jc, L, M, sc1, oc1 ); // Re(component jc) + ic4 = iptr( i, k, jc, L, M, sc1, oc1 ); // Im(component jc) + + // Form the sum and difference combinations: + c1[ ic1 ] = ch[ ih1 ] + ch[ ih3 ]; // Re(j) + Re(jc) + c1[ ic3 ] = ch[ ih2 ] - ch[ ih4 ]; // Im(j) - Im(jc) + c1[ ic2 ] = ch[ ih2 ] + ch[ ih4 ]; // Im(j) + Im(jc) + c1[ ic4 ] = ch[ ih3 ] - ch[ ih1 ]; // Re(jc) - Re(j) + } + } + } + } else { + // Loop over harmonics before sub-sequences... + for ( j = 1; j < Rph; j++ ) { + jc = Rp1 - j - 1; // "mirror" index + for ( i = 2; i < M; i += 2 ) { + for ( k = 0; k < L; k++ ) { + // Resolve ch indices for component j: + ih1 = optr( i-1, k, j, L, M, sch, och ); // Re(ch[j]) + ih2 = optr( i, k, j, L, M, sch, och ); // Im(ch[j]) + + // Resolve ch indices for conjugate component jc: + ih3 = optr( i-1, k, jc, L, M, sch, och ); // Re(ch[jc]) + ih4 = optr( i, k, jc, L, M, sch, och ); // Im(ch[jc]) + + // Resolve input-workspace indices for component j: + ic1 = iptr( i-1, k, j, L, M, sc1, oc1 ); // Re(component j) + ic2 = iptr( i, k, j, L, M, sc1, oc1 ); // Im(component j) + + // Resolve input-workspace indices for component jc: + ic3 = iptr( i-1, k, jc, L, M, sc1, oc1 ); // Re(component jc) + ic4 = iptr( i, k, jc, L, M, sc1, oc1 ); // Im(component jc) + + // Form the sum and difference combinations: + c1[ ic1 ] = ch[ ih1 ] + ch[ ih3 ]; // Re(j) + Re(jc) + c1[ ic3 ] = ch[ ih2 ] - ch[ ih4 ]; // Im(j) - Im(jc) + c1[ ic2 ] = ch[ ih2 ] + ch[ ih4 ]; // Im(j) + Im(jc) + c1[ ic4 ] = ch[ ih3 ] - ch[ ih1 ]; // Re(jc) - Re(j) + } + } + } + } + } + + // Combine the DC terms for each conjugate-symmetric component pair... + for ( j = 1; j < Rph; j++ ) { + jc = Rp1 - j - 1; // "mirror" index + for ( k = 0; k < L; k++ ) { + ih1 = optr( 0, k, j, L, M, sch, och ); + ih2 = optr( 0, k, jc, L, M, sch, och ); + + ic1 = iptr( 0, k, j, L, M, sc1, oc1 ); + ic2 = iptr( 0, k, jc, L, M, sc1, oc1 ); + + c1[ ic1 ] = ch[ ih1 ] + ch[ ih2 ]; + c1[ ic2 ] = ch[ ih2 ] - ch[ ih1 ]; + } + } + + /* + * Next, apply the radix-`R` DFT matrix. + * + * For each mirrored output pair, accumulate the contributions of the radix components in the flattened work array. + * + * The required rotation factors are generated recursively from the previous ones: + * + * W_l = W_{l-1} ⋅ W_1 + * W_{lj} = W_{l(j-1)} ⋅ W_l + * + * Here, `W_1` is the basic radix-`R` rotation, `W_{l-1}` and `W_l` are the previous and current outer factors, and `W_{l(j-1)}` and `W_{lj}` are the previous and current inner factors. + */ + ar1 = 1.0; // Re(W_0) + ai1 = 0.0; // Im(W_0) + for ( l = 1; l < Rph; l++ ) { + lc = Rp1 - l - 1; // "mirror" index + + // Advance the rotation factor by one step: W_l = W_{l-1} ⋅ W_1 + ar1h = ( dcp * ar1 ) - ( dsp * ai1 ); // Re(W_l) + ai1 = ( dcp * ai1 ) + ( dsp * ar1 ); // Im(W_l) + ar1 = ar1h; + + // Accumulate the contributions from components 1 and R-1: + for ( ik = 0; ik < ML; ik++ ) { + ich = fptr( ik, l, ML, sch2, och2 ); + ic1o = fptr( ik, lc, ML, sch2, och2 ); + + ic2o = fptr( ik, 0, ML, sc2, oc2 ); + ic1 = fptr( ik, 1, ML, sc2, oc2 ); + ic2 = fptr( ik, R-1, ML, sc2, oc2 ); + + ch2[ ich ] = c2[ ic2o ] + ( ar1 * c2[ ic1 ] ); + ch2[ ic1o ] = ai1 * c2[ ic2 ]; + } + + // Save W_l for the inner recurrence: + dc2 = ar1; + ds2 = ai1; + + // Initialize W_{l·1}: + ar2 = ar1; + ai2 = ai1; + + // Accumulate the remaining component pairs: + for ( j = 2; j < Rph; j++ ) { + jc = Rp1 - j - 1; // "mirror" index + + // Advance the nested rotation factor by one step: W_{l*j} = W_{l*(j-1)} ⋅ W_l + ar2h = ( dc2 * ar2 ) - ( ds2 * ai2 ); // Re(W_{l*j}) + ai2 = ( dc2 * ai2 ) + ( ds2 * ar2 ); // Im(W_{l*j}) + ar2 = ar2h; + + // Accumulate the contributions from components j and jc: + for ( ik = 0; ik < ML; ik++ ) { + ich = fptr( ik, l, ML, sch2, och2 ); + ic1o = fptr( ik, lc, ML, sch2, och2 ); + + ic1 = fptr( ik, j, ML, sc2, oc2 ); + ic2 = fptr( ik, jc, ML, sc2, oc2 ); + + ch2[ ich ] += ar2 * c2[ ic1 ]; + ch2[ ic1o ] += ai2 * c2[ ic2 ]; + } + } + } + + // Accumulate the nonzero radix components into the DC output... + for ( j = 1; j < Rph; j++ ) { + for ( ik = 0; ik < ML; ik++ ) { + ich = fptr( ik, 0, ML, sch2, och2 ); + ic1 = fptr( ik, j, ML, sc2, oc2 ); + ch2[ ich ] += c2[ ic1 ]; + } + } + + // Copy the first output column from ch to cc... + if ( M >= L ) { + for ( k = 0; k < L; k++ ) { + for ( i = 0; i < M; i++ ) { + icc = iptr( i, k, 0, L, M, sc, oc ); + ih1 = optr( i, k, 0, L, M, sch, och ); + cc[ icc ] = ch[ ih1 ]; + } + } + } else { + for ( i = 0; i < M; i++ ) { + for ( k = 0; k < L; k++ ) { + icc = iptr( i, k, 0, L, M, sc, oc ); + ih1 = optr( i, k, 0, L, M, sch, och ); + cc[ icc ] = ch[ ih1 ]; + } + } + } + + // Store DC harmonics for conjugate pairs in the output array... + for ( j = 1; j < Rph; j++ ) { + jc = Rp1 - j - 1; // "mirror" index + j2 = 2 * j; // output column index + for ( k = 0; k < L; k++ ) { + io1 = iptr( M-1, k, j2-1, L, M, sc, oc ); + io2 = iptr( 0, k, j2, L, M, sc, oc ); + + ih1 = optr( 0, k, j, L, M, sch, och ); + ih2 = optr( 0, k, jc, L, M, sch, och ); + + cc[ io1 ] = ch[ ih1 ]; + cc[ io2 ] = ch[ ih2 ]; + } + } + + // When M = 1, there are no non-DC harmonics to process, so we're done... + if ( M === 1 ) { + return; + } + + /* + * Finally, store the non-DC harmonics in folded format. + * + * For each mirror pair, write the non-DC harmonics from `ch` to the two output columns as the required sum and difference combinations. + */ + if ( nbd >= L ) { + // Loop over sub-sequences before harmonics... + for ( j = 1; j < Rph; j++ ) { + jc = Rp1 - j - 1; // "mirror" index + j2 = 2 * j; // output column index + for ( k = 0; k < L; k++ ) { + for ( i = 2; i < M; i += 2 ) { + im = M - i; // "mirror" harmonic index + + // Resolve ch indices for component j: + ih1 = optr( i-1, k, j, L, M, sch, och ); // Re(ch[j]) + ih2 = optr( i, k, j, L, M, sch, och ); // Im(ch[j]) + + // Resolve ch indices for conjugate component jc: + ih3 = optr( i-1, k, jc, L, M, sch, och ); // Re(ch[jc]) + ih4 = optr( i, k, jc, L, M, sch, och ); // Im(ch[jc]) + + // Resolve cc output indices: + io1 = iptr( i-1, k, j2, L, M, sc, oc ); + io2 = iptr( im-1, k, j2-1, L, M, sc, oc ); + io3 = iptr( i, k, j2, L, M, sc, oc ); + io4 = iptr( im, k, j2-1, L, M, sc, oc ); + + // Store the sum and difference combinations: + cc[ io1 ] = ch[ ih1 ] + ch[ ih3 ]; + cc[ io2 ] = ch[ ih1 ] - ch[ ih3 ]; + cc[ io3 ] = ch[ ih2 ] + ch[ ih4 ]; + cc[ io4 ] = ch[ ih4 ] - ch[ ih2 ]; + } + } + } + } else { + // Loop over harmonics before sub-sequences... + for ( j = 1; j < Rph; j++ ) { + jc = Rp1 - j - 1; // "mirror" index + j2 = 2 * j; // output column index + for ( i = 2; i < M; i += 2 ) { + im = M - i; // "mirror" harmonic index + for ( k = 0; k < L; k++ ) { + // Resolve ch indices for component j: + ih1 = optr( i-1, k, j, L, M, sch, och ); // Re(ch[j]) + ih2 = optr( i, k, j, L, M, sch, och ); // Im(ch[j]) + + // Resolve ch indices for conjugate component jc: + ih3 = optr( i-1, k, jc, L, M, sch, och ); // Re(ch[jc]) + ih4 = optr( i, k, jc, L, M, sch, och ); // Im(ch[jc]) + + // Resolve cc output indices: + io1 = iptr( i-1, k, j2, L, M, sc, oc ); + io2 = iptr( im-1, k, j2-1, L, M, sc, oc ); + io3 = iptr( i, k, j2, L, M, sc, oc ); + io4 = iptr( im, k, j2-1, L, M, sc, oc ); + + // Store the sum and difference combinations: + cc[ io1 ] = ch[ ih1 ] + ch[ ih3 ]; + cc[ io2 ] = ch[ ih1 ] - ch[ ih3 ]; + cc[ io3 ] = ch[ ih2 ] + ch[ ih4 ]; + cc[ io4 ] = ch[ ih4 ] - ch[ ih2 ]; + } + } + } + } +} + + +// EXPORTS // + +module.exports = radfg;