docs: add 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
---

Author: kgryte

lib/node_modules/@stdlib/fft/base/fftpack/lib/radf2.js

@@ -250,27 +250,20 @@ function radf2( M, L, cc, sc, oc, out, so, oo, twiddles, st, ot ) { // eslint-di
 	var i;
 	var k;
 
-	/*
+	/**
 	* First, perform the core butterfly for each sub-sequence being transformed.
 	*
-	* In the following loop, we FIXME: ????
-
-	*   Input rows  (length-M complex vectors, Re/Im interleaved):
-	*
-	*       row0 = even + odd    = [a0, b0, a1, b1, …]
-	*       row1 = even − odd    = [c0, d0, c1, d1, …]
-	*
-	* For n = 0 (elements i = 0,1) the twiddle factor is +1, so
+	* 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.
 	*
-	*   E₀ = ½(row0 + row1)   and   O₀ = ½(row0 − row1)
+	*     row0 = even + odd    (j = 0)
+	*     row1 = even - odd    (j = 1)
 	*
-	* The forward radix-2 butterfly stores two **columns**:
+	* For a radix-2 forward FFT of a real-valued signal `x` and `n = 0`,
 	*
-	*   CH(1,1,K)   ← row0 + row1           ⟶  optr( 0   , 0, k, … )
-	*   CH(IDO,2,K) ← row0 − row1           ⟶  optr( M-1 , 1, k, … )
+	*     x[0] = row0 + row1 = 2⋅even
+	*     x[M] = row0 - row1 = 2⋅odd
 	*
-	* i.e. column j = 0 holds “even + odd” starting at i = 0,
-	*      column j = 1 holds “even − odd” starting at i = M-1.
+	* 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
@@ -286,10 +279,19 @@ function radf2( M, L, cc, sc, oc, out, so, oo, twiddles, st, ot ) { // eslint-di
 	if ( M < 2 ) {
 		return;
 	}
-	/*
-	* Next, apply the general case where we need to loop through the non-trivial harmonics of each complex sub-vector.
+	/**
+	* 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
 	*
-	* // FIXME: complete extended description
+	* 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;
@@ -298,49 +300,63 @@ function radf2( M, L, cc, sc, oc, out, so, oo, twiddles, st, ot ) { // eslint-di
 		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 => ????
+				im = MP1 - i; // "mirror" index
 
-				it1 = ( (i-1)*st ) + ot; // ????
-				it2 = ( i*st ) + ot;     // ????
+				// Resolve twiddle factor indices:
+				it1 = ( (i-1)*st ) + ot; // cos(θ)
+				it2 = ( i*st ) + ot;     // sin(θ)
 
-				ip1 = iptr( i, k, 1, L, M, sc, oc ); // ????
-				ip2 = iptr( i+1, k, 1, L, M, sc, oc ); // ????
-				tr2 = ( twiddles[ it1 ] * cc[ ip1 ] ) + ( twiddles[ it2 ] * cc[ ip2 ] ); // ????
-				ti2 = ( twiddles[ it1 ] * cc[ ip2 ] ) - ( twiddles[ it2 ] * cc[ ip1 ] ); // ????
+				// 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)
 
-				ip1 = iptr( i+1, k, 0, L, M, sc, oc ); // ????
-				io = optr( i+1, 0, k, M, so, oo ); // ????
-				out[ io ] = cc[ ip1 ] + ti2; // ????
+				// 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)
 
-				io = optr( im, 1, k, M, so, oo ); // ????
-				out[ io ] = ti2 - cc[ ip1 ]; // ????
+				// 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;
 
-				ip1 = iptr( i, k, 0, L, M, sc, oc ); // ????
-				io = optr( i, 0, k, M, so, oo ); // ????
-				out[ io ] = cc[ ip1 ] + tr2; // ????
+				io = optr( im, 1, k, M, so, oo );      // Im(x[2n+1])
+				out[ io ] = ti2 - cc[ ip1 ];
 
-				io = optr( im-1, 1, k, M, so, oo ); // ????
-				out[ io ] = cc[ ip1 ] - tr2; // ????
+				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, and, thus, the central element is purely real and was handled in the very first butterfly, so we do not need to perform any further transformations...
+		// 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, FIXME: ????
+	* 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 );
+		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 ];
+		out[ io ] = -cc[ ip1 ];                // -(-Re(O_n)) = Re(O_n)
 
-		ip1 = iptr( M-1, k, 0, L, M, sc, oc );
+		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 ];
 	}