Skip to content

Commit 420ff1c

Browse files
committed
feat: add radf4, radf5
1 parent 3fbfd10 commit 420ff1c

File tree

2 files changed

+362
-0
lines changed

2 files changed

+362
-0
lines changed
Lines changed: 173 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,173 @@
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+
* ## Notice
20+
*
21+
* 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.
22+
*
23+
* ```text
24+
* Copyright (c) 2004 the University Corporation for Atmospheric
25+
* Research ("UCAR"). All rights reserved. Developed by NCAR's
26+
* Computational and Information Systems Laboratory, UCAR,
27+
* www.cisl.ucar.edu.
28+
*
29+
* Redistribution and use of the Software in source and binary forms,
30+
* with or without modification, is permitted provided that the
31+
* following conditions are met:
32+
*
33+
* - Neither the names of NCAR's Computational and Information Systems
34+
* Laboratory, the University Corporation for Atmospheric Research,
35+
* nor the names of its sponsors or contributors may be used to
36+
* endorse or promote products derived from this Software without
37+
* specific prior written permission.
38+
*
39+
* - Redistributions of source code must retain the above copyright
40+
* notices, this list of conditions, and the disclaimer below.
41+
*
42+
* - Redistributions in binary form must reproduce the above copyright
43+
* notice, this list of conditions, and the disclaimer below in the
44+
* documentation and/or other materials provided with the
45+
* distribution.
46+
*
47+
* THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
48+
* EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF
49+
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
50+
* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT
51+
* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL,
52+
* EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN
53+
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
54+
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE
55+
* SOFTWARE.
56+
* ```
57+
*/
58+
59+
/* eslint-disable max-len */
60+
61+
'use strict';
62+
63+
// MODULES //
64+
65+
var chRef = require( './ch_ref.js' );
66+
var ccRef = require( './cc_ref.js' );
67+
68+
69+
// VARIABLES //
70+
71+
var hsqt2 = 0.7071067811865475;
72+
73+
74+
// MAIN //
75+
76+
/**
77+
* Performs the forward FFT of length 4 for real-valued sequences.
78+
*
79+
* @private
80+
* @param {number} ido - number of real values for each transform
81+
* @param {number} l1 - length of the input sequences
82+
* @param {Float64Array} cc - input array containing sequences to be transformed
83+
* @param {Float64Array} ch - output array containing transformed sequences
84+
* @param {Float64Array} wa1 - first array of twiddle factors
85+
* @param {Float64Array} wa2 - second array of twiddle factors
86+
* @param {Float64Array} wa3 - third array of twiddle factors
87+
* @returns {void}
88+
*/
89+
function radf4( ido, l1, cc, ch, wa1, wa2, wa3 ) {
90+
var chOffset;
91+
var ccOffset;
92+
var idp2;
93+
var tr4;
94+
var tr3;
95+
var tr2;
96+
var tr1;
97+
var ti4;
98+
var ti3;
99+
var ti2;
100+
var ti1;
101+
var cr2;
102+
var cr3;
103+
var cr4;
104+
var ci2;
105+
var ci3;
106+
var ci4;
107+
var ic;
108+
var i;
109+
var k;
110+
111+
// Parameter adjustments...
112+
chOffset = 1 + ( ido * 5 );
113+
ccOffset = 1 + ( ido * ( 1 + l1 ) );
114+
115+
// Function body:
116+
for ( k = 1; k <= l1; ++k ) {
117+
tr1 = cc[ ccRef( 1, k, 2, l1, ido ) - ccOffset ] + cc[ ccRef( 1, k, 4, l1, ido ) - ccOffset ];
118+
tr2 = cc[ ccRef( 1, k, 1, l1, ido ) - ccOffset ] + cc[ ccRef( 1, k, 3, l1, ido ) - ccOffset ];
119+
ch[ chRef( 1, 1, k, 4, ido ) - chOffset ] = tr1 + tr2;
120+
ch[ chRef( ido, 4, k, 4, ido ) - chOffset ] = tr2 - tr1;
121+
ch[ chRef( ido, 2, k, 4, ido ) - chOffset ] = cc[ ccRef( 1, k, 1, l1, ido ) - ccOffset ] - cc[ ccRef( 1, k, 3, l1, ido ) - ccOffset ];
122+
ch[ chRef( 1, 3, k, 4, ido ) - chOffset ] = cc[ ccRef( 1, k, 4, l1, ido ) - ccOffset ] - cc[ ccRef( 1, k, 2, l1, ido ) - ccOffset ];
123+
}
124+
if ( ido < 2 ) {
125+
return;
126+
}
127+
if ( ido !== 2 ) {
128+
idp2 = ido + 2;
129+
for ( k = 1; k <= l1; ++k ) {
130+
for ( i = 3; i <= ido; i += 2 ) {
131+
ic = idp2 - i;
132+
cr2 = ( wa1[ i - 2 - 1 ] * cc[ ccRef( i - 1, k, 2, l1, ido ) - ccOffset ] ) + ( wa1[ i - 1 - 1 ] * cc[ ccRef( i, k, 2, l1, ido ) - ccOffset ] );
133+
ci2 = ( wa1[ i - 2 - 1 ] * cc[ ccRef( i, k, 2, l1, ido ) - ccOffset ] ) - ( wa1[ i - 1 - 1 ] * cc[ ccRef( i - 1, k, 2, l1, ido ) - ccOffset ] );
134+
cr3 = ( wa2[ i - 2 - 1 ] * cc[ ccRef( i - 1, k, 3, l1, ido ) - ccOffset ] ) + ( wa2[ i - 1 - 1 ] * cc[ ccRef( i, k, 3, l1, ido ) - ccOffset ] );
135+
ci3 = ( wa2[ i - 2 - 1 ] * cc[ ccRef( i, k, 3, l1, ido ) - ccOffset ] ) - ( wa2[ i - 1 - 1 ] * cc[ ccRef( i - 1, k, 3, l1, ido ) - ccOffset ] );
136+
cr4 = ( wa3[ i - 2 - 1 ] * cc[ ccRef( i - 1, k, 4, l1, ido ) - ccOffset ] ) + ( wa3[ i - 1 - 1 ] * cc[ ccRef( i, k, 4, l1, ido ) - ccOffset ] );
137+
ci4 = ( wa3[ i - 2 - 1 ] * cc[ ccRef( i, k, 4, l1, ido ) - ccOffset ] ) - ( wa3[ i - 1 - 1 ] * cc[ ccRef( i - 1, k, 4, l1, ido ) - ccOffset ] );
138+
tr1 = cr2 + cr4;
139+
tr4 = cr4 - cr2;
140+
ti1 = ci2 + ci4;
141+
ti4 = ci2 - ci4;
142+
ti2 = cc[ ccRef( i, k, 1, l1, ido ) - ccOffset ] + ci3;
143+
ti3 = cc[ ccRef( i, k, 1, l1, ido ) - ccOffset ] - ci3;
144+
tr2 = cc[ ccRef( i - 1, k, 1, l1, ido ) - ccOffset ] + cr3;
145+
tr3 = cc[ ccRef( i - 1, k, 1, l1, ido ) - ccOffset ] - cr3;
146+
ch[ chRef( i - 1, 1, k, 4, ido ) - chOffset ] = tr1 + tr2;
147+
ch[ chRef( ic - 1, 4, k, 4, ido ) - chOffset ] = tr2 - tr1;
148+
ch[ chRef( i, 1, k, 4, ido ) - chOffset ] = ti1 + ti2;
149+
ch[ chRef( ic, 4, k, 4, ido ) - chOffset ] = ti1 - ti2;
150+
ch[ chRef( i - 1, 3, k, 4, ido ) - chOffset ] = ti4 + tr3;
151+
ch[ chRef( ic - 1, 2, k, 4, ido ) - chOffset ] = tr3 - ti4;
152+
ch[ chRef( i, 3, k, 4, ido ) - chOffset ] = tr4 + ti3;
153+
ch[ chRef( ic, 2, k, 4, ido ) - chOffset ] = tr4 - ti3;
154+
}
155+
}
156+
if ( ido % 2 === 1 ) {
157+
return;
158+
}
159+
}
160+
for ( k = 1; k <= l1; ++k ) {
161+
ti1 = -hsqt2 * ( cc[ ccRef( ido, k, 2, l1, ido ) - ccOffset ] + cc[ ccRef( ido, k, 4, l1, ido ) - ccOffset ] );
162+
tr1 = hsqt2 * ( cc[ ccRef( ido, k, 2, l1, ido ) - ccOffset ] - cc[ ccRef( ido, k, 4, l1, ido ) - ccOffset ] );
163+
ch[ chRef( ido, 1, k, 4, ido ) - chOffset ] = cc[ ccRef( ido, k, 1, l1, ido ) - ccOffset ] + tr1;
164+
ch[ chRef( ido, 3, k, 4, ido ) - chOffset ] = cc[ ccRef( ido, k, 1, l1, ido ) - ccOffset ] - tr1;
165+
ch[ chRef( 1, 2, k, 4, ido ) - chOffset ] = ti1 - cc[ ccRef( ido, k, 3, l1, ido ) - ccOffset ];
166+
ch[ chRef( 1, 4, k, 4, ido ) - chOffset ] = cc[ ccRef( ido, k, 3, l1, ido ) - ccOffset ] + ti1;
167+
}
168+
}
169+
170+
171+
// EXPORTS //
172+
173+
module.exports = radf4;
Lines changed: 189 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,189 @@
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+
* ## Notice
20+
*
21+
* 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.
22+
*
23+
* ```text
24+
* Copyright (c) 2004 the University Corporation for Atmospheric
25+
* Research ("UCAR"). All rights reserved. Developed by NCAR's
26+
* Computational and Information Systems Laboratory, UCAR,
27+
* www.cisl.ucar.edu.
28+
*
29+
* Redistribution and use of the Software in source and binary forms,
30+
* with or without modification, is permitted provided that the
31+
* following conditions are met:
32+
*
33+
* - Neither the names of NCAR's Computational and Information Systems
34+
* Laboratory, the University Corporation for Atmospheric Research,
35+
* nor the names of its sponsors or contributors may be used to
36+
* endorse or promote products derived from this Software without
37+
* specific prior written permission.
38+
*
39+
* - Redistributions of source code must retain the above copyright
40+
* notices, this list of conditions, and the disclaimer below.
41+
*
42+
* - Redistributions in binary form must reproduce the above copyright
43+
* notice, this list of conditions, and the disclaimer below in the
44+
* documentation and/or other materials provided with the
45+
* distribution.
46+
*
47+
* THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
48+
* EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF
49+
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
50+
* NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT
51+
* HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL,
52+
* EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN
53+
* ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
54+
* CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE
55+
* SOFTWARE.
56+
* ```
57+
*/
58+
59+
/* eslint-disable max-len */
60+
61+
'use strict';
62+
63+
// MODULES //
64+
65+
var chRef = require( './ch_ref.js' );
66+
var ccRef = require( './cc_ref.js' );
67+
68+
69+
// VARIABLES //
70+
71+
var tr11 = 0.309016994374947;
72+
var ti11 = 0.951056516295154;
73+
var tr12 = -0.809016994374947;
74+
var ti12 = 0.587785252292473;
75+
76+
77+
// MAIN //
78+
79+
/**
80+
* Performs the forward FFT of length 5 for real-valued sequences.
81+
*
82+
* @private
83+
* @param {number} ido - number of real values for each transform
84+
* @param {number} l1 - length of the input sequences
85+
* @param {Float64Array} cc - input array containing sequences to be transformed
86+
* @param {Float64Array} ch - output array containing transformed sequences
87+
* @param {Float64Array} wa1 - first array of twiddle factors
88+
* @param {Float64Array} wa2 - second array of twiddle factors
89+
* @param {Float64Array} wa3 - third array of twiddle factors
90+
* @param {Float64Array} wa4 - fourth array of twiddle factors
91+
* @returns {void}
92+
*/
93+
function radf5( ido, l1, cc, ch, wa1, wa2, wa3, wa4 ) {
94+
var chOffset;
95+
var ccOffset;
96+
var idp2;
97+
var ci2;
98+
var di2;
99+
var ci4;
100+
var ci5;
101+
var di3;
102+
var di4;
103+
var di5;
104+
var ci3;
105+
var cr2;
106+
var cr3;
107+
var dr2;
108+
var dr3;
109+
var dr4;
110+
var dr5;
111+
var cr5;
112+
var cr4;
113+
var ti2;
114+
var ti3;
115+
var ti5;
116+
var ti4;
117+
var tr2;
118+
var tr3;
119+
var tr4;
120+
var tr5;
121+
var ic;
122+
var i;
123+
var k;
124+
125+
// Parameter adjustments...
126+
chOffset = 1 + ( ido * 6 );
127+
ccOffset = 1 + ( ido * ( 1 + l1 ) );
128+
129+
// Function body:
130+
for ( k = 1; k <= l1; ++k ) {
131+
cr2 = cc[ ccRef( 1, k, 5, l1, ido ) - ccOffset ] + cc[ ccRef( 1, k, 2, l1, ido ) - ccOffset ];
132+
ci5 = cc[ ccRef( 1, k, 5, l1, ido ) - ccOffset ] - cc[ ccRef( 1, k, 2, l1, ido ) - ccOffset ];
133+
cr3 = cc[ ccRef( 1, k, 4, l1, ido ) - ccOffset ] + cc[ ccRef( 1, k, 3, l1, ido ) - ccOffset ];
134+
ci4 = cc[ ccRef( 1, k, 4, l1, ido ) - ccOffset ] - cc[ ccRef( 1, k, 3, l1, ido ) - ccOffset ];
135+
ch[ chRef( 1, 1, k, 5, ido ) - chOffset ] = cc[ ccRef( 1, k, 1, l1, ido ) - ccOffset ] + cr2 + cr3;
136+
ch[ chRef( ido, 2, k, 5, ido ) - chOffset ] = cc[ ccRef( 1, k, 1, l1, ido ) - ccOffset ] + ( tr11 * cr2 ) + ( tr12 * cr3 );
137+
ch[ chRef( 1, 3, k, 5, ido ) - chOffset ] = ( ti11 * ci5 ) + ( ti12 * ci4 );
138+
ch[ chRef( ido, 4, k, 5, ido ) - chOffset ] = cc[ ccRef( 1, k, 1, l1, ido ) - ccOffset ] + ( tr12 * cr2 ) + ( tr11 * cr3 );
139+
ch[ chRef( 1, 5, k, 5, ido ) - chOffset ] = ( ti12 * ci5 ) - ( ti11 * ci4 );
140+
}
141+
if ( ido === 1 ) {
142+
return;
143+
}
144+
idp2 = ido + 2;
145+
for ( k = 1; k <= l1; ++k ) {
146+
for ( i = 3; i <= ido; i += 2 ) {
147+
ic = idp2 - i;
148+
dr2 = ( wa1[ i - 2 - 1 ] * cc[ ccRef( i - 1, k, 2, l1, ido ) - ccOffset ] ) + ( wa1[ i - 1 - 1 ] * cc[ ccRef( i, k, 2, l1, ido ) - ccOffset ] );
149+
di2 = ( wa1[ i - 2 - 1 ] * cc[ ccRef( i, k, 2, l1, ido ) - ccOffset ] ) - ( wa1[ i - 1 - 1 ] * cc[ ccRef( i - 1, k, 2, l1, ido ) - ccOffset ] );
150+
dr3 = ( wa2[ i - 2 - 1 ] * cc[ ccRef( i - 1, k, 3, l1, ido ) - ccOffset ] ) + ( wa2[ i - 1 - 1 ] * cc[ ccRef( i, k, 3, l1, ido ) - ccOffset ] );
151+
di3 = ( wa2[ i - 2 - 1 ] * cc[ ccRef( i, k, 3, l1, ido ) - ccOffset ] ) - ( wa2[ i - 1 - 1 ] * cc[ ccRef( i - 1, k, 3, l1, ido ) - ccOffset ] );
152+
dr4 = ( wa3[ i - 2 - 1 ] * cc[ ccRef( i - 1, k, 4, l1, ido ) - ccOffset ] ) + ( wa3[ i - 1 - 1 ] * cc[ ccRef( i, k, 4, l1, ido ) - ccOffset ] );
153+
di4 = ( wa3[ i - 2 - 1 ] * cc[ ccRef( i, k, 4, l1, ido ) - ccOffset ] ) - ( wa3[ i - 1 - 1 ] * cc[ ccRef( i - 1, k, 4, l1, ido ) - ccOffset ] );
154+
dr5 = ( wa4[ i - 2 - 1 ] * cc[ ccRef( i - 1, k, 5, l1, ido ) - ccOffset ] ) + ( wa4[ i - 1 - 1 ] * cc[ ccRef( i, k, 5, l1, ido ) - ccOffset ] );
155+
di5 = ( wa4[ i - 2 - 1 ] * cc[ ccRef( i, k, 5, l1, ido ) - ccOffset ] ) - ( wa4[ i - 1 - 1 ] * cc[ ccRef( i - 1, k, 5, l1, ido ) - ccOffset ] );
156+
cr2 = dr2 + dr5;
157+
ci5 = dr5 - dr2;
158+
cr5 = di2 - di5;
159+
ci2 = di5 + di2;
160+
cr3 = dr3 + dr4;
161+
ci4 = dr4 - dr3;
162+
cr4 = di3 - di4;
163+
ci3 = di3 + di4;
164+
ch[ chRef( i - 1, 1, k, 5, ido ) - chOffset ] = cc[ ccRef( i - 1, k, 1, l1, ido ) - ccOffset ] + cr2 + cr3;
165+
ch[ chRef( i, 1, k, 5, ido ) - chOffset ] = cc[ ccRef( i, k, 1, l1, ido ) - ccOffset ] + ci2 + ci3;
166+
tr2 = cc[ ccRef( i - 1, k, 1, l1, ido ) - ccOffset ] + ( tr11 * cr2 ) + ( tr12 * cr3 );
167+
ti2 = cc[ ccRef( i, k, 1, l1, ido ) - ccOffset ] + ( tr11 * ci2 ) + ( tr12 * ci3 );
168+
tr3 = cc[ ccRef( i - 1, k, 1, l1, ido ) - ccOffset ] + ( tr12 * cr2 ) + ( tr11 * cr3 );
169+
ti3 = cc[ ccRef( i, k, 1, l1, ido ) - ccOffset ] + ( tr12 * ci2 ) + ( tr11 * ci3 );
170+
tr5 = ( ti11 * cr5 ) + ( ti12 * cr4 );
171+
ti5 = ( ti11 * ci5 ) + ( ti12 * ci4 );
172+
tr4 = ( ti12 * cr5 ) - ( ti11 * cr4 );
173+
ti4 = ( ti12 * ci5 ) - ( ti11 * ci4 );
174+
ch[ chRef( i - 1, 3, k, 5, ido ) - chOffset ] = tr2 + tr5;
175+
ch[ chRef( ic - 1, 2, k, 5, ido ) - chOffset ] = tr2 - tr5;
176+
ch[ chRef( i, 3, k, 5, ido ) - chOffset ] = ti2 + ti5;
177+
ch[ chRef( ic, 2, k, 5, ido ) - chOffset ] = ti5 - ti2;
178+
ch[ chRef( i - 1, 5, k, 5, ido ) - chOffset ] = tr3 + tr4;
179+
ch[ chRef( ic - 1, 4, k, 5, ido ) - chOffset ] = tr3 - tr4;
180+
ch[ chRef( i, 5, k, 5, ido ) - chOffset ] = ti3 + ti4;
181+
ch[ chRef( ic, 4, k, 5, ido ) - chOffset ] = ti4 - ti3;
182+
}
183+
}
184+
}
185+
186+
187+
// EXPORTS //
188+
189+
module.exports = radf5;

0 commit comments

Comments
 (0)