RavEngine
Loading...
Searching...
No Matches
fftpack.h
1/*
2 Interface for the f2c translation of fftpack as found on http://www.netlib.org/fftpack/
3
4 FFTPACK license:
5
6 http://www.cisl.ucar.edu/css/software/fftpack5/ftpk.html
7
8 Copyright (c) 2004 the University Corporation for Atmospheric
9 Research ("UCAR"). All rights reserved. Developed by NCAR's
10 Computational and Information Systems Laboratory, UCAR,
11 www.cisl.ucar.edu.
12
13 Redistribution and use of the Software in source and binary forms,
14 with or without modification, is permitted provided that the
15 following conditions are met:
16
17 - Neither the names of NCAR's Computational and Information Systems
18 Laboratory, the University Corporation for Atmospheric Research,
19 nor the names of its sponsors or contributors may be used to
20 endorse or promote products derived from this Software without
21 specific prior written permission.
22
23 - Redistributions of source code must retain the above copyright
24 notices, this list of conditions, and the disclaimer below.
25
26 - Redistributions in binary form must reproduce the above copyright
27 notice, this list of conditions, and the disclaimer below in the
28 documentation and/or other materials provided with the
29 distribution.
30
31 THIS SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
32 EXPRESS OR IMPLIED, INCLUDING, BUT NOT LIMITED TO THE WARRANTIES OF
33 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
34 NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR COPYRIGHT
35 HOLDERS BE LIABLE FOR ANY CLAIM, INDIRECT, INCIDENTAL, SPECIAL,
36 EXEMPLARY, OR CONSEQUENTIAL DAMAGES OR OTHER LIABILITY, WHETHER IN AN
37 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
38 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS WITH THE
39 SOFTWARE.
40
41 ChangeLog:
42 2011/10/02: this is my first release of this file.
43*/
44
45#ifndef FFTPACK_H
46#define FFTPACK_H
47
48#ifdef __cplusplus
49extern "C" {
50#endif
51
52// just define FFTPACK_DOUBLE_PRECISION if you want to build it as a double precision fft
53
54#ifndef FFTPACK_DOUBLE_PRECISION
55 typedef float fftpack_real;
56 typedef int fftpack_int;
57#else
58 typedef double fftpack_real;
59 typedef int fftpack_int;
60#endif
61
62 void cffti(fftpack_int n, fftpack_real *wsave);
63
64 void cfftf(fftpack_int n, fftpack_real *c, fftpack_real *wsave);
65
66 void cfftb(fftpack_int n, fftpack_real *c, fftpack_real *wsave);
67
68 void rffti(fftpack_int n, fftpack_real *wsave);
69 void rfftf(fftpack_int n, fftpack_real *r, fftpack_real *wsave);
70 void rfftb(fftpack_int n, fftpack_real *r, fftpack_real *wsave);
71
72 void cosqi(fftpack_int n, fftpack_real *wsave);
73 void cosqf(fftpack_int n, fftpack_real *x, fftpack_real *wsave);
74 void cosqb(fftpack_int n, fftpack_real *x, fftpack_real *wsave);
75
76 void costi(fftpack_int n, fftpack_real *wsave);
77 void cost(fftpack_int n, fftpack_real *x, fftpack_real *wsave);
78
79 void sinqi(fftpack_int n, fftpack_real *wsave);
80 void sinqb(fftpack_int n, fftpack_real *x, fftpack_real *wsave);
81 void sinqf(fftpack_int n, fftpack_real *x, fftpack_real *wsave);
82
83 void sinti(fftpack_int n, fftpack_real *wsave);
84 void sint(fftpack_int n, fftpack_real *x, fftpack_real *wsave);
85
86#ifdef __cplusplus
87}
88#endif
89
90#endif /* FFTPACK_H */
91
92/*
93
94 FFTPACK
95
96* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
97
98 version 4 april 1985
99
100 a package of fortran subprograms for the fast fourier
101 transform of periodic and other symmetric sequences
102
103 by
104
105 paul n swarztrauber
106
107 national center for atmospheric research boulder,colorado 80307
108
109 which is sponsored by the national science foundation
110
111* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
112
113
114this package consists of programs which perform fast fourier
115transforms for both complex and real periodic sequences and
116certain other symmetric sequences that are listed below.
117
1181. rffti initialize rfftf and rfftb
1192. rfftf forward transform of a real periodic sequence
1203. rfftb backward transform of a real coefficient array
121
1224. ezffti initialize ezfftf and ezfftb
1235. ezfftf a simplified real periodic forward transform
1246. ezfftb a simplified real periodic backward transform
125
1267. sinti initialize sint
1278. sint sine transform of a real odd sequence
128
1299. costi initialize cost
13010. cost cosine transform of a real even sequence
131
13211. sinqi initialize sinqf and sinqb
13312. sinqf forward sine transform with odd wave numbers
13413. sinqb unnormalized inverse of sinqf
135
13614. cosqi initialize cosqf and cosqb
13715. cosqf forward cosine transform with odd wave numbers
13816. cosqb unnormalized inverse of cosqf
139
14017. cffti initialize cfftf and cfftb
14118. cfftf forward transform of a complex periodic sequence
14219. cfftb unnormalized inverse of cfftf
143
144
145******************************************************************
146
147subroutine rffti(n,wsave)
148
149 ****************************************************************
150
151subroutine rffti initializes the array wsave which is used in
152both rfftf and rfftb. the prime factorization of n together with
153a tabulation of the trigonometric functions are computed and
154stored in wsave.
155
156input parameter
157
158n the length of the sequence to be transformed.
159
160output parameter
161
162wsave a work array which must be dimensioned at least 2*n+15.
163 the same work array can be used for both rfftf and rfftb
164 as long as n remains unchanged. different wsave arrays
165 are required for different values of n. the contents of
166 wsave must not be changed between calls of rfftf or rfftb.
167
168******************************************************************
169
170subroutine rfftf(n,r,wsave)
171
172******************************************************************
173
174subroutine rfftf computes the fourier coefficients of a real
175perodic sequence (fourier analysis). the transform is defined
176below at output parameter r.
177
178input parameters
179
180n the length of the array r to be transformed. the method
181 is most efficient when n is a product of small primes.
182 n may change so long as different work arrays are provided
183
184r a real array of length n which contains the sequence
185 to be transformed
186
187wsave a work array which must be dimensioned at least 2*n+15.
188 in the program that calls rfftf. the wsave array must be
189 initialized by calling subroutine rffti(n,wsave) and a
190 different wsave array must be used for each different
191 value of n. this initialization does not have to be
192 repeated so long as n remains unchanged thus subsequent
193 transforms can be obtained faster than the first.
194 the same wsave array can be used by rfftf and rfftb.
195
196
197output parameters
198
199r r(1) = the sum from i=1 to i=n of r(i)
200
201 if n is even set l =n/2 , if n is odd set l = (n+1)/2
202
203 then for k = 2,...,l
204
205 r(2*k-2) = the sum from i = 1 to i = n of
206
207 r(i)*cos((k-1)*(i-1)*2*pi/n)
208
209 r(2*k-1) = the sum from i = 1 to i = n of
210
211 -r(i)*sin((k-1)*(i-1)*2*pi/n)
212
213 if n is even
214
215 r(n) = the sum from i = 1 to i = n of
216
217 (-1)**(i-1)*r(i)
218
219 ***** note
220 this transform is unnormalized since a call of rfftf
221 followed by a call of rfftb will multiply the input
222 sequence by n.
223
224wsave contains results which must not be destroyed between
225 calls of rfftf or rfftb.
226
227
228******************************************************************
229
230subroutine rfftb(n,r,wsave)
231
232******************************************************************
233
234subroutine rfftb computes the real perodic sequence from its
235fourier coefficients (fourier synthesis). the transform is defined
236below at output parameter r.
237
238input parameters
239
240n the length of the array r to be transformed. the method
241 is most efficient when n is a product of small primes.
242 n may change so long as different work arrays are provided
243
244r a real array of length n which contains the sequence
245 to be transformed
246
247wsave a work array which must be dimensioned at least 2*n+15.
248 in the program that calls rfftb. the wsave array must be
249 initialized by calling subroutine rffti(n,wsave) and a
250 different wsave array must be used for each different
251 value of n. this initialization does not have to be
252 repeated so long as n remains unchanged thus subsequent
253 transforms can be obtained faster than the first.
254 the same wsave array can be used by rfftf and rfftb.
255
256
257output parameters
258
259r for n even and for i = 1,...,n
260
261 r(i) = r(1)+(-1)**(i-1)*r(n)
262
263 plus the sum from k=2 to k=n/2 of
264
265 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)
266
267 -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
268
269 for n odd and for i = 1,...,n
270
271 r(i) = r(1) plus the sum from k=2 to k=(n+1)/2 of
272
273 2.*r(2*k-2)*cos((k-1)*(i-1)*2*pi/n)
274
275 -2.*r(2*k-1)*sin((k-1)*(i-1)*2*pi/n)
276
277 ***** note
278 this transform is unnormalized since a call of rfftf
279 followed by a call of rfftb will multiply the input
280 sequence by n.
281
282wsave contains results which must not be destroyed between
283 calls of rfftb or rfftf.
284
285******************************************************************
286
287subroutine sinti(n,wsave)
288
289******************************************************************
290
291subroutine sinti initializes the array wsave which is used in
292subroutine sint. the prime factorization of n together with
293a tabulation of the trigonometric functions are computed and
294stored in wsave.
295
296input parameter
297
298n the length of the sequence to be transformed. the method
299 is most efficient when n+1 is a product of small primes.
300
301output parameter
302
303wsave a work array with at least int(2.5*n+15) locations.
304 different wsave arrays are required for different values
305 of n. the contents of wsave must not be changed between
306 calls of sint.
307
308******************************************************************
309
310subroutine sint(n,x,wsave)
311
312******************************************************************
313
314subroutine sint computes the discrete fourier sine transform
315of an odd sequence x(i). the transform is defined below at
316output parameter x.
317
318sint is the unnormalized inverse of itself since a call of sint
319followed by another call of sint will multiply the input sequence
320x by 2*(n+1).
321
322the array wsave which is used by subroutine sint must be
323initialized by calling subroutine sinti(n,wsave).
324
325input parameters
326
327n the length of the sequence to be transformed. the method
328 is most efficient when n+1 is the product of small primes.
329
330x an array which contains the sequence to be transformed
331
332
333wsave a work array with dimension at least int(2.5*n+15)
334 in the program that calls sint. the wsave array must be
335 initialized by calling subroutine sinti(n,wsave) and a
336 different wsave array must be used for each different
337 value of n. this initialization does not have to be
338 repeated so long as n remains unchanged thus subsequent
339 transforms can be obtained faster than the first.
340
341output parameters
342
343x for i=1,...,n
344
345 x(i)= the sum from k=1 to k=n
346
347 2*x(k)*sin(k*i*pi/(n+1))
348
349 a call of sint followed by another call of
350 sint will multiply the sequence x by 2*(n+1).
351 hence sint is the unnormalized inverse
352 of itself.
353
354wsave contains initialization calculations which must not be
355 destroyed between calls of sint.
356
357******************************************************************
358
359subroutine costi(n,wsave)
360
361******************************************************************
362
363subroutine costi initializes the array wsave which is used in
364subroutine cost. the prime factorization of n together with
365a tabulation of the trigonometric functions are computed and
366stored in wsave.
367
368input parameter
369
370n the length of the sequence to be transformed. the method
371 is most efficient when n-1 is a product of small primes.
372
373output parameter
374
375wsave a work array which must be dimensioned at least 3*n+15.
376 different wsave arrays are required for different values
377 of n. the contents of wsave must not be changed between
378 calls of cost.
379
380******************************************************************
381
382subroutine cost(n,x,wsave)
383
384******************************************************************
385
386subroutine cost computes the discrete fourier cosine transform
387of an even sequence x(i). the transform is defined below at output
388parameter x.
389
390cost is the unnormalized inverse of itself since a call of cost
391followed by another call of cost will multiply the input sequence
392x by 2*(n-1). the transform is defined below at output parameter x
393
394the array wsave which is used by subroutine cost must be
395initialized by calling subroutine costi(n,wsave).
396
397input parameters
398
399n the length of the sequence x. n must be greater than 1.
400 the method is most efficient when n-1 is a product of
401 small primes.
402
403x an array which contains the sequence to be transformed
404
405wsave a work array which must be dimensioned at least 3*n+15
406 in the program that calls cost. the wsave array must be
407 initialized by calling subroutine costi(n,wsave) and a
408 different wsave array must be used for each different
409 value of n. this initialization does not have to be
410 repeated so long as n remains unchanged thus subsequent
411 transforms can be obtained faster than the first.
412
413output parameters
414
415x for i=1,...,n
416
417 x(i) = x(1)+(-1)**(i-1)*x(n)
418
419 + the sum from k=2 to k=n-1
420
421 2*x(k)*cos((k-1)*(i-1)*pi/(n-1))
422
423 a call of cost followed by another call of
424 cost will multiply the sequence x by 2*(n-1)
425 hence cost is the unnormalized inverse
426 of itself.
427
428wsave contains initialization calculations which must not be
429 destroyed between calls of cost.
430
431******************************************************************
432
433subroutine sinqi(n,wsave)
434
435******************************************************************
436
437subroutine sinqi initializes the array wsave which is used in
438both sinqf and sinqb. the prime factorization of n together with
439a tabulation of the trigonometric functions are computed and
440stored in wsave.
441
442input parameter
443
444n the length of the sequence to be transformed. the method
445 is most efficient when n is a product of small primes.
446
447output parameter
448
449wsave a work array which must be dimensioned at least 3*n+15.
450 the same work array can be used for both sinqf and sinqb
451 as long as n remains unchanged. different wsave arrays
452 are required for different values of n. the contents of
453 wsave must not be changed between calls of sinqf or sinqb.
454
455******************************************************************
456
457subroutine sinqf(n,x,wsave)
458
459******************************************************************
460
461subroutine sinqf computes the fast fourier transform of quarter
462wave data. that is , sinqf computes the coefficients in a sine
463series representation with only odd wave numbers. the transform
464is defined below at output parameter x.
465
466sinqb is the unnormalized inverse of sinqf since a call of sinqf
467followed by a call of sinqb will multiply the input sequence x
468by 4*n.
469
470the array wsave which is used by subroutine sinqf must be
471initialized by calling subroutine sinqi(n,wsave).
472
473
474input parameters
475
476n the length of the array x to be transformed. the method
477 is most efficient when n is a product of small primes.
478
479x an array which contains the sequence to be transformed
480
481wsave a work array which must be dimensioned at least 3*n+15.
482 in the program that calls sinqf. the wsave array must be
483 initialized by calling subroutine sinqi(n,wsave) and a
484 different wsave array must be used for each different
485 value of n. this initialization does not have to be
486 repeated so long as n remains unchanged thus subsequent
487 transforms can be obtained faster than the first.
488
489output parameters
490
491x for i=1,...,n
492
493 x(i) = (-1)**(i-1)*x(n)
494
495 + the sum from k=1 to k=n-1 of
496
497 2*x(k)*sin((2*i-1)*k*pi/(2*n))
498
499 a call of sinqf followed by a call of
500 sinqb will multiply the sequence x by 4*n.
501 therefore sinqb is the unnormalized inverse
502 of sinqf.
503
504wsave contains initialization calculations which must not
505 be destroyed between calls of sinqf or sinqb.
506
507******************************************************************
508
509subroutine sinqb(n,x,wsave)
510
511******************************************************************
512
513subroutine sinqb computes the fast fourier transform of quarter
514wave data. that is , sinqb computes a sequence from its
515representation in terms of a sine series with odd wave numbers.
516the transform is defined below at output parameter x.
517
518sinqf is the unnormalized inverse of sinqb since a call of sinqb
519followed by a call of sinqf will multiply the input sequence x
520by 4*n.
521
522the array wsave which is used by subroutine sinqb must be
523initialized by calling subroutine sinqi(n,wsave).
524
525
526input parameters
527
528n the length of the array x to be transformed. the method
529 is most efficient when n is a product of small primes.
530
531x an array which contains the sequence to be transformed
532
533wsave a work array which must be dimensioned at least 3*n+15.
534 in the program that calls sinqb. the wsave array must be
535 initialized by calling subroutine sinqi(n,wsave) and a
536 different wsave array must be used for each different
537 value of n. this initialization does not have to be
538 repeated so long as n remains unchanged thus subsequent
539 transforms can be obtained faster than the first.
540
541output parameters
542
543x for i=1,...,n
544
545 x(i)= the sum from k=1 to k=n of
546
547 4*x(k)*sin((2k-1)*i*pi/(2*n))
548
549 a call of sinqb followed by a call of
550 sinqf will multiply the sequence x by 4*n.
551 therefore sinqf is the unnormalized inverse
552 of sinqb.
553
554wsave contains initialization calculations which must not
555 be destroyed between calls of sinqb or sinqf.
556
557******************************************************************
558
559subroutine cosqi(n,wsave)
560
561******************************************************************
562
563subroutine cosqi initializes the array wsave which is used in
564both cosqf and cosqb. the prime factorization of n together with
565a tabulation of the trigonometric functions are computed and
566stored in wsave.
567
568input parameter
569
570n the length of the array to be transformed. the method
571 is most efficient when n is a product of small primes.
572
573output parameter
574
575wsave a work array which must be dimensioned at least 3*n+15.
576 the same work array can be used for both cosqf and cosqb
577 as long as n remains unchanged. different wsave arrays
578 are required for different values of n. the contents of
579 wsave must not be changed between calls of cosqf or cosqb.
580
581******************************************************************
582
583subroutine cosqf(n,x,wsave)
584
585******************************************************************
586
587subroutine cosqf computes the fast fourier transform of quarter
588wave data. that is , cosqf computes the coefficients in a cosine
589series representation with only odd wave numbers. the transform
590is defined below at output parameter x
591
592cosqf is the unnormalized inverse of cosqb since a call of cosqf
593followed by a call of cosqb will multiply the input sequence x
594by 4*n.
595
596the array wsave which is used by subroutine cosqf must be
597initialized by calling subroutine cosqi(n,wsave).
598
599
600input parameters
601
602n the length of the array x to be transformed. the method
603 is most efficient when n is a product of small primes.
604
605x an array which contains the sequence to be transformed
606
607wsave a work array which must be dimensioned at least 3*n+15
608 in the program that calls cosqf. the wsave array must be
609 initialized by calling subroutine cosqi(n,wsave) and a
610 different wsave array must be used for each different
611 value of n. this initialization does not have to be
612 repeated so long as n remains unchanged thus subsequent
613 transforms can be obtained faster than the first.
614
615output parameters
616
617x for i=1,...,n
618
619 x(i) = x(1) plus the sum from k=2 to k=n of
620
621 2*x(k)*cos((2*i-1)*(k-1)*pi/(2*n))
622
623 a call of cosqf followed by a call of
624 cosqb will multiply the sequence x by 4*n.
625 therefore cosqb is the unnormalized inverse
626 of cosqf.
627
628wsave contains initialization calculations which must not
629 be destroyed between calls of cosqf or cosqb.
630
631******************************************************************
632
633subroutine cosqb(n,x,wsave)
634
635******************************************************************
636
637subroutine cosqb computes the fast fourier transform of quarter
638wave data. that is , cosqb computes a sequence from its
639representation in terms of a cosine series with odd wave numbers.
640the transform is defined below at output parameter x.
641
642cosqb is the unnormalized inverse of cosqf since a call of cosqb
643followed by a call of cosqf will multiply the input sequence x
644by 4*n.
645
646the array wsave which is used by subroutine cosqb must be
647initialized by calling subroutine cosqi(n,wsave).
648
649
650input parameters
651
652n the length of the array x to be transformed. the method
653 is most efficient when n is a product of small primes.
654
655x an array which contains the sequence to be transformed
656
657wsave a work array that must be dimensioned at least 3*n+15
658 in the program that calls cosqb. the wsave array must be
659 initialized by calling subroutine cosqi(n,wsave) and a
660 different wsave array must be used for each different
661 value of n. this initialization does not have to be
662 repeated so long as n remains unchanged thus subsequent
663 transforms can be obtained faster than the first.
664
665output parameters
666
667x for i=1,...,n
668
669 x(i)= the sum from k=1 to k=n of
670
671 4*x(k)*cos((2*k-1)*(i-1)*pi/(2*n))
672
673 a call of cosqb followed by a call of
674 cosqf will multiply the sequence x by 4*n.
675 therefore cosqf is the unnormalized inverse
676 of cosqb.
677
678wsave contains initialization calculations which must not
679 be destroyed between calls of cosqb or cosqf.
680
681******************************************************************
682
683subroutine cffti(n,wsave)
684
685******************************************************************
686
687subroutine cffti initializes the array wsave which is used in
688both cfftf and cfftb. the prime factorization of n together with
689a tabulation of the trigonometric functions are computed and
690stored in wsave.
691
692input parameter
693
694n the length of the sequence to be transformed
695
696output parameter
697
698wsave a work array which must be dimensioned at least 4*n+15
699 the same work array can be used for both cfftf and cfftb
700 as long as n remains unchanged. different wsave arrays
701 are required for different values of n. the contents of
702 wsave must not be changed between calls of cfftf or cfftb.
703
704******************************************************************
705
706subroutine cfftf(n,c,wsave)
707
708******************************************************************
709
710subroutine cfftf computes the forward complex discrete fourier
711transform (the fourier analysis). equivalently , cfftf computes
712the fourier coefficients of a complex periodic sequence.
713the transform is defined below at output parameter c.
714
715the transform is not normalized. to obtain a normalized transform
716the output must be divided by n. otherwise a call of cfftf
717followed by a call of cfftb will multiply the sequence by n.
718
719the array wsave which is used by subroutine cfftf must be
720initialized by calling subroutine cffti(n,wsave).
721
722input parameters
723
724
725n the length of the complex sequence c. the method is
726 more efficient when n is the product of small primes. n
727
728c a complex array of length n which contains the sequence
729
730wsave a real work array which must be dimensioned at least 4n+15
731 in the program that calls cfftf. the wsave array must be
732 initialized by calling subroutine cffti(n,wsave) and a
733 different wsave array must be used for each different
734 value of n. this initialization does not have to be
735 repeated so long as n remains unchanged thus subsequent
736 transforms can be obtained faster than the first.
737 the same wsave array can be used by cfftf and cfftb.
738
739output parameters
740
741c for j=1,...,n
742
743 c(j)=the sum from k=1,...,n of
744
745 c(k)*exp(-i*(j-1)*(k-1)*2*pi/n)
746
747 where i=sqrt(-1)
748
749wsave contains initialization calculations which must not be
750 destroyed between calls of subroutine cfftf or cfftb
751
752******************************************************************
753
754subroutine cfftb(n,c,wsave)
755
756******************************************************************
757
758subroutine cfftb computes the backward complex discrete fourier
759transform (the fourier synthesis). equivalently , cfftb computes
760a complex periodic sequence from its fourier coefficients.
761the transform is defined below at output parameter c.
762
763a call of cfftf followed by a call of cfftb will multiply the
764sequence by n.
765
766the array wsave which is used by subroutine cfftb must be
767initialized by calling subroutine cffti(n,wsave).
768
769input parameters
770
771
772n the length of the complex sequence c. the method is
773 more efficient when n is the product of small primes.
774
775c a complex array of length n which contains the sequence
776
777wsave a real work array which must be dimensioned at least 4n+15
778 in the program that calls cfftb. the wsave array must be
779 initialized by calling subroutine cffti(n,wsave) and a
780 different wsave array must be used for each different
781 value of n. this initialization does not have to be
782 repeated so long as n remains unchanged thus subsequent
783 transforms can be obtained faster than the first.
784 the same wsave array can be used by cfftf and cfftb.
785
786output parameters
787
788c for j=1,...,n
789
790 c(j)=the sum from k=1,...,n of
791
792 c(k)*exp(i*(j-1)*(k-1)*2*pi/n)
793
794 where i=sqrt(-1)
795
796wsave contains initialization calculations which must not be
797 destroyed between calls of subroutine cfftf or cfftb
798
799*/
glm::sint
signed int sint
Definition integer.hpp:55