NCEPLIBS-ip 5.3.0
•All Data Structures Namespaces Files Functions Variables Pages
fftpack.F
Go to the documentation of this file.
1C> @file
2C> @brief A concatenation of several subroutines from the
3C> <a target="_blank" href="https://netlib.org/fftpack/">FFTPACK</a> collection.
4C> FFTPACK is a package of Fortran subprograms for the fast Fourier
5C> transform of periodic and other symmetric sequences. It includes
6C> complex, real, sine, cosine, and quarter-wave transforms.
7C>
8C> See
9C> <a target="_blank" href="https://netlib.org/fftpack/doc">FFTPACK documentation</a>,
10C> <a target="_blank" href="https://help.graphica.com.au/irix-6.5.30/man/3S/dzfft">IRIX man pages</a>, and
11C> <a target="_blank" href="https://www.ibm.com/docs/en/essl/6.2?topic=subroutines-scrft-dcrft-complex-real-fourier-transform">IBM ESSL documentation</a>
12C> for further historical context for various subroutines included in this file.
13C>
14C> Reference:
15C> - P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
16C> (G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83.
17C>
18C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
19
20C> @brief Computes a set of `m` real discrete `n`-point Fourier transforms of complex conjugate even data
21C>
22C> The 'd' stands for double precision, but this is historical and may be ignored,
23C> and this subroutine is identical to `scrft()`.
24C>
25C> @param init Initialization flag; if nonzero, initializes the transformation
26C> @param x Input array of real values, dimensions (`2*ldx`, `m`)
27C> @param ldx Leading dimension of the input array `x`
28C> @param y Output array storing the transformed values, dimensions (`ldy`, `m`)
29C> @param ldy Leading dimension of the output array `y`
30C> @param n Number of data points in each transform (first dimension of `y`)
31C> @param m Number of transforms to be computed (second dimension of `y`)
32C> @param isign Sign indicator for the transform direction (not explicitly used)
33C> @param scale Scaling factor applied to the transformed values after computation
34C> @param table Work array of size 44002, precomputed during initialization
35C> @param n1 Auxiliary parameter
36C> @param wrk Work array used for intermediate storage during computation
37C> @param n2 Auxiliary parameter
38C> @param z Optional dummy parameter for compatibility
39C> @param nz Optional dummy parameter for compatibility
40C>
41C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
42 SUBROUTINE dcrft(init,x,ldx,y,ldy,n,m,isign,scale,
43 & table,n1,wrk,n2,z,nz)
44
45 implicit none
46 integer init,ldx,ldy,n,m,isign,n1,n2,i,j
47 real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk(*)
48 real, optional :: z
49 integer, optional :: nz
50
51 IF (init.ne.0) THEN
52 CALL rffti(n,table)
53 ELSE
54!OCL NOVREC
55 DO j=1,m
56 y(1,j)=x(1,j)
57 DO i=2,n
58 y(i,j)=x(i+1,j)
59 ENDDO
60 CALL rfftb(n,y(1,j),table)
61 DO i=1,n
62 y(i,j)=scale*y(i,j)
63 ENDDO
64 ENDDO
65 ENDIF
66
67 RETURN
68 END
69
70C> @brief Compute a set of m real discrete n-point Fourier transforms of complex conjugate even data
71C>
72C> The 's' stands for single precision, but this is historical and may be ignored,
73C> and this subroutine is identical to `dcrft()`.
74C>
75C> @param init Initialization flag; if nonzero, initializes the transformation
76C> @param x Input array of real values, dimensions (`2*ldx`, `m`)
77C> @param ldx Leading dimension of the input array `x`
78C> @param y Output array storing the transformed values, dimensions (`ldy`, `m`)
79C> @param ldy Leading dimension of the output array `y`
80C> @param n Number of data points in each transform (first dimension of `y`)
81C> @param m Number of transforms to be computed (second dimension of `y`)
82C> @param isign Sign indicator for the transform direction (not explicitly used)
83C> @param scale Scaling factor applied to the transformed values after computation
84C> @param table Work array of size 44002, precomputed during initialization
85C> @param n1 Auxiliary parameter
86C> @param wrk Work array used for intermediate storage during computation
87C> @param n2 Auxiliary parameter
88C> @param z Optional dummy parameter for compatibility
89C> @param nz Optional dummy parameter for compatibility
90C>
91C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
92 SUBROUTINE scrft(init,x,ldx,y,ldy,n,m,isign,scale,
93 & table,n1,wrk,n2,z,nz)
94
95 implicit none
96 integer init,ldx,ldy,n,m,isign,n1,n2,i,j
97 real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk(*)
98 real, optional :: z
99 integer, optional :: nz
100
101 IF (init.ne.0) THEN
102 CALL rffti(n,table)
103 ELSE
104!OCL NOVREC
105 DO j=1,m
106 y(1,j)=x(1,j)
107 DO i=2,n
108 y(i,j)=x(i+1,j)
109 ENDDO
110 CALL rfftb(n,y(1,j),table)
111 DO i=1,n
112 y(i,j)=scale*y(i,j)
113 ENDDO
114 ENDDO
115 ENDIF
116
117 RETURN
118 END
119
120C> @brief Compute a complex discrete Fourier transform of real data
121C>
122C> This subroutine performs a Fourier transform on real input data.
123C> It supports initialization of the trigonometric coefficient
124C> table and forward or inverse transformations. The `isign` parameter determines
125C> the operation mode.
126C>
127C> @param isign Operation mode: 0 for initialization, 1 for inverse transform
128C> @param n Number of data points in the transform
129C> @param scale Scaling factor applied to the transformed values
130C> @param x Input array of real values, length at least `n+1`
131C> @param y Output array storing the transformed values, length at least `n`
132C> @param table Work array used for storing trigonometric coefficients
133C> @param work Work array used for intermediate computations
134C> @param isys Dummy parameter
135C>
136C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
137 SUBROUTINE csfft(isign,n,scale,x,y,table,work,isys)
138
139 implicit none
140 integer isign,n,isys,i
141 real scale,x(*),y(*),table(*),work(*)
142
143 IF (isign.eq.0) THEN
144 CALL rffti(n,table)
145 ENDIF
146 IF (isign.eq.1) THEN
147 y(1)=x(1)
148 DO i=2,n
149 y(i)=x(i+1)
150 ENDDO
151 CALL rfftb(n,y,table)
152 DO i=1,n
153 y(i)=scale*y(i)
154 ENDDO
155 ENDIF
156
157 RETURN
158 END
159
160C> @brief Compute a set of m complex discrete n-point Fourier transforms of real data
161C>
162C> This subroutine performs multiple (`m`) discrete Fourier transforms of length `n`
163C> on real input data. The transforms are computed using the real-to-complex Fast Fourier
164C> Transform (FFT) approach. The subroutine requires two invocations: once with `init` set
165C> to a non-zero value to initialize work storage, the second with `init=0` to perform the DFT.
166C>
167C> The 'd' stands for double precision, but this is historical and may be ignored,
168C> and this subroutine is identical to `srcft()`.
169C>
170C> @param init Initialization flag; if nonzero, initializes the transformation
171C> @param x Input array of real values, dimensions (`ldx`, `m`)
172C> @param ldx Leading dimension of the input array `x`
173C> @param y Output array storing the transformed values, dimensions (`2*ldy`, `m`)
174C> @param ldy Leading dimension of the output array `y`
175C> @param n Number of data points in each transform (first dimension of `y`)
176C> @param m Number of independent transforms to be computed (second dimension of `y`)
177C> @param isign Sign indicator for the transform direction (not explicitly used)
178C> @param scale Scaling factor applied to the transformed values
179C> @param table Work array of size 44002
180C> @param n1 Auxiliary parameter
181C> @param wrk Work array used for intermediate storage during computation
182C> @param n2 Auxiliary parameter
183C> @param z Optional dummy parameter for compatibility
184C> @param nz Optional dummy parameter for compatibility
185C>
186C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
187 SUBROUTINE drcft(init,x,ldx,y,ldy,n,m,isign,scale,
188 & table,n1,wrk,n2,z,nz)
189
190 implicit none
191 integer init,ldx,ldy,n,m,isign,n1,n2,i,j
192 real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk(*)
193 real, optional :: z
194 integer, optional :: nz
195
196 IF (init.ne.0) THEN
197 CALL rffti(n,table)
198 ELSE
199 DO j=1,m
200 DO i=1,n
201 y(i,j)=x(i,j)
202 ENDDO
203 CALL rfftf(n,y(1,j),table)
204 DO i=1,n
205 y(i,j)=scale*y(i,j)
206 ENDDO
207 DO i=n,2,-1
208 y(i+1,j)=y(i,j)
209 ENDDO
210 y(2,j)=0.
211C 01/17/2013 vvvvvvvvvvvvv E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
212 y(n+2,j) = 0.
213 ENDDO
214 ENDIF
215
216 RETURN
217 END
218
219C> @brief Compute a set of m complex discrete n-point Fourier transforms of real data
220C>
221C> This subroutine performs multiple (`m`) discrete Fourier transforms of length `n`
222C> on real input data. The transforms are computed using the real-to-complex Fast Fourier
223C> Transform (FFT) approach. The subroutine requires two invocations: once with `init` set
224C> to a non-zero value to initialize work storage, the second with `init=0` to perform the DFT.
225C>
226C> The 's' stands for single precision, but this is historical and may be ignored,
227C> and this subroutine is identical to `drcft()`.
228C>
229C> @param init Initialization flag; if nonzero, initializes the transformation
230C> @param x Input array of real values, dimensions (`ldx`, `m`)
231C> @param ldx Leading dimension of the input array `x`
232C> @param y Output array storing the transformed values, dimensions (`2*ldy`, `m`)
233C> @param ldy Leading dimension of the output array `y`
234C> @param n Number of data points in each transform (first dimension of `y`)
235C> @param m Number of independent transforms to be computed (second dimension of `y`)
236C> @param isign Sign indicator for the transform direction (not explicitly used)
237C> @param scale Scaling factor applied to the transformed values
238C> @param table Work array of size 44002
239C> @param n1 Auxiliary parameter
240C> @param wrk Work array used for intermediate storage during computation
241C> @param n2 Auxiliary parameter
242C> @param z Optional dummy parameter for compatibility
243C> @param nz Optional dummy parameter for compatibility
244C>
245C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
246 SUBROUTINE srcft(init,x,ldx,y,ldy,n,m,isign,scale,
247 & table,n1,wrk,n2,z,nz)
248
249 implicit none
250 integer init,ldx,ldy,n,m,isign,n1,n2,i,j
251 real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk(*)
252 real, optional :: z
253 integer, optional :: nz
254
255 IF (init.ne.0) THEN
256 CALL rffti(n,table)
257 ELSE
258 DO j=1,m
259 DO i=1,n
260 y(i,j)=x(i,j)
261 ENDDO
262 CALL rfftf(n,y(1,j),table)
263 DO i=1,n
264 y(i,j)=scale*y(i,j)
265 ENDDO
266 DO i=n,2,-1
267 y(i+1,j)=y(i,j)
268 ENDDO
269 y(2,j)=0.
270 y(n+2,j) = 0.
271C 01/17/2013 ^^^^^^^^^^E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
272 ENDDO
273 ENDIF
274
275 RETURN
276 END
277
278C> @brief Compute a real-to-complex discrete Fourier transform
279C>
280C> This subroutine performs a Fourier transform on real input data.
281C> It supports initialization of the trigonometric coefficient table and forward transformation.
282C>
283C> @param isign Operation mode: 0 for initialization, -1 for forward transform
284C> @param n Number of data points in the transform
285C> @param scale Scaling factor applied to the transformed values
286C> @param x Input array of real values, length at least `n`
287C> @param y Output array storing the transformed values, length at least `n+1`
288C> @param table Work array used for storing trigonometric coefficients
289C> @param work Work array used for intermediate computations
290C> @param isys Dummy parameter
291C>
292C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
293 SUBROUTINE scfft(isign,n,scale,x,y,table,work,isys)
294
295 implicit none
296 integer isign,n,isys,i
297 real scale,x(*),y(*),table(*),work(*)
298
299 IF (isign.eq.0) THEN
300 CALL rffti(n,table)
301 ENDIF
302 IF (isign.eq.-1) THEN
303 DO i=1,n
304 y(i)=x(i)
305 ENDDO
306 CALL rfftf(n,y,table)
307 DO i=1,n
308 y(i)=scale*y(i)
309 ENDDO
310 DO i=n,2,-1
311 y(i+1)=y(i)
312 ENDDO
313 y(2)=0.
314 ENDIF
315
316 RETURN
317 END
318
319C> @brief Compute a forward transform of a real periodic sequence
320C>
321C> The wsave array must be initialized by calling subroutine rffti(n,wsave) and a
322C> different wsave array must be used for each different value of n. This
323C> initialization does not have to be repeated so long as n remains unchanged thus
324C> subsequent transforms can be obtained faster than the first.
325C>
326C> @note This transform is unnormalized since a call of `rfftf()` followed by a
327C> call of `rfftb()` will multiply the input sequence by n.
328C>
329C> @param n The length of the array `r` to be transformed
330C> @param r A real array of length `n` which contains the sequence to be transformed
331C> @param wsave A work array which must be dimensioned at least `2*n+15`; must not be destroyed between calls of `rfftf()` or `rfftb()`
332C>
333C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
334 SUBROUTINE rfftf (N,R,WSAVE)
335 dimension r(*) ,wsave(*)
336 IF (n .EQ. 1) RETURN
337 CALL rfftf1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
338 RETURN
339 END
340
341C> @brief Compute the real periodic sequence from its Fourier coefs
342C>
343C> @param n The length of the array `r` to be transformed
344C> @param r A real array of length `n` which contains the sequence to be transformed
345C> @param wsave A work array which must be dimensioned at least `2*n+15`; must not be destroyed between calls of `rfftf()` or `rfftb()`
346C>
347C> @note This transform is unnormalized since a call of `rfftf()` followed by a
348C> call of `rfftb()` will multiply the input sequence by n.
349C>
350C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
351 SUBROUTINE rfftb (N,R,WSAVE)
352 dimension r(*) ,wsave(*)
353 IF (n .EQ. 1) RETURN
354 CALL rfftb1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
355 RETURN
356 END
357
358C> @brief Initialize the array wsave which is used in both `rfftf()` and `rfftb()`.
359C>
360C> @param n The length of the sequence to be transformed
361C> @param wsave A work array which must be dimensioned at least `2*n+15`; must not be destroyed between calls of `rfftf()` or `rfftb()`
362C>
363C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
364 SUBROUTINE rffti (N,WSAVE)
365 dimension wsave(*)
366 IF (n .EQ. 1) RETURN
367 CALL rffti1 (n,wsave(n+1),wsave(2*n+1))
368 RETURN
369 END
370
371C> @brief Compute the inverse fast Fourier transform (IFFT) using a mixed-radix algorithm
372C>
373C> This subroutine is a low-level component of the real-to-complex IFFT computation.
374C> It applies a series of radix-based backward FFTs
375C> (`radb2()`, `radb3()`, `radb4()`, `radb5()`, and `radbg()`)
376C> to transform real input data back from the frequency domain to the time domain.
377C>
378C> @param n Number of data points in the transform
379C> @param c Input/output array of real values; on input, it contains transformed data,
380C> and on output, it stores the inverse-transformed real data
381C> @param ch Work array used for intermediate calculations
382C> @param wa Sine and cosine trigonometric table used for the transformation
383C> @param ifac Integer factorization array used to determine the radix decomposition
384C>
385C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
386 SUBROUTINE rfftb1 (N,C,CH,WA,IFAC)
387 REAL CH(*) ,C(*) ,WA(*) ,IFAC(*)
388 NF = int(ifac(2))
389 na = 0
390 l1 = 1
391 iw = 1
392 DO 116 k1=1,nf
393 ip = int(ifac(k1+2))
394 l2 = ip*l1
395 ido = n/l2
396 idl1 = ido*l1
397 IF (ip .NE. 4) GO TO 103
398 ix2 = iw+ido
399 ix3 = ix2+ido
400 IF (na .NE. 0) GO TO 101
401 CALL radb4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
402 GO TO 102
403 101 CALL radb4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
404 102 na = 1-na
405 GO TO 115
406 103 IF (ip .NE. 2) GO TO 106
407 IF (na .NE. 0) GO TO 104
408 CALL radb2 (ido,l1,c,ch,wa(iw))
409 GO TO 105
410 104 CALL radb2 (ido,l1,ch,c,wa(iw))
411 105 na = 1-na
412 GO TO 115
413 106 IF (ip .NE. 3) GO TO 109
414 ix2 = iw+ido
415 IF (na .NE. 0) GO TO 107
416 CALL radb3 (ido,l1,c,ch,wa(iw),wa(ix2))
417 GO TO 108
418 107 CALL radb3 (ido,l1,ch,c,wa(iw),wa(ix2))
419 108 na = 1-na
420 GO TO 115
421 109 IF (ip .NE. 5) GO TO 112
422 ix2 = iw+ido
423 ix3 = ix2+ido
424 ix4 = ix3+ido
425 IF (na .NE. 0) GO TO 110
426 CALL radb5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
427 GO TO 111
428 110 CALL radb5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
429 111 na = 1-na
430 GO TO 115
431 112 IF (na .NE. 0) GO TO 113
432 CALL radbg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
433 GO TO 114
434 113 CALL radbg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
435 114 IF (ido .EQ. 1) na = 1-na
436 115 l1 = l2
437 iw = iw+(ip-1)*ido
438 116 CONTINUE
439 IF (na .EQ. 0) RETURN
440 DO 117 i=1,n
441 c(i) = ch(i)
442 117 CONTINUE
443 RETURN
444 END
445
446C> @brief Compute the real forward fast Fourier transform (FFT) using a mixed-radix algorithm
447C>
448C> This subroutine performs the FFT by applying a sequence of radix-based forward
449C> transformations (`radf2()`, `radf3()`, `radf4()`, `radf5()`, and `radfg()`).
450C> It is a low-level routine used internally for computing the real-input FFT.
451C>
452C> @param n Number of data points in the transform
453C> @param c Input/output array of real values; on input, it contains time-domain data,
454C> and on output, it stores the transformed frequency-domain data
455C> @param ch Work array used for intermediate calculations
456C> @param wa Sine and cosine trigonometric table used for the transformation
457C> @param ifac Integer factorization array used to determine the radix decomposition
458C>
459C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
460 SUBROUTINE rfftf1 (N,C,CH,WA,IFAC)
461 REAL CH(*) ,C(*) ,WA(*) ,IFAC(*)
462 NF = int(ifac(2))
463 na = 1
464 l2 = n
465 iw = n
466 DO 111 k1=1,nf
467 kh = nf-k1
468 ip = int(ifac(kh+3))
469 l1 = l2/ip
470 ido = n/l2
471 idl1 = ido*l1
472 iw = iw-(ip-1)*ido
473 na = 1-na
474 IF (ip .NE. 4) GO TO 102
475 ix2 = iw+ido
476 ix3 = ix2+ido
477 IF (na .NE. 0) GO TO 101
478 CALL radf4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
479 GO TO 110
480 101 CALL radf4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
481 GO TO 110
482 102 IF (ip .NE. 2) GO TO 104
483 IF (na .NE. 0) GO TO 103
484 CALL radf2 (ido,l1,c,ch,wa(iw))
485 GO TO 110
486 103 CALL radf2 (ido,l1,ch,c,wa(iw))
487 GO TO 110
488 104 IF (ip .NE. 3) GO TO 106
489 ix2 = iw+ido
490 IF (na .NE. 0) GO TO 105
491 CALL radf3 (ido,l1,c,ch,wa(iw),wa(ix2))
492 GO TO 110
493 105 CALL radf3 (ido,l1,ch,c,wa(iw),wa(ix2))
494 GO TO 110
495 106 IF (ip .NE. 5) GO TO 108
496 ix2 = iw+ido
497 ix3 = ix2+ido
498 ix4 = ix3+ido
499 IF (na .NE. 0) GO TO 107
500 CALL radf5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
501 GO TO 110
502 107 CALL radf5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
503 GO TO 110
504 108 IF (ido .EQ. 1) na = 1-na
505 IF (na .NE. 0) GO TO 109
506 CALL radfg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
507 na = 1
508 GO TO 110
509 109 CALL radfg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
510 na = 0
511 110 l2 = l1
512 111 CONTINUE
513 IF (na .EQ. 1) RETURN
514 DO 112 i=1,n
515 c(i) = ch(i)
516 112 CONTINUE
517 RETURN
518 END
519
520C> @brief Compute the initialization factors for the real-input FFT
521C>
522C> This subroutine prepares the trigonometric table and factorization
523C> for the real-input Fast Fourier Transform (FFT). It determines the
524C> prime factorization of N and computes the necessary twiddle factors
525C> for efficient computation.
526C>
527C> @param n The length of the input data array
528C> @param wa Output array containing computed trigonometric factors
529C> @param ifac Output integer array storing the factorization of `n` for FFT processing
530C>
531C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
532 SUBROUTINE rffti1 (N,WA,IFAC)
533 REAL WA(*) ,IFAC(*) ,NTRYH(4)
534 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
535 NL = n
536 nf = 0
537 j = 0
538 101 j = j+1
539 IF ((j-4).LE.0) THEN
540 GO TO 102
541 ELSE
542 GO TO 103
543 ENDIF
544 102 ntry = int(ntryh(j))
545 GO TO 104
546 103 ntry = ntry+2
547 104 nq = nl/ntry
548 nr = nl-ntry*nq
549 IF (nr.EQ.0) THEN
550 GO TO 105
551 ELSE
552 GO TO 101
553 ENDIF
554 105 nf = nf+1
555 ifac(nf+2) = ntry
556 nl = nq
557 IF (ntry .NE. 2) GO TO 107
558 IF (nf .EQ. 1) GO TO 107
559 DO 106 i=2,nf
560 ib = nf-i+2
561 ifac(ib+2) = ifac(ib+1)
562 106 CONTINUE
563 ifac(3) = 2
564 107 IF (nl .NE. 1) GO TO 104
565 ifac(1) = n
566 ifac(2) = nf
567 tpi = 6.28318530717959
568 argh = tpi/float(n)
569 is = 0
570 nfm1 = nf-1
571 l1 = 1
572 IF (nfm1 .EQ. 0) RETURN
573!OCL NOVREC
574 DO 110 k1=1,nfm1
575 ip = int(ifac(k1+2))
576 ld = 0
577 l2 = l1*ip
578 ido = n/l2
579 ipm = ip-1
580 DO 109 j=1,ipm
581 ld = ld+l1
582 i = is
583 argld = float(ld)*argh
584 fi = 0
585!OCL SCALAR
586 DO 108 ii=3,ido,2
587 i = i+2
588 fi = fi+1
589 arg = fi*argld
590 wa(i-1) = cos(arg)
591 wa(i) = sin(arg)
592 108 CONTINUE
593 is = is+ido
594 109 CONTINUE
595 l1 = l2
596 110 CONTINUE
597 RETURN
598 END
599
600C> @brief Perform the backward radix-2 FFT step
601C>
602C> This subroutine computes one stage of the backward Fast Fourier Transform
603C> (FFT) using a radix-2 decomposition. It is used as part of the larger
604C> FFT process and applies twiddle factors to transform the input data.
605C>
606C> @param ido The leading dimension, representing the number of data points in a transform section
607C> @param l1 The number of radix-2 stages processed in this step
608C> @param cc Input complex array containing intermediate FFT data
609C> @param ch Output complex array storing transformed data
610C> @param wa1 Twiddle factors array for weighting the FFT computation
611C>
612C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
613 SUBROUTINE radb2 (IDO,L1,CC,CH,WA1)
614 dimension cc(ido,2,l1) ,ch(ido,l1,2) ,
615 1 wa1(*)
616 DO 101 k=1,l1
617 ch(1,k,1) = cc(1,1,k)+cc(ido,2,k)
618 ch(1,k,2) = cc(1,1,k)-cc(ido,2,k)
619 101 CONTINUE
620 IF (ido.LT.2) THEN
621 GO TO 107
622 ELSEIF (ido.EQ.2) THEN
623 GO TO 105
624 ELSE
625 GO TO 102
626 ENDIF
627 102 idp2 = ido+2
628!OCL NOVREC
629 DO 104 k=1,l1
630 DO 103 i=3,ido,2
631 ic = idp2-i
632 ch(i-1,k,1) = cc(i-1,1,k)+cc(ic-1,2,k)
633 tr2 = cc(i-1,1,k)-cc(ic-1,2,k)
634 ch(i,k,1) = cc(i,1,k)-cc(ic,2,k)
635 ti2 = cc(i,1,k)+cc(ic,2,k)
636 ch(i-1,k,2) = wa1(i-2)*tr2-wa1(i-1)*ti2
637 ch(i,k,2) = wa1(i-2)*ti2+wa1(i-1)*tr2
638 103 CONTINUE
639 104 CONTINUE
640 IF (mod(ido,2) .EQ. 1) RETURN
641 105 DO 106 k=1,l1
642 ch(ido,k,1) = cc(ido,1,k)+cc(ido,1,k)
643 ch(ido,k,2) = -(cc(1,2,k)+cc(1,2,k))
644 106 CONTINUE
645 107 RETURN
646 END
647
648C> @brief Perform the backward radix-3 FFT step
649C>
650C> This subroutine computes one stage of the backward Fast Fourier Transform
651C> (FFT) using a radix-3 decomposition. It applies twiddle factors and combines
652C> intermediate results to reconstruct transformed data.
653C>
654C> @param ido The leading dimension, representing the number of data points in a transform section
655C> @param l1 The number of radix-3 stages processed in this step
656C> @param cc Input complex array containing intermediate FFT data
657C> @param ch Output complex array storing transformed data
658C> @param wa1 Twiddle factors array for the first rotation
659C> @param wa2 Twiddle factors array for the second rotation
660C>
661C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
662 SUBROUTINE radb3 (IDO,L1,CC,CH,WA1,WA2)
663 dimension cc(ido,3,l1) ,ch(ido,l1,3) ,
664 1 wa1(*) ,wa2(*)
665 DATA taur,taui /-.5,.866025403784439/
666 DO 101 k=1,l1
667 tr2 = cc(ido,2,k)+cc(ido,2,k)
668 cr2 = cc(1,1,k)+taur*tr2
669 ch(1,k,1) = cc(1,1,k)+tr2
670 ci3 = taui*(cc(1,3,k)+cc(1,3,k))
671 ch(1,k,2) = cr2-ci3
672 ch(1,k,3) = cr2+ci3
673 101 CONTINUE
674 IF (ido .EQ. 1) RETURN
675 idp2 = ido+2
676!OCL NOVREC
677 DO 103 k=1,l1
678 DO 102 i=3,ido,2
679 ic = idp2-i
680 tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
681 cr2 = cc(i-1,1,k)+taur*tr2
682 ch(i-1,k,1) = cc(i-1,1,k)+tr2
683 ti2 = cc(i,3,k)-cc(ic,2,k)
684 ci2 = cc(i,1,k)+taur*ti2
685 ch(i,k,1) = cc(i,1,k)+ti2
686 cr3 = taui*(cc(i-1,3,k)-cc(ic-1,2,k))
687 ci3 = taui*(cc(i,3,k)+cc(ic,2,k))
688 dr2 = cr2-ci3
689 dr3 = cr2+ci3
690 di2 = ci2+cr3
691 di3 = ci2-cr3
692 ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
693 ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
694 ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
695 ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
696 102 CONTINUE
697 103 CONTINUE
698 RETURN
699 END
700
701C> @brief Perform the backward radix-4 FFT step
702C>
703C> This subroutine computes one stage of the backward Fast Fourier Transform
704C> (FFT) using a radix-4 decomposition. It applies the FFT formulae by processing
705C> the input data and performing necessary rotations and twiddle factor multiplications.
706C> The result is stored in the output array `ch`.
707C>
708C> @param ido The leading dimension, representing the number of data points in a transform section
709C> @param l1 The number of radix-4 stages processed in this step
710C> @param cc Input complex array containing intermediate FFT data
711C> @param ch Output complex array storing transformed data
712C> @param wa1 Twiddle factors array for the first rotation
713C> @param wa2 Twiddle factors array for the second rotation
714C> @param wa3 Twiddle factors array for the third rotation
715C>
716C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
717 SUBROUTINE radb4 (IDO,L1,CC,CH,WA1,WA2,WA3)
718 dimension cc(ido,4,l1) ,ch(ido,l1,4) ,
719 1 wa1(*) ,wa2(*) ,wa3(*)
720 DATA sqrt2 /1.414213562373095/
721 DO 101 k=1,l1
722 tr1 = cc(1,1,k)-cc(ido,4,k)
723 tr2 = cc(1,1,k)+cc(ido,4,k)
724 tr3 = cc(ido,2,k)+cc(ido,2,k)
725 tr4 = cc(1,3,k)+cc(1,3,k)
726 ch(1,k,1) = tr2+tr3
727 ch(1,k,2) = tr1-tr4
728 ch(1,k,3) = tr2-tr3
729 ch(1,k,4) = tr1+tr4
730 101 CONTINUE
731 IF (ido.LT.2) THEN
732 GO TO 107
733 ELSEIF (ido.EQ.2) THEN
734 GO TO 105
735 ELSE
736 GO TO 102
737 ENDIF
738 102 idp2 = ido+2
739!OCL NOVREC
740 DO 104 k=1,l1
741 DO 103 i=3,ido,2
742 ic = idp2-i
743 ti1 = cc(i,1,k)+cc(ic,4,k)
744 ti2 = cc(i,1,k)-cc(ic,4,k)
745 ti3 = cc(i,3,k)-cc(ic,2,k)
746 tr4 = cc(i,3,k)+cc(ic,2,k)
747 tr1 = cc(i-1,1,k)-cc(ic-1,4,k)
748 tr2 = cc(i-1,1,k)+cc(ic-1,4,k)
749 ti4 = cc(i-1,3,k)-cc(ic-1,2,k)
750 tr3 = cc(i-1,3,k)+cc(ic-1,2,k)
751 ch(i-1,k,1) = tr2+tr3
752 cr3 = tr2-tr3
753 ch(i,k,1) = ti2+ti3
754 ci3 = ti2-ti3
755 cr2 = tr1-tr4
756 cr4 = tr1+tr4
757 ci2 = ti1+ti4
758 ci4 = ti1-ti4
759 ch(i-1,k,2) = wa1(i-2)*cr2-wa1(i-1)*ci2
760 ch(i,k,2) = wa1(i-2)*ci2+wa1(i-1)*cr2
761 ch(i-1,k,3) = wa2(i-2)*cr3-wa2(i-1)*ci3
762 ch(i,k,3) = wa2(i-2)*ci3+wa2(i-1)*cr3
763 ch(i-1,k,4) = wa3(i-2)*cr4-wa3(i-1)*ci4
764 ch(i,k,4) = wa3(i-2)*ci4+wa3(i-1)*cr4
765 103 CONTINUE
766 104 CONTINUE
767 IF (mod(ido,2) .EQ. 1) RETURN
768 105 CONTINUE
769 DO 106 k=1,l1
770 ti1 = cc(1,2,k)+cc(1,4,k)
771 ti2 = cc(1,4,k)-cc(1,2,k)
772 tr1 = cc(ido,1,k)-cc(ido,3,k)
773 tr2 = cc(ido,1,k)+cc(ido,3,k)
774 ch(ido,k,1) = tr2+tr2
775 ch(ido,k,2) = sqrt2*(tr1-ti1)
776 ch(ido,k,3) = ti2+ti2
777 ch(ido,k,4) = -sqrt2*(tr1+ti1)
778 106 CONTINUE
779 107 RETURN
780 END
781
782C> @brief Perform the backward radix-5 FFT step
783C>
784C> This subroutine computes one stage of the backward Fast Fourier Transform
785C> (FFT) using a radix-5 decomposition. It processes input data in sections,
786C> applies necessary rotations and twiddle factors, and stores the transformed
787C> results in the output array `ch`. The operation involves calculating
788C> intermediate results for each data point and applying appropriate scaling and
789C> twiddle factor multiplication.
790C>
791C> @param ido The leading dimension, representing the number of data points in a transform section
792C> @param l1 The number of radix-5 stages processed in this step
793C> @param cc Input complex array containing intermediate FFT data
794C> @param ch Output complex array storing transformed data
795C> @param wa1 Twiddle factors array for the first rotation
796C> @param wa2 Twiddle factors array for the second rotation
797C> @param wa3 Twiddle factors array for the third rotation
798C> @param wa4 Twiddle factors array for the fourth rotation
799C>
800C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
801 SUBROUTINE radb5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
802 dimension cc(ido,5,l1) ,ch(ido,l1,5) ,
803 1 wa1(*) ,wa2(*) ,wa3(*) ,wa4(*)
804 DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
805 1-.809016994374947,.587785252292473/
806 DO 101 k=1,l1
807 ti5 = cc(1,3,k)+cc(1,3,k)
808 ti4 = cc(1,5,k)+cc(1,5,k)
809 tr2 = cc(ido,2,k)+cc(ido,2,k)
810 tr3 = cc(ido,4,k)+cc(ido,4,k)
811 ch(1,k,1) = cc(1,1,k)+tr2+tr3
812 cr2 = cc(1,1,k)+tr11*tr2+tr12*tr3
813 cr3 = cc(1,1,k)+tr12*tr2+tr11*tr3
814 ci5 = ti11*ti5+ti12*ti4
815 ci4 = ti12*ti5-ti11*ti4
816 ch(1,k,2) = cr2-ci5
817 ch(1,k,3) = cr3-ci4
818 ch(1,k,4) = cr3+ci4
819 ch(1,k,5) = cr2+ci5
820 101 CONTINUE
821 IF (ido .EQ. 1) RETURN
822 idp2 = ido+2
823 DO 103 k=1,l1
824 DO 102 i=3,ido,2
825 ic = idp2-i
826 ti5 = cc(i,3,k)+cc(ic,2,k)
827 ti2 = cc(i,3,k)-cc(ic,2,k)
828 ti4 = cc(i,5,k)+cc(ic,4,k)
829 ti3 = cc(i,5,k)-cc(ic,4,k)
830 tr5 = cc(i-1,3,k)-cc(ic-1,2,k)
831 tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
832 tr4 = cc(i-1,5,k)-cc(ic-1,4,k)
833 tr3 = cc(i-1,5,k)+cc(ic-1,4,k)
834 ch(i-1,k,1) = cc(i-1,1,k)+tr2+tr3
835 ch(i,k,1) = cc(i,1,k)+ti2+ti3
836 cr2 = cc(i-1,1,k)+tr11*tr2+tr12*tr3
837 ci2 = cc(i,1,k)+tr11*ti2+tr12*ti3
838 cr3 = cc(i-1,1,k)+tr12*tr2+tr11*tr3
839 ci3 = cc(i,1,k)+tr12*ti2+tr11*ti3
840 cr5 = ti11*tr5+ti12*tr4
841 ci5 = ti11*ti5+ti12*ti4
842 cr4 = ti12*tr5-ti11*tr4
843 ci4 = ti12*ti5-ti11*ti4
844 dr3 = cr3-ci4
845 dr4 = cr3+ci4
846 di3 = ci3+cr4
847 di4 = ci3-cr4
848 dr5 = cr2+ci5
849 dr2 = cr2-ci5
850 di5 = ci2-cr5
851 di2 = ci2+cr5
852 ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
853 ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
854 ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
855 ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
856 ch(i-1,k,4) = wa3(i-2)*dr4-wa3(i-1)*di4
857 ch(i,k,4) = wa3(i-2)*di4+wa3(i-1)*dr4
858 ch(i-1,k,5) = wa4(i-2)*dr5-wa4(i-1)*di5
859 ch(i,k,5) = wa4(i-2)*di5+wa4(i-1)*dr5
860 102 CONTINUE
861 103 CONTINUE
862 RETURN
863 END
864
865C> @brief Computes the backward FFT stage using a generalized FFT algorithm.
866C>
867C> This subroutine performs a backward Fast Fourier Transform (FFT) step with
868C> a generalized radix-`ip` approach. It computes the necessary intermediate
869C> steps, applying the inverse FFT transformation to the complex input data
870C> stored in `cc`, and stores the result in the output arrays `ch` and `c2`.
871C> Twiddle factors, stored in `wa`, are used for scaling and rotations at
872C> each stage of the computation.
873C>
874C> @param ido The leading dimension for data points in a transform section
875C> @param ip The number of subgroups in the FFT (radix)
876C> @param l1 The number of radix stages to process
877C> @param idl1 The leading dimension for another input array
878C> @param cc Input complex array containing intermediate FFT data
879C> @param c1 Array containing the first set of intermediate FFT results
880C> @param c2 Array to store additional intermediate FFT results
881C> @param ch Output complex array to store the final transformed data
882C> @param ch2 Output complex array storing intermediate results for the
883C> second stage of transformation
884C> @param wa Array containing the twiddle factors used in the FFT
885C>
886C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
887 SUBROUTINE radbg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
888 dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
889 1 c1(ido,l1,ip) ,c2(idl1,ip),
890 2 ch2(idl1,ip) ,wa(*)
891 DATA tpi/6.28318530717959/
892 arg = tpi/float(ip)
893 dcp = cos(arg)
894 dsp = sin(arg)
895 idp2 = ido+2
896 nbd = (ido-1)/2
897 ipp2 = ip+2
898 ipph = (ip+1)/2
899 IF (ido .LT. l1) GO TO 103
900 DO 102 k=1,l1
901 DO 101 i=1,ido
902 ch(i,k,1) = cc(i,1,k)
903 101 CONTINUE
904 102 CONTINUE
905 GO TO 106
906 103 DO 105 i=1,ido
907 DO 104 k=1,l1
908 ch(i,k,1) = cc(i,1,k)
909 104 CONTINUE
910 105 CONTINUE
911!OCL NOVREC
912 106 DO 108 j=2,ipph
913 jc = ipp2-j
914 j2 = j+j
915 DO 107 k=1,l1
916 ch(1,k,j) = cc(ido,j2-2,k)+cc(ido,j2-2,k)
917 ch(1,k,jc) = cc(1,j2-1,k)+cc(1,j2-1,k)
918 107 CONTINUE
919 108 CONTINUE
920 IF (ido .EQ. 1) GO TO 116
921 IF (nbd .LT. l1) GO TO 112
922!OCL NOVREC
923 DO 111 j=2,ipph
924 jc = ipp2-j
925 DO 110 k=1,l1
926 DO 109 i=3,ido,2
927 ic = idp2-i
928 ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
929 ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
930 ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
931 ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
932 109 CONTINUE
933 110 CONTINUE
934 111 CONTINUE
935 GO TO 116
936 112 DO 115 j=2,ipph
937 jc = ipp2-j
938 DO 114 i=3,ido,2
939 ic = idp2-i
940 DO 113 k=1,l1
941 ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
942 ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
943 ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
944 ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
945 113 CONTINUE
946 114 CONTINUE
947 115 CONTINUE
948 116 ar1 = 1.
949 ai1 = 0.
950!OCL NOVREC
951 DO 120 l=2,ipph
952 lc = ipp2-l
953 ar1h = dcp*ar1-dsp*ai1
954 ai1 = dcp*ai1+dsp*ar1
955 ar1 = ar1h
956 DO 117 ik=1,idl1
957 c2(ik,l) = ch2(ik,1)+ar1*ch2(ik,2)
958 c2(ik,lc) = ai1*ch2(ik,ip)
959 117 CONTINUE
960 dc2 = ar1
961 ds2 = ai1
962 ar2 = ar1
963 ai2 = ai1
964!OCL NOVREC
965 DO 119 j=3,ipph
966 jc = ipp2-j
967 ar2h = dc2*ar2-ds2*ai2
968 ai2 = dc2*ai2+ds2*ar2
969 ar2 = ar2h
970 DO 118 ik=1,idl1
971 c2(ik,l) = c2(ik,l)+ar2*ch2(ik,j)
972 c2(ik,lc) = c2(ik,lc)+ai2*ch2(ik,jc)
973 118 CONTINUE
974 119 CONTINUE
975 120 CONTINUE
976!OCL NOVREC
977 DO 122 j=2,ipph
978 DO 121 ik=1,idl1
979 ch2(ik,1) = ch2(ik,1)+ch2(ik,j)
980 121 CONTINUE
981 122 CONTINUE
982!OCL NOVREC
983 DO 124 j=2,ipph
984 jc = ipp2-j
985 DO 123 k=1,l1
986 ch(1,k,j) = c1(1,k,j)-c1(1,k,jc)
987 ch(1,k,jc) = c1(1,k,j)+c1(1,k,jc)
988 123 CONTINUE
989 124 CONTINUE
990 IF (ido .EQ. 1) GO TO 132
991 IF (nbd .LT. l1) GO TO 128
992!OCL NOVREC
993 DO 127 j=2,ipph
994 jc = ipp2-j
995 DO 126 k=1,l1
996 DO 125 i=3,ido,2
997 ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
998 ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
999 ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
1000 ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
1001 125 CONTINUE
1002 126 CONTINUE
1003 127 CONTINUE
1004 GO TO 132
1005 128 DO 131 j=2,ipph
1006 jc = ipp2-j
1007 DO 130 i=3,ido,2
1008 DO 129 k=1,l1
1009 ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
1010 ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
1011 ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
1012 ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
1013 129 CONTINUE
1014 130 CONTINUE
1015 131 CONTINUE
1016 132 CONTINUE
1017 IF (ido .EQ. 1) RETURN
1018 DO 133 ik=1,idl1
1019 c2(ik,1) = ch2(ik,1)
1020 133 CONTINUE
1021 DO 135 j=2,ip
1022 DO 134 k=1,l1
1023 c1(1,k,j) = ch(1,k,j)
1024 134 CONTINUE
1025 135 CONTINUE
1026 IF (nbd .GT. l1) GO TO 139
1027 is = -ido
1028 DO 138 j=2,ip
1029 is = is+ido
1030 idij = is
1031 DO 137 i=3,ido,2
1032 idij = idij+2
1033 DO 136 k=1,l1
1034 c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
1035 c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
1036 136 CONTINUE
1037 137 CONTINUE
1038 138 CONTINUE
1039 GO TO 143
1040 139 is = -ido
1041!OCL NOVREC
1042 DO 142 j=2,ip
1043 is = is+ido
1044 DO 141 k=1,l1
1045 idij = is
1046 DO 140 i=3,ido,2
1047 idij = idij+2
1048 c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
1049 c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
1050 140 CONTINUE
1051 141 CONTINUE
1052 142 CONTINUE
1053 143 RETURN
1054 END
1055
1056C> @brief Perform a radix-2 FFT step for even-odd decomposition of data
1057C>
1058C> This subroutine computes a single stage of the radix-2 Fast Fourier Transform
1059C> (FFT) by separating the even and odd indexed values from the input data
1060C> array `cc`, and stores the results in the output array `ch`. It processes
1061C> intermediate data with the twiddle factors from `wa1` for scaling and rotation.
1062C>
1063C> @param ido The size of the FFT sub-segment (number of data points to process)
1064C> @param l1 The number of stages to process in the FFT
1065C> @param cc Input complex array of size `ido x l1 x 2`, containing the data
1066C> for the FFT transformation
1067C> @param ch Output complex array of size `ido x 2 x l1`, where the
1068C> transformed data will be stored
1069C> @param wa1 Array of twiddle factors used to scale and rotate the data
1070C> during the FFT calculation
1071C>
1072C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1073 SUBROUTINE radf2 (IDO,L1,CC,CH,WA1)
1074 dimension ch(ido,2,l1) ,cc(ido,l1,2) ,
1075 1 wa1(*)
1076 DO 101 k=1,l1
1077 ch(1,1,k) = cc(1,k,1)+cc(1,k,2)
1078 ch(ido,2,k) = cc(1,k,1)-cc(1,k,2)
1079 101 CONTINUE
1080 IF (ido.LT.2) THEN
1081 GO TO 107
1082 ELSEIF (ido.EQ.2) THEN
1083 GO TO 105
1084 ELSE
1085 GO TO 102
1086 ENDIF
1087 102 idp2 = ido+2
1088 DO 104 k=1,l1
1089 DO 103 i=3,ido,2
1090 ic = idp2-i
1091 tr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1092 ti2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1093 ch(i,1,k) = cc(i,k,1)+ti2
1094 ch(ic,2,k) = ti2-cc(i,k,1)
1095 ch(i-1,1,k) = cc(i-1,k,1)+tr2
1096 ch(ic-1,2,k) = cc(i-1,k,1)-tr2
1097 103 CONTINUE
1098 104 CONTINUE
1099 IF (mod(ido,2) .EQ. 1) RETURN
1100 105 DO 106 k=1,l1
1101 ch(1,2,k) = -cc(ido,k,2)
1102 ch(ido,1,k) = cc(ido,k,1)
1103 106 CONTINUE
1104 107 RETURN
1105 END
1106
1107C> @brief Performs a radix-3 FFT step for decomposing the data into three components.
1108C>
1109C> This subroutine computes a single stage of the radix-3 Fast Fourier Transform
1110C> (FFT) by separating the input data `cc` into three components (one for each
1111C> of the three terms in the FFT). The output data is stored in the array `ch`,
1112C> with twiddle factors applied from the arrays `wa1` and `wa2`. This subroutine
1113C> operates for various sizes of FFT sub-segments, controlled by the parameter
1114C> `ido`.
1115C>
1116C> The subroutine handles different configurations based on the value of `ido`,
1117C> which affects how the data is processed.
1118C>
1119C> @param ido The size of the FFT sub-segment (number of data points to process)
1120C> @param l1 The number of stages to process in the FFT
1121C> @param cc Input complex array of size `ido x l1 x 3`, containing the data
1122C> for the FFT transformation
1123C> @param ch Output complex array of size `ido x 3 x l1`, where the transformed
1124C> data will be stored
1125C> @param wa1 Array of twiddle factors for the first part of the FFT calculation
1126C> @param wa2 Array of twiddle factors for the second part of the FFT calculation
1127C>
1128C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1129 SUBROUTINE radf3 (IDO,L1,CC,CH,WA1,WA2)
1130 dimension ch(ido,3,l1) ,cc(ido,l1,3) ,
1131 1 wa1(*) ,wa2(*)
1132 DATA taur,taui /-.5,.866025403784439/
1133 DO 101 k=1,l1
1134 cr2 = cc(1,k,2)+cc(1,k,3)
1135 ch(1,1,k) = cc(1,k,1)+cr2
1136 ch(1,3,k) = taui*(cc(1,k,3)-cc(1,k,2))
1137 ch(ido,2,k) = cc(1,k,1)+taur*cr2
1138 101 CONTINUE
1139 IF (ido .EQ. 1) RETURN
1140 idp2 = ido+2
1141 DO 103 k=1,l1
1142 DO 102 i=3,ido,2
1143 ic = idp2-i
1144 dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1145 di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1146 dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1147 di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1148 cr2 = dr2+dr3
1149 ci2 = di2+di3
1150 ch(i-1,1,k) = cc(i-1,k,1)+cr2
1151 ch(i,1,k) = cc(i,k,1)+ci2
1152 tr2 = cc(i-1,k,1)+taur*cr2
1153 ti2 = cc(i,k,1)+taur*ci2
1154 tr3 = taui*(di2-di3)
1155 ti3 = taui*(dr3-dr2)
1156 ch(i-1,3,k) = tr2+tr3
1157 ch(ic-1,2,k) = tr2-tr3
1158 ch(i,3,k) = ti2+ti3
1159 ch(ic,2,k) = ti3-ti2
1160 102 CONTINUE
1161 103 CONTINUE
1162 RETURN
1163 END
1164
1165C> @brief Perform a radix-4 FFT step for decomposing the data into four components
1166C>
1167C> This subroutine computes a single stage of the radix-4 Fast Fourier Transform
1168C> (FFT) by separating the input data `cc` into four components (one for each
1169C> of the four terms in the FFT). The output data is stored in the array `ch`,
1170C> with twiddle factors applied from the arrays `wa1`, `wa2`, and `wa3`. This
1171C> subroutine operates for various sizes of FFT sub-segments, controlled by the
1172C> parameter `ido`.
1173C>
1174C> @param ido The size of the FFT sub-segment (number of data points to process)
1175C> @param l1 The number of stages to process in the FFT
1176C> @param cc Input complex array of size `ido x l1 x 4`, containing the data
1177C> for the FFT transformation
1178C> @param ch Output complex array of size `ido x 4 x l1`, where the transformed
1179C> data will be stored
1180C> @param wa1 Array of twiddle factors for the first part of the FFT calculation
1181C> @param wa2 Array of twiddle factors for the second part of the FFT calculation
1182C> @param wa3 Array of twiddle factors for the third part of the FFT calculation
1183C>
1184C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1185 SUBROUTINE radf4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1186 dimension cc(ido,l1,4) ,ch(ido,4,l1) ,
1187 1 wa1(*) ,wa2(*) ,wa3(*)
1188 DATA hsqt2 /.7071067811865475/
1189 DO 101 k=1,l1
1190 tr1 = cc(1,k,2)+cc(1,k,4)
1191 tr2 = cc(1,k,1)+cc(1,k,3)
1192 ch(1,1,k) = tr1+tr2
1193 ch(ido,4,k) = tr2-tr1
1194 ch(ido,2,k) = cc(1,k,1)-cc(1,k,3)
1195 ch(1,3,k) = cc(1,k,4)-cc(1,k,2)
1196 101 CONTINUE
1197 IF (ido.LT.2) THEN
1198 GO TO 107
1199 ELSEIF (ido.EQ.2) THEN
1200 GO TO 105
1201 ELSE
1202 GO TO 102
1203 ENDIF
1204 102 idp2 = ido+2
1205!OCL NOVREC
1206 DO 104 k=1,l1
1207 DO 103 i=3,ido,2
1208 ic = idp2-i
1209 cr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1210 ci2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1211 cr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1212 ci3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1213 cr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
1214 ci4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
1215 tr1 = cr2+cr4
1216 tr4 = cr4-cr2
1217 ti1 = ci2+ci4
1218 ti4 = ci2-ci4
1219 ti2 = cc(i,k,1)+ci3
1220 ti3 = cc(i,k,1)-ci3
1221 tr2 = cc(i-1,k,1)+cr3
1222 tr3 = cc(i-1,k,1)-cr3
1223 ch(i-1,1,k) = tr1+tr2
1224 ch(ic-1,4,k) = tr2-tr1
1225 ch(i,1,k) = ti1+ti2
1226 ch(ic,4,k) = ti1-ti2
1227 ch(i-1,3,k) = ti4+tr3
1228 ch(ic-1,2,k) = tr3-ti4
1229 ch(i,3,k) = tr4+ti3
1230 ch(ic,2,k) = tr4-ti3
1231 103 CONTINUE
1232 104 CONTINUE
1233 IF (mod(ido,2) .EQ. 1) RETURN
1234 105 CONTINUE
1235 DO 106 k=1,l1
1236 ti1 = -hsqt2*(cc(ido,k,2)+cc(ido,k,4))
1237 tr1 = hsqt2*(cc(ido,k,2)-cc(ido,k,4))
1238 ch(ido,1,k) = tr1+cc(ido,k,1)
1239 ch(ido,3,k) = cc(ido,k,1)-tr1
1240 ch(1,2,k) = ti1-cc(ido,k,3)
1241 ch(1,4,k) = ti1+cc(ido,k,3)
1242 106 CONTINUE
1243 107 RETURN
1244 END
1245
1246C> @brief Perform a radix-5 FFT step for decomposing the data into five components
1247C>
1248C> This subroutine computes a single stage of the radix-5 Fast Fourier Transform
1249C> (FFT) by separating the input data `cc` into five components (one for each
1250C> of the five terms in the FFT). The output data is stored in the array `ch`,
1251C> with twiddle factors applied from the arrays `wa1`, `wa2`, `wa3`, and `wa4`.
1252C> This subroutine operates for various sizes of FFT sub-segments, controlled by
1253C> the parameter `ido`.
1254C>
1255C> @param ido The size of the FFT sub-segment (number of data points to process)
1256C> @param l1 The number of stages to process in the FFT
1257C> @param cc Input complex array of size `ido x l1 x 5`, containing the data
1258C> for the FFT transformation
1259C> @param ch Output complex array of size `ido x 5 x l1`, where the transformed
1260C> data will be stored
1261C> @param wa1 Array of twiddle factors for the first part of the FFT calculation
1262C> @param wa2 Array of twiddle factors for the second part of the FFT calculation
1263C> @param wa3 Array of twiddle factors for the third part of the FFT calculation
1264C> @param wa4 Array of twiddle factors for the fourth part of the FFT calculation
1265C>
1266C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1267 SUBROUTINE radf5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1268 dimension cc(ido,l1,5) ,ch(ido,5,l1) ,
1269 1 wa1(*) ,wa2(*) ,wa3(*) ,wa4(*)
1270 DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
1271 1-.809016994374947,.587785252292473/
1272 DO 101 k=1,l1
1273 cr2 = cc(1,k,5)+cc(1,k,2)
1274 ci5 = cc(1,k,5)-cc(1,k,2)
1275 cr3 = cc(1,k,4)+cc(1,k,3)
1276 ci4 = cc(1,k,4)-cc(1,k,3)
1277 ch(1,1,k) = cc(1,k,1)+cr2+cr3
1278 ch(ido,2,k) = cc(1,k,1)+tr11*cr2+tr12*cr3
1279 ch(1,3,k) = ti11*ci5+ti12*ci4
1280 ch(ido,4,k) = cc(1,k,1)+tr12*cr2+tr11*cr3
1281 ch(1,5,k) = ti12*ci5-ti11*ci4
1282 101 CONTINUE
1283 IF (ido .EQ. 1) RETURN
1284 idp2 = ido+2
1285 DO 103 k=1,l1
1286 DO 102 i=3,ido,2
1287 ic = idp2-i
1288 dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1289 di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1290 dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1291 di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1292 dr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
1293 di4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
1294 dr5 = wa4(i-2)*cc(i-1,k,5)+wa4(i-1)*cc(i,k,5)
1295 di5 = wa4(i-2)*cc(i,k,5)-wa4(i-1)*cc(i-1,k,5)
1296 cr2 = dr2+dr5
1297 ci5 = dr5-dr2
1298 cr5 = di2-di5
1299 ci2 = di2+di5
1300 cr3 = dr3+dr4
1301 ci4 = dr4-dr3
1302 cr4 = di3-di4
1303 ci3 = di3+di4
1304 ch(i-1,1,k) = cc(i-1,k,1)+cr2+cr3
1305 ch(i,1,k) = cc(i,k,1)+ci2+ci3
1306 tr2 = cc(i-1,k,1)+tr11*cr2+tr12*cr3
1307 ti2 = cc(i,k,1)+tr11*ci2+tr12*ci3
1308 tr3 = cc(i-1,k,1)+tr12*cr2+tr11*cr3
1309 ti3 = cc(i,k,1)+tr12*ci2+tr11*ci3
1310 tr5 = ti11*cr5+ti12*cr4
1311 ti5 = ti11*ci5+ti12*ci4
1312 tr4 = ti12*cr5-ti11*cr4
1313 ti4 = ti12*ci5-ti11*ci4
1314 ch(i-1,3,k) = tr2+tr5
1315 ch(ic-1,2,k) = tr2-tr5
1316 ch(i,3,k) = ti2+ti5
1317 ch(ic,2,k) = ti5-ti2
1318 ch(i-1,5,k) = tr3+tr4
1319 ch(ic-1,4,k) = tr3-tr4
1320 ch(i,5,k) = ti3+ti4
1321 ch(ic,4,k) = ti4-ti3
1322 102 CONTINUE
1323 103 CONTINUE
1324 RETURN
1325 END
1326
1327C> @brief Compute a general radix-based FFT step
1328C>
1329C> This subroutine implements a mixed-radix Fast Fourier Transform (FFT) step,
1330C> decomposing the input data for an arbitrary radix `ip`. It reorganizes the
1331C> input data `cc` and applies twiddle factors from `wa` to generate the output
1332C> in `ch` and `ch2`.
1333C>
1334C> The transformation is influenced by `ido`, which controls the number of
1335C> data points per FFT sub-segment, and `l1`, which determines the number of
1336C> processing stages.
1337C>
1338C> @param ido The number of data points per sub-segment
1339C> @param ip The radix of the FFT step
1340C> @param l1 The number of stages to process
1341C> @param idl1 The transformed data length parameter, dependent on `l1`
1342C> @param cc Input complex array of size `ido x ip x l1`, containing the data to transform
1343C> @param c1 Temporary working array for intermediate calculations
1344C> @param c2 Temporary working array for intermediate calculations
1345C> @param ch Output complex array of size `ido x l1 x ip`, where transformed data is stored
1346C> @param ch2 Another working output array used for intermediate transformations
1347C> @param wa Twiddle factor array used to apply phase shifts in the FFT calculation
1348C>
1349C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1350 SUBROUTINE radfg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1351 dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
1352 1 c1(ido,l1,ip) ,c2(idl1,ip),
1353 2 ch2(idl1,ip) ,wa(*)
1354 DATA tpi/6.28318530717959/
1355 arg = tpi/float(ip)
1356 dcp = cos(arg)
1357 dsp = sin(arg)
1358 ipph = (ip+1)/2
1359 ipp2 = ip+2
1360 idp2 = ido+2
1361 nbd = (ido-1)/2
1362 IF (ido .EQ. 1) GO TO 119
1363 DO 101 ik=1,idl1
1364 ch2(ik,1) = c2(ik,1)
1365 101 CONTINUE
1366 DO 103 j=2,ip
1367 DO 102 k=1,l1
1368 ch(1,k,j) = c1(1,k,j)
1369 102 CONTINUE
1370 103 CONTINUE
1371 IF (nbd .GT. l1) GO TO 107
1372 is = -ido
1373 DO 106 j=2,ip
1374 is = is+ido
1375 idij = is
1376 DO 105 i=3,ido,2
1377 idij = idij+2
1378 DO 104 k=1,l1
1379 ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
1380 ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
1381 104 CONTINUE
1382 105 CONTINUE
1383 106 CONTINUE
1384 GO TO 111
1385 107 is = -ido
1386 DO 110 j=2,ip
1387 is = is+ido
1388 DO 109 k=1,l1
1389 idij = is
1390 DO 108 i=3,ido,2
1391 idij = idij+2
1392 ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
1393 ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
1394 108 CONTINUE
1395 109 CONTINUE
1396 110 CONTINUE
1397 111 IF (nbd .LT. l1) GO TO 115
1398 DO 114 j=2,ipph
1399 jc = ipp2-j
1400 DO 113 k=1,l1
1401 DO 112 i=3,ido,2
1402 c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
1403 c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
1404 c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
1405 c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
1406 112 CONTINUE
1407 113 CONTINUE
1408 114 CONTINUE
1409 GO TO 121
1410 115 DO 118 j=2,ipph
1411 jc = ipp2-j
1412 DO 117 i=3,ido,2
1413 DO 116 k=1,l1
1414 c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
1415 c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
1416 c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
1417 c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
1418 116 CONTINUE
1419 117 CONTINUE
1420 118 CONTINUE
1421 GO TO 121
1422 119 DO 120 ik=1,idl1
1423 c2(ik,1) = ch2(ik,1)
1424 120 CONTINUE
1425 121 DO 123 j=2,ipph
1426 jc = ipp2-j
1427 DO 122 k=1,l1
1428 c1(1,k,j) = ch(1,k,j)+ch(1,k,jc)
1429 c1(1,k,jc) = ch(1,k,jc)-ch(1,k,j)
1430 122 CONTINUE
1431 123 CONTINUE
1432C
1433 ar1 = 1.
1434 ai1 = 0.
1435 DO 127 l=2,ipph
1436 lc = ipp2-l
1437 ar1h = dcp*ar1-dsp*ai1
1438 ai1 = dcp*ai1+dsp*ar1
1439 ar1 = ar1h
1440 DO 124 ik=1,idl1
1441 ch2(ik,l) = c2(ik,1)+ar1*c2(ik,2)
1442 ch2(ik,lc) = ai1*c2(ik,ip)
1443 124 CONTINUE
1444 dc2 = ar1
1445 ds2 = ai1
1446 ar2 = ar1
1447 ai2 = ai1
1448 DO 126 j=3,ipph
1449 jc = ipp2-j
1450 ar2h = dc2*ar2-ds2*ai2
1451 ai2 = dc2*ai2+ds2*ar2
1452 ar2 = ar2h
1453 DO 125 ik=1,idl1
1454 ch2(ik,l) = ch2(ik,l)+ar2*c2(ik,j)
1455 ch2(ik,lc) = ch2(ik,lc)+ai2*c2(ik,jc)
1456 125 CONTINUE
1457 126 CONTINUE
1458 127 CONTINUE
1459 DO 129 j=2,ipph
1460 DO 128 ik=1,idl1
1461 ch2(ik,1) = ch2(ik,1)+c2(ik,j)
1462 128 CONTINUE
1463 129 CONTINUE
1464C
1465 IF (ido .LT. l1) GO TO 132
1466 DO 131 k=1,l1
1467 DO 130 i=1,ido
1468 cc(i,1,k) = ch(i,k,1)
1469 130 CONTINUE
1470 131 CONTINUE
1471 GO TO 135
1472 132 DO 134 i=1,ido
1473 DO 133 k=1,l1
1474 cc(i,1,k) = ch(i,k,1)
1475 133 CONTINUE
1476 134 CONTINUE
1477 135 DO 137 j=2,ipph
1478 jc = ipp2-j
1479 j2 = j+j
1480 DO 136 k=1,l1
1481 cc(ido,j2-2,k) = ch(1,k,j)
1482 cc(1,j2-1,k) = ch(1,k,jc)
1483 136 CONTINUE
1484 137 CONTINUE
1485 IF (ido .EQ. 1) RETURN
1486 IF (nbd .LT. l1) GO TO 141
1487 DO 140 j=2,ipph
1488 jc = ipp2-j
1489 j2 = j+j
1490 DO 139 k=1,l1
1491 DO 138 i=3,ido,2
1492 ic = idp2-i
1493 cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1494 cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1495 cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1496 cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1497 138 CONTINUE
1498 139 CONTINUE
1499 140 CONTINUE
1500 RETURN
1501 141 DO 144 j=2,ipph
1502 jc = ipp2-j
1503 j2 = j+j
1504 DO 143 i=3,ido,2
1505 ic = idp2-i
1506 DO 142 k=1,l1
1507 cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1508 cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1509 cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1510 cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1511 142 CONTINUE
1512 143 CONTINUE
1513 144 CONTINUE
1514 RETURN
1515 END
radf5
subroutine radf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4)
Perform a radix-5 FFT step for decomposing the data into five components.
Definition fftpack.F:1268
radbg
subroutine radbg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa)
Computes the backward FFT stage using a generalized FFT algorithm.
Definition fftpack.F:888
rffti1
subroutine rffti1(n, wa, ifac)
Compute the initialization factors for the real-input FFT.
Definition fftpack.F:533
scfft
subroutine scfft(isign, n, scale, x, y, table, work, isys)
Compute a real-to-complex discrete Fourier transform.
Definition fftpack.F:294
rfftb
subroutine rfftb(n, r, wsave)
Compute the real periodic sequence from its Fourier coefs.
Definition fftpack.F:352
rfftb1
subroutine rfftb1(n, c, ch, wa, ifac)
Compute the inverse fast Fourier transform (IFFT) using a mixed-radix algorithm.
Definition fftpack.F:387
drcft
subroutine drcft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
Compute a set of m complex discrete n-point Fourier transforms of real data.
Definition fftpack.F:189
rfftf1
subroutine rfftf1(n, c, ch, wa, ifac)
Compute the real forward fast Fourier transform (FFT) using a mixed-radix algorithm.
Definition fftpack.F:461
dcrft
subroutine dcrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
Computes a set of m real discrete n-point Fourier transforms of complex conjugate even data.
Definition fftpack.F:44
radf4
subroutine radf4(ido, l1, cc, ch, wa1, wa2, wa3)
Perform a radix-4 FFT step for decomposing the data into four components.
Definition fftpack.F:1186
rfftf
subroutine rfftf(n, r, wsave)
Compute a forward transform of a real periodic sequence.
Definition fftpack.F:335
radf3
subroutine radf3(ido, l1, cc, ch, wa1, wa2)
Performs a radix-3 FFT step for decomposing the data into three components.
Definition fftpack.F:1130
radfg
subroutine radfg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa)
Compute a general radix-based FFT step.
Definition fftpack.F:1351
radb5
subroutine radb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4)
Perform the backward radix-5 FFT step.
Definition fftpack.F:802
radf2
subroutine radf2(ido, l1, cc, ch, wa1)
Perform a radix-2 FFT step for even-odd decomposition of data.
Definition fftpack.F:1074
radb3
subroutine radb3(ido, l1, cc, ch, wa1, wa2)
Perform the backward radix-3 FFT step.
Definition fftpack.F:663
radb2
subroutine radb2(ido, l1, cc, ch, wa1)
Perform the backward radix-2 FFT step.
Definition fftpack.F:614
scrft
subroutine scrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
Compute a set of m real discrete n-point Fourier transforms of complex conjugate even data.
Definition fftpack.F:94
srcft
subroutine srcft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
Compute a set of m complex discrete n-point Fourier transforms of real data.
Definition fftpack.F:248
csfft
subroutine csfft(isign, n, scale, x, y, table, work, isys)
Compute a complex discrete Fourier transform of real data.
Definition fftpack.F:138
rffti
subroutine rffti(n, wsave)
Initialize the array wsave which is used in both rfftf() and rfftb().
Definition fftpack.F:365
radb4
subroutine radb4(ido, l1, cc, ch, wa1, wa2, wa3)
Perform the backward radix-4 FFT step.
Definition fftpack.F:718