Skip to content

Commit 5d5bb6a

Browse files
kgrytekgryte
committed
refactor: initial stride support and additional comments
--- 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: na - task: lint_typescript_tests status: na - task: lint_license_headers status: passed ---
1 parent 72a102e commit 5d5bb6a

File tree

1 file changed

+122
-66
lines changed
  • lib/node_modules/@stdlib/fft/base/fftpack/lib

1 file changed

+122
-66
lines changed

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

Lines changed: 122 additions & 66 deletions
Original file line numberDiff line numberDiff line change
@@ -63,54 +63,97 @@
6363
// FUNCTIONS //
6464

6565
/**
66-
* Resolves an index into the input array.
66+
* Resolves an index in the input array.
67+
*
68+
* ## Notes
69+
*
70+
* When transforming 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 "columns" corresponding to the real and imaginary parts of a folded complex vector and where each "column" has `M` elements.
71+
*
72+
* Accordingly, the following is a logical view of an input array (zero-based indexing) which contains `L = 3` transforms and in which each "column" sub-sequence has length `M = 4`:
73+
*
74+
* ```text
75+
* j = 0 ("even" column) j = 1 ("odd" column)
76+
* k = 0 ─┬────────────────────────────────────┬────────────────────────────────────┐
77+
* │ cc(0,0,0) cc(1,0,0) ... cc(3,0,0) │ cc(0,1,0) cc(1,1,0) ... cc(3,1,0) │
78+
* └────────────────────────────────────┴────────────────────────────────────┤
79+
* k = 1 ─┬────────────────────────────────────┬────────────────────────────────────┤
80+
* │ cc(0,0,1) cc(1,0,1) ... cc(3,0,1) │ cc(0,1,1) cc(1,1,1) ... cc(3,1,1) │
81+
* └────────────────────────────────────┴────────────────────────────────────┤
82+
* k = 2 ─┬────────────────────────────────────┬────────────────────────────────────┤
83+
* │ cc(0,0,2) cc(1,0,2) ... cc(3,0,2) │ cc(0,1,2) cc(1,1,2) ... cc(3,1,2) │
84+
* └────────────────────────────────────┴────────────────────────────────────┘
85+
* ↑ ↑ ↑ ↑ ↑ ↑
86+
* i = 0 1 M-1 0 1 M-1
87+
* ```
88+
*
89+
* In the above,
90+
*
91+
* - `i` is the fastest varying index, which walks within one short "column" sub-sequence.
92+
* - `j` selects between the even and odd half-spectra (i.e., the real or imaginary part of a folded complex vector).
93+
* - `k` specifies the index of one of the `L` independent transforms we are processing.
94+
*
95+
* In linear memory, the three-dimensional logical view is arranged as follows:
96+
*
97+
* ```text
98+
* | cc(0,0,0) ... cc(3,0,0) | cc(0,1,0) ... cc(3,1,0) | cc(0,0,1) ... cc(3,0,1) | ... | cc(0,1,2) ... cc(3,1,2) |
99+
* ```
67100
*
68101
* @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
102+
* @param {NonNegativeInteger} i - index of an element within a sub-sequence
103+
* @param {NonNegativeInteger} j - index specifying which of the two complex halves we are in (either `0` or `1`)
104+
* @param {NonNegativeInteger} k - index of the sub-sequence being transformed
105+
* @param {NonNegativeInteger} M - sub-sequence length
106+
* @param {integer} stride - stride length of the input array
107+
* @param {NonNegativeInteger} offset - index specifying the first indexed element in the input array
108+
* @returns {NonNegativeInteger} computed index
75109
*/
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
110+
function iptr( i, j, k, M, stride, offset ) {
111+
var n = i + ( ( j+(k*2) ) * M );
112+
return ( n*stride ) + offset;
78113
}
79114

