Press n or j to go to the next uncovered block, b, p or k for the previous block.
| 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 | 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 4x 4x 4x 4x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 4x 4x 4x 4x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 1x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 2x 1x 1x 1x 1x 1x | /**
* @license Apache-2.0
*
* Copyright (c) 2025 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.
* ```
*/
/* 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 of the previous stage.
* - `j` selects between the `even + odd` and `even - odd` row of the previous stage.
* - `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 MP1;
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_N^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 + 1 - i` selects the conjugate-symmetric partner, thus allowing the routine to read each symmetry pair only once.
*/
if ( M >= 3 ) {
MP1 = M + 1;
// Loop over each sub-subsequence 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 = MP1 - i; // "mirror" index
// Resolve twiddle factor indices:
it1 = ( (i-1)*st ) + ot; // cos(θ)
it2 = ( i*st ) + ot; // sin(θ)
// Load the `even-odd` row...
ip1 = iptr( i, k, 1, L, M, sc, oc ); // c = Re(row1)
ip2 = iptr( i+1, k, 1, L, M, sc, oc ); // d = Im(row1)
// 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+1, k, 0, L, M, sc, oc ); // b = Im(row0)
io = optr( i+1, 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, k, 0, L, M, sc, oc ); // a = Re(row0)
io = optr( i, 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;
|