Skip to content

Commit 72a102e

Browse files
kgrytekgryte
committed
chore: fix offset bugs in implementation and refactor and fix normalization in tests
--- 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: passed - 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: na - task: lint_typescript_tests status: na - task: lint_license_headers status: passed ---
1 parent 5082d3c commit 72a102e

File tree

14 files changed

+450
-200
lines changed

14 files changed

+450
-200
lines changed

‎lib/node_modules/@stdlib/fft/base/fftpack/lib/cfftb.js‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@ var cfftb1 = require( './cfftb1.js' );
7474
*
7575
* - Invoking `cfftf` prior to invoking this function will multiply an input sequence by `N`.
7676
*
77-
* - Prior to invoking this function, `wsave` must be initialized by calling `cffti(N, wsave)`.
77+
* - Prior to invoking this function, `workspace` must be initialized by calling `cffti(N, workspace)`.
7878
*
7979
* - This function is more efficient when `N` is the product of small prime numbers.
8080
*
@@ -99,14 +99,14 @@ var cfftb1 = require( './cfftb1.js' );
9999
* @param {NonNegativeInteger} N - length of the sequence to transform
100100
* @param {Float64Array} c - real-valued input array containing interleaved real and imaginary components corresponding to the sequence to be transformed
101101
* @param {NonNegativeInteger} offsetC - starting index for `c`
102-
* @param {Float64Array} wsave - real-valued workspace array containing pre-computed values
103-
* @param {NonNegativeInteger} offsetW - starting index for `wsave`
102+
* @param {Float64Array} workspace - real-valued workspace array containing pre-computed values
103+
* @param {NonNegativeInteger} offsetW - starting index for `workspace`
104104
* @returns {void}
105105
*
106106
* @example
107107
* // TODO
108108
*/
109-
function cfftb( N, c, offsetC, wsave, offsetW ) {
109+
function cfftb( N, c, offsetC, workspace, offsetW ) {
110110
var offsetT;
111111
var offsetF;
112112

@@ -146,10 +146,10 @@ function cfftb( N, c, offsetC, wsave, offsetW ) {
146146
*
147147
* In short, twiddle factors are cached roots of unity that allow each stage of the algorithm to rotate data quickly and predictably.
148148
*/
149-
offsetT = N * 2; // index offset for twiddle factors
150-
offsetF = offsetT * 2; // index offset for factors describing the sub-transforms
149+
offsetT = offsetW + ( 2*N ); // index offset for twiddle factors
150+
offsetF = offsetT + ( 2*N ); // index offset for factors describing the sub-transforms
151151

152-
cfftb1( N, c, offsetC, wsave, offsetW, wsave, offsetW+offsetT, wsave, offsetW+offsetF ); // eslint-disable-line max-len
152+
cfftb1( N, c, offsetC, workspace, offsetW, workspace, offsetT, workspace, offsetF ); // eslint-disable-line max-len
153153
}
154154

155155

‎lib/node_modules/@stdlib/fft/base/fftpack/lib/cfftb1.js‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -83,7 +83,7 @@ var passfb = require( './passfb.js' );
8383
* @param {NonNegativeInteger} chOffset - starting index for `ch`
8484
* @param {Float64Array} wa - workspace array for storing twiddle factors
8585
* @param {NonNegativeInteger} waOffset - starting index for `wa`
86-
* @param {Int32Array} ifac - workspace array for storing factorization information
86+
* @param {Int32Array} ifac - workspace array for storing factorization results
8787
* @param {NonNegativeInteger} ifacOffset - starting index for `ifac`
8888
* @returns {void}
8989
*/
@@ -156,7 +156,7 @@ function cfftb1( N, c, cOffset, ch, chOffset, wa, waOffset, ifac, ifacOffset ) {
156156
}
157157
// Now that we've finished computing the transforms, copy over the final results to the input array...
158158
for ( i = 0; i < N*2; i++ ) {
159-
c[ i+cOffset ] = ch[ i+chOffset ];
159+
c[ i+cOffset ] = ch[ i+chOffset ]; // FIXME: use `blas/base/dcopy`
160160
}
161161
}
162162

‎lib/node_modules/@stdlib/fft/base/fftpack/lib/ch_ref.js‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -67,12 +67,12 @@
6767
* @param {NonNegativeInteger} a1 - index of first dimension
6868
* @param {NonNegativeInteger} a2 - index of second dimension
6969
* @param {NonNegativeInteger} a3 - index of third dimension
70-
* @param {NonNegativeInteger} l1 - length parameter related to the FFT stage
70+
* @param {NonNegativeInteger} radix - length parameter related to the FFT stage
7171
* @param {NonNegativeInteger} ido - dimension order
7272
* @returns {NonNegativeInteger} - calculated index
7373
*/
74-
function chRef( a1, a2, a3, l1, ido ) {
75-
return ( ( (a3*l1)+a2 ) * ido ) + a1;
74+
function chRef( a1, a2, a3, radix, ido ) {
75+
return ( ( (a3*radix)+a2 ) * ido ) + a1;
7676
}
7777

7878

‎lib/node_modules/@stdlib/fft/base/fftpack/lib/decompose.js‎

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -97,11 +97,28 @@ var floor = require( '@stdlib/math/base/special/floor' );
9797
* var factors = [ 0, 0, 0, 0, 0, 0, 0 ];
9898
*
9999
* // Factorize the sequence length into its prime factors:
100-
* var numFactors = decompose( 630, initial, factors, 1, 0 );
100+
* var numFactors = decompose( N, initial, factors, 1, 0 );
101101
* // returns 5
102102
*
103103
* var f = factors.slice();
104104
* // returns [ 630, 5, 2, 3, 3, 5, 7 ]
105+
*
106+
* @example
107+
* // Specify an initial list of potential divisors:
108+
* var initial = [ 3, 4, 2, 5 ]; // as found in FFTPACK
109+
*
110+
* // Define a sequence length:
111+
* var N = 8;
112+
*
113+
* // Initialize an array for storing factorization results:
114+
* var factors = [ 0, 0, 0, 0 ];
115+
*
116+
* // Factorize the sequence length into its prime factors:
117+
* var numFactors = decompose( N, initial, factors, 1, 0 );
118+
* // returns 2
119+
*
120+
* var f = factors.slice();
121+
* // returns [ 8, 2, 2, 4 ]
105122
*/
106123
function decompose( N, initial, out, stride, offset ) {
107124
var divisor;

‎lib/node_modules/@stdlib/fft/base/fftpack/lib/print_real_workspace.js‎

Lines changed: 64 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,64 @@
1+
/**
2+
* @license Apache-2.0
3+
*
4+
* Copyright (c) 2025 The Stdlib Authors.
5+
*
6+
* Licensed under the Apache License, Version 2.0 (the "License");
7+
* you may not use this file except in compliance with the License.
8+
* You may obtain a copy of the License at
9+
*
10+
* http://www.apache.org/licenses/LICENSE-2.0
11+
*
12+
* Unless required by applicable law or agreed to in writing, software
13+
* distributed under the License is distributed on an "AS IS" BASIS,
14+
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15+
* See the License for the specific language governing permissions and
16+
* limitations under the License.
17+
*/
18+
19+
'use strict';
20+
21+
// MODULES //
22+
23+
var strided2array = require( '@stdlib/array/base/from-strided' );
24+
25+
26+
// MAIN //
27+
28+
/**
29+
* Prints the contents of a workspace array for a real-valued Fourier transform.
30+
*
31+
* @private
32+
* @param {NonNegativeInteger} N - sequence length
33+
* @param {Float64Array} workspace - workspace array
34+
* @param {integer} strideW - stride length for `workspace`
35+
* @param {NonNegativeInteger} offsetW - starting index for `workspace`
36+
* @returns {void}
37+
*/
38+
function printWorkspace( N, workspace, strideW, offsetW ) {
39+
var offsetT;
40+
var offsetF;
41+
var tmp;
42+
43+
tmp = strided2array( N, workspace, strideW, offsetW );
44+
console.log( 'INTERMEDIATE RESULTS:' );
45+
console.log( tmp );
46+
console.log( ' ' );
47+
48+
offsetT = offsetW + ( N*strideW );
49+
tmp = strided2array( N, workspace, strideW, offsetT );
50+
console.log( 'TWIDDLE FACTORS:' );
51+
console.log( tmp );
52+
console.log( ' ' );
53+
54+
offsetF = offsetT + ( N*strideW );
55+
tmp = strided2array( N, workspace, strideW, offsetF );
56+
console.log( 'FACTORIZATION:' );
57+
console.log( tmp );
58+
console.log( ' ' );
59+
}
60+
61+
62+
// EXPORTS //
63+
64+
module.exports = printWorkspace;

‎lib/node_modules/@stdlib/fft/base/fftpack/lib/radb2.js‎

Lines changed: 95 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -60,69 +60,134 @@
6060

6161
'use strict';
6262

63-
// MODULES //
63+
// FUNCTIONS //
6464

65-
var chRef = require( './ch_ref.js' );
66-
var ccRef = require( './cc_ref.js' );
65+
/**
66+
* Resolves an index into the input array.
67+
*
68+
* @private
69+
* @param {NonNegativeInteger} a1 - first parameter
70+
* @param {NonNegativeInteger} a2 - second parameter
71+
* @param {NonNegativeInteger} a3 - third parameter
72+
* @param {NonNegativeInteger} ido - fourth parameter
73+
* @param {NonNegativeInteger} offset - index offset
74+
* @returns {NonNegativeInteger} index
75+
*/
76+
function iptr( a1, a2, a3, ido, offset ) { // FIXME: update with more descriptive parameter descriptions
77+
return ( ( (a3*2) + a2) * ido ) + a1 + offset; // FIXME: support strides
78+
}
79+
80+
/**
81+
* Resolves an index into the output array.
82+
*
83+
* @private
84+
* @param {NonNegativeInteger} a1 - first parameter
85+
* @param {NonNegativeInteger} a2 - second parameter
86+
* @param {NonNegativeInteger} a3 - third parameter
87+
* @param {NonNegativeInteger} l1 - fourth parameter
88+
* @param {NonNegativeInteger} ido - fifth parameter
89+
* @param {NonNegativeInteger} offset - index offset
90+
* @returns {NonNegativeInteger} index
91+
*/
92+
function optr( a1, a2, a3, l1, ido, offset ) { // FIXME: update with more descriptive parameter descriptions
93+
return ( ( (a3*l1) + a2) * ido ) + a1 + offset; // FIXME: support strides
94+
}
6795

6896

6997
// MAIN //
7098

7199
/**
72-
* Performs the backward FFT of length 2 for real-valued sequences.
100+
* Performs the backward Fourier transform of length 2 for real-valued sequences.
73101
*
74102
* @private
75-
* @param {integer} ido - number of real values for each transform
76-
* @param {integer} l1 - length of the input sequences
77-
* @param {Float64Array} cc - input array containing sequences to be transformed
78-
* @param {number} ccOffset - offset for the input array
79-
* @param {Float64Array} ch - output array containing transformed sequences
80-
* @param {number} chOffset - offset for the output array
81-
* @param {Float64Array} wa1 - array of twiddle factors
82-
* @param {number} wa1Offset - offset for the twiddle factors array
103+
* @param {NonNegativeInteger} ido - number of real values for each transform
104+
* @param {NonNegativeInteger} l1 - length of the input sequences
105+
* @param {Float64Array} cc - input array containing the sequences to be transformed
106+
* @param {NonNegativeInteger} offsetCC - offset for `cc`
107+
* @param {Float64Array} out - output array containing transformed sequences
108+
* @param {NonNegativeInteger} offsetOut - offset for `out`
109+
* @param {Float64Array} twiddles - array of twiddle factors
110+
* @param {NonNegativeInteger} offsetT - offset for `twiddles`
83111
* @returns {void}
84112
*/
85-
function radb2( ido, l1, cc, ccOffset, ch, chOffset, wa1, wa1Offset ) {
113+
function radb2( ido, l1, cc, offsetCC, out, offsetOut, twiddles, offsetT ) { // FIXME: support strides
86114
var idp2;
87-
var ti2;
88115
var tr2;
116+
var ti2;
117+
var ip1;
118+
var ip2;
119+
var it1;
120+
var it2;
121+
var io;
89122
var ic;
90123
var i;
91124
var k;
92125

93-
// Parameter adjustments...
94-
chOffset -= 1 + ( ido * ( 1 + l1 ) );
95-
ccOffset -= 1 + ( ido * 3 );
96-
wa1Offset -= 1;
126+
// Adjust offsets: // FIXME: why are these adjustments necessary? Describe here what is happening.
127+
offsetOut -= 1 + ( ido*(1+l1) );
128+
offsetCC -= 1 + ( ido*3 );
97129

98-
// Function body:
130+
// FIXME: describe what is happening below
99131
for ( k = 1; k <= l1; k++ ) {
100-
ch[ chRef( 1, k, 1, l1, ido ) + chOffset ] = cc[ ccRef( 1, 1, k, 2, ido ) + ccOffset ] + cc[ ccRef( ido, 2, k, 2, ido ) + ccOffset ];
101-
ch[ chRef( 1, k, 2, l1, ido ) + chOffset ] = cc[ ccRef( 1, 1, k, 2, ido ) + ccOffset ] - cc[ ccRef( ido, 2, k, 2, ido ) + ccOffset ];
132+
ip1 = iptr( 1, 1, k, ido, offsetCC );
133+
ip2 = iptr( ido, 2, k, ido, offsetCC );
134+
135+
io = optr( 1, k, 1, l1, ido, offsetOut );
136+
out[ io ] = cc[ ip1 ] + cc[ ip2 ];
137+
138+
io = optr( 1, k, 2, l1, ido, offsetOut );
139+
out[ io ] = cc[ ip1 ] - cc[ ip2 ];
102140
}
141+
// FIXME: describe why this check is necessary
103142
if ( ido < 2 ) {
104143
return;
105144
}
106145
if ( ido !== 2 ) {
107146
idp2 = ido + 2;
147+
148+
// FIXME: describe what is happening in the loops below
108149
for ( k = 1; k <= l1; k++ ) {
109150
for ( i = 3; i <= ido; i += 2 ) {
110151
ic = idp2 - i;
111-
ch[ chRef( i - 1, k, 1, l1, ido ) + chOffset ] = cc[ ccRef( i - 1, 1, k, 2, ido ) + ccOffset ] + cc[ ccRef( ic - 1, 2, k, 2, ido ) + ccOffset ];
112-
tr2 = cc[ ccRef( i - 1, 1, k, 2, ido ) + ccOffset ] - cc[ ccRef( ic - 1, 2, k, 2, ido ) + ccOffset ];
113-
ch[ chRef( i, k, 1, l1, ido ) + chOffset ] = cc[ ccRef( i, 1, k, 2, ido ) + ccOffset ] - cc[ ccRef( ic, 2, k, 2, ido ) + ccOffset ];
114-
ti2 = cc[ ccRef( i, 1, k, 2, ido ) + ccOffset ] + cc[ ccRef( ic, 2, k, 2, ido ) + ccOffset ];
115-
ch[ chRef( i - 1, k, 2, l1, ido ) + chOffset ] = ( wa1[ i - 2 + wa1Offset ] * tr2 ) - ( wa1[ i - 1 + wa1Offset ] * ti2 );
116-
ch[ chRef( i, k, 2, l1, ido ) + chOffset ] = ( wa1[ i - 2 + wa1Offset ] * ti2 ) + ( wa1[ i - 1 + wa1Offset ] * tr2 );
152+
153+
it1 = i - 2 + offsetT;
154+
it2 = i - 1 + offsetT;
155+
156+
ip1 = iptr( i-1, 1, k, ido, offsetCC );
157+
ip2 = iptr( ic-1, 2, k, ido, offsetCC );
158+
tr2 = cc[ ip1 ] - cc[ ip2 ];
159+
160+
io = optr( i-1, k, 1, l1, ido, offsetOut );
161+
out[ io ] = cc[ ip1 ] + cc[ ip2 ];
162+
163+
ip1 = iptr( i, 1, k, ido, offsetCC );
164+
ip2 = iptr( ic, 2, k, ido, offsetCC );
165+
ti2 = cc[ ip1 ] + cc[ ip2 ];
166+
167+
io = optr( i, k, 1, l1, ido, offsetOut );
168+
out[ io ] = cc[ ip1 ] - cc[ ip2 ];
169+
170+
io = optr( i-1, k, 2, l1, ido, offsetOut );
171+
out[ io ] = ( twiddles[ it1 ] * tr2 ) - ( twiddles[ it2 ] * ti2 );
172+
173+
io = optr( i, k, 2, l1, ido, offsetOut );
174+
out[ io ] = ( twiddles[ it1 ] * ti2 ) + ( twiddles[ it2 ] * tr2 );
117175
}
118176
}
119-
if ( ido % 2 === 1 ) {
177+
// FIXME: explain why we are checking whether `ido` is odd
178+
if ( ido%2 === 1 ) {
120179
return;
121180
}
122181
}
182+
// FIXME: describe what is happening here
123183
for ( k = 1; k <= l1; k++ ) {
124-
ch[ chRef( ido, k, 1, l1, ido ) + chOffset ] = cc[ ccRef( ido, 1, k, 2, ido ) + ccOffset ] + cc[ ccRef( ido, 1, k, 2, ido ) + ccOffset ];
125-
ch[ chRef( ido, k, 2, l1, ido ) + chOffset ] = -( cc[ ccRef( 1, 2, k, 2, ido ) + ccOffset ] + cc[ ccRef( 1, 2, k, 2, ido ) + ccOffset ] );
184+
ip1 = iptr( ido, 1, k, ido, offsetCC );
185+
io = optr( ido, k, 1, l1, ido, offsetOut );
186+
out[ io ] = cc[ ip1 ] + cc[ ip1 ];
187+
188+
ip1 = iptr( 1, 2, k, ido, offsetCC );
189+
io = optr( ido, k, 2, l1, ido, offsetOut );
190+
out[ io ] = -( cc[ ip1 ] + cc[ ip1 ] );
126191
}
127192
}
128193

0 commit comments

Comments
 (0)