80115
/**
81-
* Resolves an index into the output array.
116+
* Resolves an index in the output array.
117+
*
118+
* ## Notes
119+
*
120+
* 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 real and imaginary parts of a folded complex vector and where each "column" has `M` elements.
121+
*
122+
* Accordingly, the following is a logical view of an input array (zero-based indexing) which contains `L = 3` transforms and in which each "column" sub-sequence has length `M = 4`:
82123
*
83124
* @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
125+
* @param {NonNegativeInteger} i - index of an element within a sub-sequence
126+
* @param {NonNegativeInteger} k - index of the sub-sequence being transformed
127+
* @param {NonNegativeInteger} j - output row
128+
* @param {NonNegativeInteger} L - number of sub-sequences
129+
* @param {NonNegativeInteger} M - sub-sequence length
130+
* @param {integer} stride - stride length of the output array
131+
* @param {NonNegativeInteger} offset - index specifying the first indexed element in the output array
132+
* @returns {NonNegativeInteger} computed index
91133
*/
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
134+
function optr( i, k, j, L, M, stride, offset ) {
135+
var n = i + ( ( k+(j*L) ) * M );
136+
return ( n*stride ) + offset;
94137
}
95138

96139

97140
// MAIN //
98141

99142
/**
100-
* Performs the backward Fourier transform of length 2 for real-valued sequences.
143+
* Performs the backward Fourier radix-2 transform for a real-valued sequence.
101144
*
102145
* @private
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`
146+
* @param {NonNegativeInteger} M - number of elements in each sub-sequence to be transformed
147+
* @param {NonNegativeInteger} L - number of sub-sequences to be transformed
148+
* @param {Float64Array} cc - input array containing the sub-sequences to be transformed
149+
* @param {NonNegativeInteger} oc - offset for `cc`
150+
* @param {Float64Array} out - output array containing transformed sub-sequences
151+
* @param {NonNegativeInteger} oo - offset for `out`
109152
* @param {Float64Array} twiddles - array of twiddle factors
110-
* @param {NonNegativeInteger} offsetT - offset for `twiddles`
153+
* @param {NonNegativeInteger} ot - offset for `twiddles`
111154
* @returns {void}
112155
*/
113-
function radb2( ido, l1, cc, offsetCC, out, offsetOut, twiddles, offsetT ) { // FIXME: support strides
156+
function radb2( M, L, cc, oc, out, oo, twiddles, ot ) { // FIXME: support strides
114157
var idp2;
115158
var tr2;
116159
var ti2;
@@ -120,73 +163,86 @@ function radb2( ido, l1, cc, offsetCC, out, offsetOut, twiddles, offsetT ) { //
120163
var it2;
121164
var io;
122165
var ic;
166+
var sc = 1;
167+
var so = 1;
123168
var i;
124-
var k;
125-
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 );
129-
130-
// FIXME: describe what is happening below
131-
for ( k = 1; k <= l1; k++ ) {
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 ];
169+
var k; // FIXME: make a parameter
170+
171+
/*
172+
* Perform the core butterfly for each sub-sequence being transformed.
173+
*
174+
* In the following loop, we combine two half-length complex vectors.
175+
*
176+
* cc(:, 0, k) holds Re(X_even) + j Im(X_even)
177+
* cc(:, 1, k) holds Re(X_odd) + j Im(X_odd)
178+
*
179+
* For a radix-2 backward FFT, the twiddle factor for the first column is +1, so we can compute to real outputs with just
180+
*
181+
* out(even) = even + odd
182+
* out(odd) = even - odd
183+
*/
184+
for ( k = 0; k < L; k++ ) {
185+
ip1 = iptr( 0, 0, k, M, sc, oc );
186+
ip2 = iptr( M-1, 1, k, M, sc, oc );
187+
188+
io = optr( 0, k, 0, L, M, so, oo );
189+
out[ io ] = cc[ ip1 ] + cc[ ip2 ]; // out(even) = even + odd
190+
191+
io = optr( 0, k, 1, L, M, so, oo );
192+
out[ io ] = cc[ ip1 ] - cc[ ip2 ]; // out(odd) = even - odd
140193
}
141194
// FIXME: describe why this check is necessary
142-
if ( ido < 2 ) {
195+
if ( M < 2 ) {
143196
return;
144197
}
145-
if ( ido !== 2 ) {
146-
idp2 = ido + 2;
198+
199+
// FIXME: everything below needs updating as `iptr` and `optr` have been converted to zero-based indexing
200+
201+
if ( M !== 2 ) {
202+
idp2 = M + 2;
147203

148204
// FIXME: describe what is happening in the loops below
149-
for ( k = 1; k <= l1; k++ ) {
150-
for ( i = 3; i <= ido; i += 2 ) {
205+
for ( k = 0; k < L; k++ ) {
206+
for ( i = 3; i <= M; i += 2 ) {
151207
ic = idp2 - i;
152208

153-
it1 = i - 2 + offsetT;
154-
it2 = i - 1 + offsetT;
209+
it1 = i - 2 + ot;
210+
it2 = i - 1 + ot;
155211

156-
ip1 = iptr( i-1, 1, k, ido, offsetCC );
157-
ip2 = iptr( ic-1, 2, k, ido, offsetCC );
212+
ip1 = iptr( i-1, 1, k, M, oc );
213+
ip2 = iptr( ic-1, 2, k, M, oc );
158214
tr2 = cc[ ip1 ] - cc[ ip2 ];
159215

160-
io = optr( i-1, k, 1, l1, ido, offsetOut );
216+
io = optr( i-1, k, 1, L, M, oo );
161217
out[ io ] = cc[ ip1 ] + cc[ ip2 ];
162218

163-
ip1 = iptr( i, 1, k, ido, offsetCC );
164-
ip2 = iptr( ic, 2, k, ido, offsetCC );
219+
ip1 = iptr( i, 1, k, M, oc );
220+
ip2 = iptr( ic, 2, k, M, oc );
165221
ti2 = cc[ ip1 ] + cc[ ip2 ];
166222

167-
io = optr( i, k, 1, l1, ido, offsetOut );
223+
io = optr( i, k, 1, L, M, oo );
168224
out[ io ] = cc[ ip1 ] - cc[ ip2 ];
169225

170-
io = optr( i-1, k, 2, l1, ido, offsetOut );
226+
io = optr( i-1, k, 2, L, M, oo );
171227
out[ io ] = ( twiddles[ it1 ] * tr2 ) - ( twiddles[ it2 ] * ti2 );
172228

173-
io = optr( i, k, 2, l1, ido, offsetOut );
229+
io = optr( i, k, 2, L, M, oo );
174230
out[ io ] = ( twiddles[ it1 ] * ti2 ) + ( twiddles[ it2 ] * tr2 );
175231
}
176232
}
177-
// FIXME: explain why we are checking whether `ido` is odd
178-
if ( ido%2 === 1 ) {
233+
// FIXME: explain why we are checking whether `M` is odd
234+
if ( M%2 === 1 ) {
179235
return;
180236
}
181237
}
182238
// FIXME: describe what is happening here
183-
for ( k = 1; k <= l1; k++ ) {
184-
ip1 = iptr( ido, 1, k, ido, offsetCC );
185-
io = optr( ido, k, 1, l1, ido, offsetOut );
239+
for ( k = 0; k < L; k++ ) {
240+
ip1 = iptr( M, 1, k, M, oc );
241+
io = optr( M, k, 1, L, M, oo );
186242
out[ io ] = cc[ ip1 ] + cc[ ip1 ];
187243

188-
ip1 = iptr( 1, 2, k, ido, offsetCC );
189-
io = optr( ido, k, 2, l1, ido, offsetOut );
244+
ip1 = iptr( 1, 2, k, M, oc );
245+
io = optr( M, k, 2, L, M, oo );
190246
out[ io ] = -( cc[ ip1 ] + cc[ ip1 ] );
191247
}
192248
}

0 commit comments

Comments
 (0)