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 49 extern "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 114 this package consists of programs which perform fast fourier 115 transforms for both complex and real periodic sequences and 116 certain other symmetric sequences that are listed below. 117 118 1. rffti initialize rfftf and rfftb 119 2. rfftf forward transform of a real periodic sequence 120 3. rfftb backward transform of a real coefficient array 121 122 4. ezffti initialize ezfftf and ezfftb 123 5. ezfftf a simplified real periodic forward transform 124 6. ezfftb a simplified real periodic backward transform 125 126 7. sinti initialize sint 127 8. sint sine transform of a real odd sequence 128 129 9. costi initialize cost 130 10. cost cosine transform of a real even sequence 131 132 11. sinqi initialize sinqf and sinqb 133 12. sinqf forward sine transform with odd wave numbers 134 13. sinqb unnormalized inverse of sinqf 135 136 14. cosqi initialize cosqf and cosqb 137 15. cosqf forward cosine transform with odd wave numbers 138 16. cosqb unnormalized inverse of cosqf 139 140 17. cffti initialize cfftf and cfftb 141 18. cfftf forward transform of a complex periodic sequence 142 19. cfftb unnormalized inverse of cfftf 143 144 145 ****************************************************************** 146 147 subroutine rffti(n,wsave) 148 149 **************************************************************** 150 151 subroutine rffti initializes the array wsave which is used in 152 both rfftf and rfftb. the prime factorization of n together with 153 a tabulation of the trigonometric functions are computed and 154 stored in wsave. 155 156 input parameter 157 158 n the length of the sequence to be transformed. 159 160 output parameter 161 162 wsave 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 170 subroutine rfftf(n,r,wsave) 171 172 ****************************************************************** 173 174 subroutine rfftf computes the fourier coefficients of a real 175 perodic sequence (fourier analysis). the transform is defined 176 below at output parameter r. 177 178 input parameters 179 180 n 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 184 r a real array of length n which contains the sequence 185 to be transformed 186 187 wsave 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 197 output parameters 198 199 r 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 224 wsave contains results which must not be destroyed between 225 calls of rfftf or rfftb. 226 227 228 ****************************************************************** 229 230 subroutine rfftb(n,r,wsave) 231 232 ****************************************************************** 233 234 subroutine rfftb computes the real perodic sequence from its 235 fourier coefficients (fourier synthesis). the transform is defined 236 below at output parameter r. 237 238 input parameters 239 240 n 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 244 r a real array of length n which contains the sequence 245 to be transformed 246 247 wsave 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 257 output parameters 258 259 r 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 282 wsave contains results which must not be destroyed between 283 calls of rfftb or rfftf. 284 285 ****************************************************************** 286 287 subroutine sinti(n,wsave) 288 289 ****************************************************************** 290 291 subroutine sinti initializes the array wsave which is used in 292 subroutine sint. the prime factorization of n together with 293 a tabulation of the trigonometric functions are computed and 294 stored in wsave. 295 296 input parameter 297 298 n the length of the sequence to be transformed. the method 299 is most efficient when n+1 is a product of small primes. 300 301 output parameter 302 303 wsave 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 310 subroutine sint(n,x,wsave) 311 312 ****************************************************************** 313 314 subroutine sint computes the discrete fourier sine transform 315 of an odd sequence x(i). the transform is defined below at 316 output parameter x. 317 318 sint is the unnormalized inverse of itself since a call of sint 319 followed by another call of sint will multiply the input sequence 320 x by 2*(n+1). 321 322 the array wsave which is used by subroutine sint must be 323 initialized by calling subroutine sinti(n,wsave). 324 325 input parameters 326 327 n the length of the sequence to be transformed. the method 328 is most efficient when n+1 is the product of small primes. 329 330 x an array which contains the sequence to be transformed 331 332 333 wsave 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 341 output parameters 342 343 x 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 354 wsave contains initialization calculations which must not be 355 destroyed between calls of sint. 356 357 ****************************************************************** 358 359 subroutine costi(n,wsave) 360 361 ****************************************************************** 362 363 subroutine costi initializes the array wsave which is used in 364 subroutine cost. the prime factorization of n together with 365 a tabulation of the trigonometric functions are computed and 366 stored in wsave. 367 368 input parameter 369 370 n the length of the sequence to be transformed. the method 371 is most efficient when n-1 is a product of small primes. 372 373 output parameter 374 375 wsave 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 382 subroutine cost(n,x,wsave) 383 384 ****************************************************************** 385 386 subroutine cost computes the discrete fourier cosine transform 387 of an even sequence x(i). the transform is defined below at output 388 parameter x. 389 390 cost is the unnormalized inverse of itself since a call of cost 391 followed by another call of cost will multiply the input sequence 392 x by 2*(n-1). the transform is defined below at output parameter x 393 394 the array wsave which is used by subroutine cost must be 395 initialized by calling subroutine costi(n,wsave). 396 397 input parameters 398 399 n 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 403 x an array which contains the sequence to be transformed 404 405 wsave 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 413 output parameters 414 415 x 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 428 wsave contains initialization calculations which must not be 429 destroyed between calls of cost. 430 431 ****************************************************************** 432 433 subroutine sinqi(n,wsave) 434 435 ****************************************************************** 436 437 subroutine sinqi initializes the array wsave which is used in 438 both sinqf and sinqb. the prime factorization of n together with 439 a tabulation of the trigonometric functions are computed and 440 stored in wsave. 441 442 input parameter 443 444 n the length of the sequence to be transformed. the method 445 is most efficient when n is a product of small primes. 446 447 output parameter 448 449 wsave 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 457 subroutine sinqf(n,x,wsave) 458 459 ****************************************************************** 460 461 subroutine sinqf computes the fast fourier transform of quarter 462 wave data. that is , sinqf computes the coefficients in a sine 463 series representation with only odd wave numbers. the transform 464 is defined below at output parameter x. 465 466 sinqb is the unnormalized inverse of sinqf since a call of sinqf 467 followed by a call of sinqb will multiply the input sequence x 468 by 4*n. 469 470 the array wsave which is used by subroutine sinqf must be 471 initialized by calling subroutine sinqi(n,wsave). 472 473 474 input parameters 475 476 n the length of the array x to be transformed. the method 477 is most efficient when n is a product of small primes. 478 479 x an array which contains the sequence to be transformed 480 481 wsave 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 489 output parameters 490 491 x 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 504 wsave contains initialization calculations which must not 505 be destroyed between calls of sinqf or sinqb. 506 507 ****************************************************************** 508 509 subroutine sinqb(n,x,wsave) 510 511 ****************************************************************** 512 513 subroutine sinqb computes the fast fourier transform of quarter 514 wave data. that is , sinqb computes a sequence from its 515 representation in terms of a sine series with odd wave numbers. 516 the transform is defined below at output parameter x. 517 518 sinqf is the unnormalized inverse of sinqb since a call of sinqb 519 followed by a call of sinqf will multiply the input sequence x 520 by 4*n. 521 522 the array wsave which is used by subroutine sinqb must be 523 initialized by calling subroutine sinqi(n,wsave). 524 525 526 input parameters 527 528 n the length of the array x to be transformed. the method 529 is most efficient when n is a product of small primes. 530 531 x an array which contains the sequence to be transformed 532 533 wsave 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 541 output parameters 542 543 x 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 554 wsave contains initialization calculations which must not 555 be destroyed between calls of sinqb or sinqf. 556 557 ****************************************************************** 558 559 subroutine cosqi(n,wsave) 560 561 ****************************************************************** 562 563 subroutine cosqi initializes the array wsave which is used in 564 both cosqf and cosqb. the prime factorization of n together with 565 a tabulation of the trigonometric functions are computed and 566 stored in wsave. 567 568 input parameter 569 570 n the length of the array to be transformed. the method 571 is most efficient when n is a product of small primes. 572 573 output parameter 574 575 wsave 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 583 subroutine cosqf(n,x,wsave) 584 585 ****************************************************************** 586 587 subroutine cosqf computes the fast fourier transform of quarter 588 wave data. that is , cosqf computes the coefficients in a cosine 589 series representation with only odd wave numbers. the transform 590 is defined below at output parameter x 591 592 cosqf is the unnormalized inverse of cosqb since a call of cosqf 593 followed by a call of cosqb will multiply the input sequence x 594 by 4*n. 595 596 the array wsave which is used by subroutine cosqf must be 597 initialized by calling subroutine cosqi(n,wsave). 598 599 600 input parameters 601 602 n the length of the array x to be transformed. the method 603 is most efficient when n is a product of small primes. 604 605 x an array which contains the sequence to be transformed 606 607 wsave 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 615 output parameters 616 617 x 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 628 wsave contains initialization calculations which must not 629 be destroyed between calls of cosqf or cosqb. 630 631 ****************************************************************** 632 633 subroutine cosqb(n,x,wsave) 634 635 ****************************************************************** 636 637 subroutine cosqb computes the fast fourier transform of quarter 638 wave data. that is , cosqb computes a sequence from its 639 representation in terms of a cosine series with odd wave numbers. 640 the transform is defined below at output parameter x. 641 642 cosqb is the unnormalized inverse of cosqf since a call of cosqb 643 followed by a call of cosqf will multiply the input sequence x 644 by 4*n. 645 646 the array wsave which is used by subroutine cosqb must be 647 initialized by calling subroutine cosqi(n,wsave). 648 649 650 input parameters 651 652 n the length of the array x to be transformed. the method 653 is most efficient when n is a product of small primes. 654 655 x an array which contains the sequence to be transformed 656 657 wsave 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 665 output parameters 666 667 x 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 678 wsave contains initialization calculations which must not 679 be destroyed between calls of cosqb or cosqf. 680 681 ****************************************************************** 682 683 subroutine cffti(n,wsave) 684 685 ****************************************************************** 686 687 subroutine cffti initializes the array wsave which is used in 688 both cfftf and cfftb. the prime factorization of n together with 689 a tabulation of the trigonometric functions are computed and 690 stored in wsave. 691 692 input parameter 693 694 n the length of the sequence to be transformed 695 696 output parameter 697 698 wsave 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 706 subroutine cfftf(n,c,wsave) 707 708 ****************************************************************** 709 710 subroutine cfftf computes the forward complex discrete fourier 711 transform (the fourier analysis). equivalently , cfftf computes 712 the fourier coefficients of a complex periodic sequence. 713 the transform is defined below at output parameter c. 714 715 the transform is not normalized. to obtain a normalized transform 716 the output must be divided by n. otherwise a call of cfftf 717 followed by a call of cfftb will multiply the sequence by n. 718 719 the array wsave which is used by subroutine cfftf must be 720 initialized by calling subroutine cffti(n,wsave). 721 722 input parameters 723 724 725 n the length of the complex sequence c. the method is 726 more efficient when n is the product of small primes. n 727 728 c a complex array of length n which contains the sequence 729 730 wsave 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 739 output parameters 740 741 c 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 749 wsave contains initialization calculations which must not be 750 destroyed between calls of subroutine cfftf or cfftb 751 752 ****************************************************************** 753 754 subroutine cfftb(n,c,wsave) 755 756 ****************************************************************** 757 758 subroutine cfftb computes the backward complex discrete fourier 759 transform (the fourier synthesis). equivalently , cfftb computes 760 a complex periodic sequence from its fourier coefficients. 761 the transform is defined below at output parameter c. 762 763 a call of cfftf followed by a call of cfftb will multiply the 764 sequence by n. 765 766 the array wsave which is used by subroutine cfftb must be 767 initialized by calling subroutine cffti(n,wsave). 768 769 input parameters 770 771 772 n the length of the complex sequence c. the method is 773 more efficient when n is the product of small primes. 774 775 c a complex array of length n which contains the sequence 776 777 wsave 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 786 output parameters 787 788 c 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 796 wsave contains initialization calculations which must not be 797 destroyed between calls of subroutine cfftf or cfftb 798 799 */ 800