Julien Pommier | 370d209 | 2011-11-19 18:04:25 +0100 | [diff] [blame] | 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 | |
hayati ayguen | c974c1d | 2020-03-29 03:39:30 +0200 | [diff] [blame] | 52 | /* just define FFTPACK_DOUBLE_PRECISION if you want to build it as a double precision fft */ |
Julien Pommier | 370d209 | 2011-11-19 18:04:25 +0100 | [diff] [blame] | 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 | */ |