1! Fast Fourier/Cosine/Sine Transform 2! dimension :one 3! data length :power of 2 4! decimation :frequency 5! radix :split-radix 6! data :inplace 7! table :use 8! subroutines 9! cdft: Complex Discrete Fourier Transform 10! rdft: Real Discrete Fourier Transform 11! ddct: Discrete Cosine Transform 12! ddst: Discrete Sine Transform 13! dfct: Cosine Transform of RDFT (Real Symmetric DFT) 14! dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) 15! 16! 17! -------- Complex DFT (Discrete Fourier Transform) -------- 18! [definition] 19! <case1> 20! X(k) = sum_j=0^n-1 x(j)*exp(2*pi*i*j*k/n), 0<=k<n 21! <case2> 22! X(k) = sum_j=0^n-1 x(j)*exp(-2*pi*i*j*k/n), 0<=k<n 23! (notes: sum_j=0^n-1 is a summation from j=0 to n-1) 24! [usage] 25! <case1> 26! ip(0) = 0 ! first time only 27! call cdft(2*n, 1, a, ip, w) 28! <case2> 29! ip(0) = 0 ! first time only 30! call cdft(2*n, -1, a, ip, w) 31! [parameters] 32! 2*n :data length (integer) 33! n >= 1, n = power of 2 34! a(0:2*n-1) :input/output data (real*8) 35! input data 36! a(2*j) = Re(x(j)), 37! a(2*j+1) = Im(x(j)), 0<=j<n 38! output data 39! a(2*k) = Re(X(k)), 40! a(2*k+1) = Im(X(k)), 0<=k<n 41! ip(0:*) :work area for bit reversal (integer) 42! length of ip >= 2+sqrt(n) 43! strictly, 44! length of ip >= 45! 2+2**(int(log(n+0.5)/log(2.0))/2). 46! ip(0),ip(1) are pointers of the cos/sin table. 47! w(0:n/2-1) :cos/sin table (real*8) 48! w(),ip() are initialized if ip(0) = 0. 49! [remark] 50! Inverse of 51! call cdft(2*n, -1, a, ip, w) 52! is 53! call cdft(2*n, 1, a, ip, w) 54! do j = 0, 2 * n - 1 55! a(j) = a(j) / n 56! end do 57! . 58! 59! 60! -------- Real DFT / Inverse of Real DFT -------- 61! [definition] 62! <case1> RDFT 63! R(k) = sum_j=0^n-1 a(j)*cos(2*pi*j*k/n), 0<=k<=n/2 64! I(k) = sum_j=0^n-1 a(j)*sin(2*pi*j*k/n), 0<k<n/2 65! <case2> IRDFT (excluding scale) 66! a(k) = (R(0) + R(n/2)*cos(pi*k))/2 + 67! sum_j=1^n/2-1 R(j)*cos(2*pi*j*k/n) + 68! sum_j=1^n/2-1 I(j)*sin(2*pi*j*k/n), 0<=k<n 69! [usage] 70! <case1> 71! ip(0) = 0 ! first time only 72! call rdft(n, 1, a, ip, w) 73! <case2> 74! ip(0) = 0 ! first time only 75! call rdft(n, -1, a, ip, w) 76! [parameters] 77! n :data length (integer) 78! n >= 2, n = power of 2 79! a(0:n-1) :input/output data (real*8) 80! <case1> 81! output data 82! a(2*k) = R(k), 0<=k<n/2 83! a(2*k+1) = I(k), 0<k<n/2 84! a(1) = R(n/2) 85! <case2> 86! input data 87! a(2*j) = R(j), 0<=j<n/2 88! a(2*j+1) = I(j), 0<j<n/2 89! a(1) = R(n/2) 90! ip(0:*) :work area for bit reversal (integer) 91! length of ip >= 2+sqrt(n/2) 92! strictly, 93! length of ip >= 94! 2+2**(int(log(n/2+0.5)/log(2.0))/2). 95! ip(0),ip(1) are pointers of the cos/sin table. 96! w(0:n/2-1) :cos/sin table (real*8) 97! w(),ip() are initialized if ip(0) = 0. 98! [remark] 99! Inverse of 100! call rdft(n, 1, a, ip, w) 101! is 102! call rdft(n, -1, a, ip, w) 103! do j = 0, n - 1 104! a(j) = a(j) * 2 / n 105! end do 106! . 107! 108! 109! -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- 110! [definition] 111! <case1> IDCT (excluding scale) 112! C(k) = sum_j=0^n-1 a(j)*cos(pi*j*(k+1/2)/n), 0<=k<n 113! <case2> DCT 114! C(k) = sum_j=0^n-1 a(j)*cos(pi*(j+1/2)*k/n), 0<=k<n 115! [usage] 116! <case1> 117! ip(0) = 0 ! first time only 118! call ddct(n, 1, a, ip, w) 119! <case2> 120! ip(0) = 0 ! first time only 121! call ddct(n, -1, a, ip, w) 122! [parameters] 123! n :data length (integer) 124! n >= 2, n = power of 2 125! a(0:n-1) :input/output data (real*8) 126! output data 127! a(k) = C(k), 0<=k<n 128! ip(0:*) :work area for bit reversal (integer) 129! length of ip >= 2+sqrt(n/2) 130! strictly, 131! length of ip >= 132! 2+2**(int(log(n/2+0.5)/log(2.0))/2). 133! ip(0),ip(1) are pointers of the cos/sin table. 134! w(0:n*5/4-1) :cos/sin table (real*8) 135! w(),ip() are initialized if ip(0) = 0. 136! [remark] 137! Inverse of 138! call ddct(n, -1, a, ip, w) 139! is 140! a(0) = a(0) / 2 141! call ddct(n, 1, a, ip, w) 142! do j = 0, n - 1 143! a(j) = a(j) * 2 / n 144! end do 145! . 146! 147! 148! -------- DST (Discrete Sine Transform) / Inverse of DST -------- 149! [definition] 150! <case1> IDST (excluding scale) 151! S(k) = sum_j=1^n A(j)*sin(pi*j*(k+1/2)/n), 0<=k<n 152! <case2> DST 153! S(k) = sum_j=0^n-1 a(j)*sin(pi*(j+1/2)*k/n), 0<k<=n 154! [usage] 155! <case1> 156! ip(0) = 0 ! first time only 157! call ddst(n, 1, a, ip, w) 158! <case2> 159! ip(0) = 0 ! first time only 160! call ddst(n, -1, a, ip, w) 161! [parameters] 162! n :data length (integer) 163! n >= 2, n = power of 2 164! a(0:n-1) :input/output data (real*8) 165! <case1> 166! input data 167! a(j) = A(j), 0<j<n 168! a(0) = A(n) 169! output data 170! a(k) = S(k), 0<=k<n 171! <case2> 172! output data 173! a(k) = S(k), 0<k<n 174! a(0) = S(n) 175! ip(0:*) :work area for bit reversal (integer) 176! length of ip >= 2+sqrt(n/2) 177! strictly, 178! length of ip >= 179! 2+2**(int(log(n/2+0.5)/log(2.0))/2). 180! ip(0),ip(1) are pointers of the cos/sin table. 181! w(0:n*5/4-1) :cos/sin table (real*8) 182! w(),ip() are initialized if ip(0) = 0. 183! [remark] 184! Inverse of 185! call ddst(n, -1, a, ip, w) 186! is 187! a(0) = a(0) / 2 188! call ddst(n, 1, a, ip, w) 189! do j = 0, n - 1 190! a(j) = a(j) * 2 / n 191! end do 192! . 193! 194! 195! -------- Cosine Transform of RDFT (Real Symmetric DFT) -------- 196! [definition] 197! C(k) = sum_j=0^n a(j)*cos(pi*j*k/n), 0<=k<=n 198! [usage] 199! ip(0) = 0 ! first time only 200! call dfct(n, a, t, ip, w) 201! [parameters] 202! n :data length - 1 (integer) 203! n >= 2, n = power of 2 204! a(0:n) :input/output data (real*8) 205! output data 206! a(k) = C(k), 0<=k<=n 207! t(0:n/2) :work area (real*8) 208! ip(0:*) :work area for bit reversal (integer) 209! length of ip >= 2+sqrt(n/4) 210! strictly, 211! length of ip >= 212! 2+2**(int(log(n/4+0.5)/log(2.0))/2). 213! ip(0),ip(1) are pointers of the cos/sin table. 214! w(0:n*5/8-1) :cos/sin table (real*8) 215! w(),ip() are initialized if ip(0) = 0. 216! [remark] 217! Inverse of 218! a(0) = a(0) / 2 219! a(n) = a(n) / 2 220! call dfct(n, a, t, ip, w) 221! is 222! a(0) = a(0) / 2 223! a(n) = a(n) / 2 224! call dfct(n, a, t, ip, w) 225! do j = 0, n 226! a(j) = a(j) * 2 / n 227! end do 228! . 229! 230! 231! -------- Sine Transform of RDFT (Real Anti-symmetric DFT) -------- 232! [definition] 233! S(k) = sum_j=1^n-1 a(j)*sin(pi*j*k/n), 0<k<n 234! [usage] 235! ip(0) = 0 ! first time only 236! call dfst(n, a, t, ip, w) 237! [parameters] 238! n :data length + 1 (integer) 239! n >= 2, n = power of 2 240! a(0:n-1) :input/output data (real*8) 241! output data 242! a(k) = S(k), 0<k<n 243! (a(0) is used for work area) 244! t(0:n/2-1) :work area (real*8) 245! ip(0:*) :work area for bit reversal (integer) 246! length of ip >= 2+sqrt(n/4) 247! strictly, 248! length of ip >= 249! 2+2**(int(log(n/4+0.5)/log(2.0))/2). 250! ip(0),ip(1) are pointers of the cos/sin table. 251! w(0:n*5/8-1) :cos/sin table (real*8) 252! w(),ip() are initialized if ip(0) = 0. 253! [remark] 254! Inverse of 255! call dfst(n, a, t, ip, w) 256! is 257! call dfst(n, a, t, ip, w) 258! do j = 1, n - 1 259! a(j) = a(j) * 2 / n 260! end do 261! . 262! 263! 264! Appendix : 265! The cos/sin table is recalculated when the larger table required. 266! w() and ip() are compatible with all routines. 267! 268! 269 subroutine cdft(n, isgn, a, ip, w) 270 integer n, isgn, ip(0 : *), nw 271 real*8 a(0 : n - 1), w(0 : *) 272 nw = ip(0) 273 if (n .gt. 4 * nw) then 274 nw = n / 4 275 call makewt(nw, ip, w) 276 end if 277 if (isgn .ge. 0) then 278 call cftfsub(n, a, ip, nw, w) 279 else 280 call cftbsub(n, a, ip, nw, w) 281 end if 282 end 283! 284 subroutine rdft(n, isgn, a, ip, w) 285 integer n, isgn, ip(0 : *), nw, nc 286 real*8 a(0 : n - 1), w(0 : *), xi 287 nw = ip(0) 288 if (n .gt. 4 * nw) then 289 nw = n / 4 290 call makewt(nw, ip, w) 291 end if 292 nc = ip(1) 293 if (n .gt. 4 * nc) then 294 nc = n / 4 295 call makect(nc, ip, w(nw)) 296 end if 297 if (isgn .ge. 0) then 298 if (n .gt. 4) then 299 call cftfsub(n, a, ip, nw, w) 300 call rftfsub(n, a, nc, w(nw)) 301 else if (n .eq. 4) then 302 call cftfsub(n, a, ip, nw, w) 303 end if 304 xi = a(0) - a(1) 305 a(0) = a(0) + a(1) 306 a(1) = xi 307 else 308 a(1) = 0.5d0 * (a(0) - a(1)) 309 a(0) = a(0) - a(1) 310 if (n .gt. 4) then 311 call rftbsub(n, a, nc, w(nw)) 312 call cftbsub(n, a, ip, nw, w) 313 else if (n .eq. 4) then 314 call cftbsub(n, a, ip, nw, w) 315 end if 316 end if 317 end 318! 319 subroutine ddct(n, isgn, a, ip, w) 320 integer n, isgn, ip(0 : *), j, nw, nc 321 real*8 a(0 : n - 1), w(0 : *), xr 322 nw = ip(0) 323 if (n .gt. 4 * nw) then 324 nw = n / 4 325 call makewt(nw, ip, w) 326 end if 327 nc = ip(1) 328 if (n .gt. nc) then 329 nc = n 330 call makect(nc, ip, w(nw)) 331 end if 332 if (isgn .lt. 0) then 333 xr = a(n - 1) 334 do j = n - 2, 2, -2 335 a(j + 1) = a(j) - a(j - 1) 336 a(j) = a(j) + a(j - 1) 337 end do 338 a(1) = a(0) - xr 339 a(0) = a(0) + xr 340 if (n .gt. 4) then 341 call rftbsub(n, a, nc, w(nw)) 342 call cftbsub(n, a, ip, nw, w) 343 else if (n .eq. 4) then 344 call cftbsub(n, a, ip, nw, w) 345 end if 346 end if 347 call dctsub(n, a, nc, w(nw)) 348 if (isgn .ge. 0) then 349 if (n .gt. 4) then 350 call cftfsub(n, a, ip, nw, w) 351 call rftfsub(n, a, nc, w(nw)) 352 else if (n .eq. 4) then 353 call cftfsub(n, a, ip, nw, w) 354 end if 355 xr = a(0) - a(1) 356 a(0) = a(0) + a(1) 357 do j = 2, n - 2, 2 358 a(j - 1) = a(j) - a(j + 1) 359 a(j) = a(j) + a(j + 1) 360 end do 361 a(n - 1) = xr 362 end if 363 end 364! 365 subroutine ddst(n, isgn, a, ip, w) 366 integer n, isgn, ip(0 : *), j, nw, nc 367 real*8 a(0 : n - 1), w(0 : *), xr 368 nw = ip(0) 369 if (n .gt. 4 * nw) then 370 nw = n / 4 371 call makewt(nw, ip, w) 372 end if 373 nc = ip(1) 374 if (n .gt. nc) then 375 nc = n 376 call makect(nc, ip, w(nw)) 377 end if 378 if (isgn .lt. 0) then 379 xr = a(n - 1) 380 do j = n - 2, 2, -2 381 a(j + 1) = -a(j) - a(j - 1) 382 a(j) = a(j) - a(j - 1) 383 end do 384 a(1) = a(0) + xr 385 a(0) = a(0) - xr 386 if (n .gt. 4) then 387 call rftbsub(n, a, nc, w(nw)) 388 call cftbsub(n, a, ip, nw, w) 389 else if (n .eq. 4) then 390 call cftbsub(n, a, ip, nw, w) 391 end if 392 end if 393 call dstsub(n, a, nc, w(nw)) 394 if (isgn .ge. 0) then 395 if (n .gt. 4) then 396 call cftfsub(n, a, ip, nw, w) 397 call rftfsub(n, a, nc, w(nw)) 398 else if (n .eq. 4) then 399 call cftfsub(n, a, ip, nw, w) 400 end if 401 xr = a(0) - a(1) 402 a(0) = a(0) + a(1) 403 do j = 2, n - 2, 2 404 a(j - 1) = -a(j) - a(j + 1) 405 a(j) = a(j) - a(j + 1) 406 end do 407 a(n - 1) = -xr 408 end if 409 end 410! 411 subroutine dfct(n, a, t, ip, w) 412 integer n, ip(0 : *), j, k, l, m, mh, nw, nc 413 real*8 a(0 : n), t(0 : n / 2), w(0 : *), xr, xi, yr, yi 414 nw = ip(0) 415 if (n .gt. 8 * nw) then 416 nw = n / 8 417 call makewt(nw, ip, w) 418 end if 419 nc = ip(1) 420 if (n .gt. 2 * nc) then 421 nc = n / 2 422 call makect(nc, ip, w(nw)) 423 end if 424 m = n / 2 425 yi = a(m) 426 xi = a(0) + a(n) 427 a(0) = a(0) - a(n) 428 t(0) = xi - yi 429 t(m) = xi + yi 430 if (n .gt. 2) then 431 mh = m / 2 432 do j = 1, mh - 1 433 k = m - j 434 xr = a(j) - a(n - j) 435 xi = a(j) + a(n - j) 436 yr = a(k) - a(n - k) 437 yi = a(k) + a(n - k) 438 a(j) = xr 439 a(k) = yr 440 t(j) = xi - yi 441 t(k) = xi + yi 442 end do 443 t(mh) = a(mh) + a(n - mh) 444 a(mh) = a(mh) - a(n - mh) 445 call dctsub(m, a, nc, w(nw)) 446 if (m .gt. 4) then 447 call cftfsub(m, a, ip, nw, w) 448 call rftfsub(m, a, nc, w(nw)) 449 else if (m .eq. 4) then 450 call cftfsub(m, a, ip, nw, w) 451 end if 452 a(n - 1) = a(0) - a(1) 453 a(1) = a(0) + a(1) 454 do j = m - 2, 2, -2 455 a(2 * j + 1) = a(j) + a(j + 1) 456 a(2 * j - 1) = a(j) - a(j + 1) 457 end do 458 l = 2 459 m = mh 460 do while (m .ge. 2) 461 call dctsub(m, t, nc, w(nw)) 462 if (m .gt. 4) then 463 call cftfsub(m, t, ip, nw, w) 464 call rftfsub(m, t, nc, w(nw)) 465 else if (m .eq. 4) then 466 call cftfsub(m, t, ip, nw, w) 467 end if 468 a(n - l) = t(0) - t(1) 469 a(l) = t(0) + t(1) 470 k = 0 471 do j = 2, m - 2, 2 472 k = k + 4 * l 473 a(k - l) = t(j) - t(j + 1) 474 a(k + l) = t(j) + t(j + 1) 475 end do 476 l = 2 * l 477 mh = m / 2 478 do j = 0, mh - 1 479 k = m - j 480 t(j) = t(m + k) - t(m + j) 481 t(k) = t(m + k) + t(m + j) 482 end do 483 t(mh) = t(m + mh) 484 m = mh 485 end do 486 a(l) = t(0) 487 a(n) = t(2) - t(1) 488 a(0) = t(2) + t(1) 489 else 490 a(1) = a(0) 491 a(2) = t(0) 492 a(0) = t(1) 493 end if 494 end 495! 496 subroutine dfst(n, a, t, ip, w) 497 integer n, ip(0 : *), j, k, l, m, mh, nw, nc 498 real*8 a(0 : n - 1), t(0 : n / 2 - 1), w(0 : *), xr, xi, yr, yi 499 nw = ip(0) 500 if (n .gt. 8 * nw) then 501 nw = n / 8 502 call makewt(nw, ip, w) 503 end if 504 nc = ip(1) 505 if (n .gt. 2 * nc) then 506 nc = n / 2 507 call makect(nc, ip, w(nw)) 508 end if 509 if (n .gt. 2) then 510 m = n / 2 511 mh = m / 2 512 do j = 1, mh - 1 513 k = m - j 514 xr = a(j) + a(n - j) 515 xi = a(j) - a(n - j) 516 yr = a(k) + a(n - k) 517 yi = a(k) - a(n - k) 518 a(j) = xr 519 a(k) = yr 520 t(j) = xi + yi 521 t(k) = xi - yi 522 end do 523 t(0) = a(mh) - a(n - mh) 524 a(mh) = a(mh) + a(n - mh) 525 a(0) = a(m) 526 call dstsub(m, a, nc, w(nw)) 527 if (m .gt. 4) then 528 call cftfsub(m, a, ip, nw, w) 529 call rftfsub(m, a, nc, w(nw)) 530 else if (m .eq. 4) then 531 call cftfsub(m, a, ip, nw, w) 532 end if 533 a(n - 1) = a(1) - a(0) 534 a(1) = a(0) + a(1) 535 do j = m - 2, 2, -2 536 a(2 * j + 1) = a(j) - a(j + 1) 537 a(2 * j - 1) = -a(j) - a(j + 1) 538 end do 539 l = 2 540 m = mh 541 do while (m .ge. 2) 542 call dstsub(m, t, nc, w(nw)) 543 if (m .gt. 4) then 544 call cftfsub(m, t, ip, nw, w) 545 call rftfsub(m, t, nc, w(nw)) 546 else if (m .eq. 4) then 547 call cftfsub(m, t, ip, nw, w) 548 end if 549 a(n - l) = t(1) - t(0) 550 a(l) = t(0) + t(1) 551 k = 0 552 do j = 2, m - 2, 2 553 k = k + 4 * l 554 a(k - l) = -t(j) - t(j + 1) 555 a(k + l) = t(j) - t(j + 1) 556 end do 557 l = 2 * l 558 mh = m / 2 559 do j = 1, mh - 1 560 k = m - j 561 t(j) = t(m + k) + t(m + j) 562 t(k) = t(m + k) - t(m + j) 563 end do 564 t(0) = t(m + mh) 565 m = mh 566 end do 567 a(l) = t(0) 568 end if 569 a(0) = 0 570 end 571! 572! -------- initializing routines -------- 573! 574 subroutine makewt(nw, ip, w) 575 integer nw, ip(0 : *), j, nwh, nw0, nw1 576 real*8 w(0 : nw - 1), delta, wn4r, wk1r, wk1i, wk3r, wk3i 577 ip(0) = nw 578 ip(1) = 1 579 if (nw .gt. 2) then 580 nwh = nw / 2 581 delta = atan(1.0d0) / nwh 582 wn4r = cos(delta * nwh) 583 w(0) = 1 584 w(1) = wn4r 585 if (nwh .eq. 4) then 586 w(2) = cos(delta * 2) 587 w(3) = sin(delta * 2) 588 else if (nwh .gt. 4) then 589 call makeipt(nw, ip) 590 w(2) = 0.5d0 / cos(delta * 2) 591 w(3) = 0.5d0 / cos(delta * 6) 592 do j = 4, nwh - 4, 4 593 w(j) = cos(delta * j) 594 w(j + 1) = sin(delta * j) 595 w(j + 2) = cos(3 * delta * j) 596 w(j + 3) = -sin(3 * delta * j) 597 end do 598 end if 599 nw0 = 0 600 do while (nwh .gt. 2) 601 nw1 = nw0 + nwh 602 nwh = nwh / 2 603 w(nw1) = 1 604 w(nw1 + 1) = wn4r 605 if (nwh .eq. 4) then 606 wk1r = w(nw0 + 4) 607 wk1i = w(nw0 + 5) 608 w(nw1 + 2) = wk1r 609 w(nw1 + 3) = wk1i 610 else if (nwh .gt. 4) then 611 wk1r = w(nw0 + 4) 612 wk3r = w(nw0 + 6) 613 w(nw1 + 2) = 0.5d0 / wk1r 614 w(nw1 + 3) = 0.5d0 / wk3r 615 do j = 4, nwh - 4, 4 616 wk1r = w(nw0 + 2 * j) 617 wk1i = w(nw0 + 2 * j + 1) 618 wk3r = w(nw0 + 2 * j + 2) 619 wk3i = w(nw0 + 2 * j + 3) 620 w(nw1 + j) = wk1r 621 w(nw1 + j + 1) = wk1i 622 w(nw1 + j + 2) = wk3r 623 w(nw1 + j + 3) = wk3i 624 end do 625 end if 626 nw0 = nw1 627 end do 628 end if 629 end 630! 631 subroutine makeipt(nw, ip) 632 integer nw, ip(0 : *), j, l, m, m2, p, q 633 ip(2) = 0 634 ip(3) = 16 635 m = 2 636 l = nw 637 do while (l .gt. 32) 638 m2 = 2 * m 639 q = 8 * m2 640 do j = m, m2 - 1 641 p = 4 * ip(j) 642 ip(m + j) = p 643 ip(m2 + j) = p + q 644 end do 645 m = m2 646 l = l / 4 647 end do 648 end 649! 650 subroutine makect(nc, ip, c) 651 integer nc, ip(0 : *), j, nch 652 real*8 c(0 : nc - 1), delta 653 ip(1) = nc 654 if (nc .gt. 1) then 655 nch = nc / 2 656 delta = atan(1.0d0) / nch 657 c(0) = cos(delta * nch) 658 c(nch) = 0.5d0 * c(0) 659 do j = 1, nch - 1 660 c(j) = 0.5d0 * cos(delta * j) 661 c(nc - j) = 0.5d0 * sin(delta * j) 662 end do 663 end if 664 end 665! 666! -------- child routines -------- 667! 668 subroutine cftfsub(n, a, ip, nw, w) 669 integer n, ip(0 : *), nw 670 real*8 a(0 : n - 1), w(0 : nw - 1) 671 if (n .gt. 8) then 672 if (n .gt. 32) then 673 call cftf1st(n, a, w(nw - n / 4)) 674 if (n .gt. 512) then 675 call cftrec4(n, a, nw, w) 676 else if (n .gt. 128) then 677 call cftleaf(n, 1, a, nw, w) 678 else 679 call cftfx41(n, a, nw, w) 680 end if 681 call bitrv2(n, ip, a) 682 else if (n .eq. 32) then 683 call cftf161(a, w(nw - 8)) 684 call bitrv216(a) 685 else 686 call cftf081(a, w) 687 call bitrv208(a) 688 end if 689 else if (n .eq. 8) then 690 call cftf040(a) 691 else if (n .eq. 4) then 692 call cftx020(a) 693 end if 694 end 695! 696 subroutine cftbsub(n, a, ip, nw, w) 697 integer n, ip(0 : *), nw 698 real*8 a(0 : n - 1), w(0 : nw - 1) 699 if (n .gt. 8) then 700 if (n .gt. 32) then 701 call cftb1st(n, a, w(nw - n / 4)) 702 if (n .gt. 512) then 703 call cftrec4(n, a, nw, w) 704 else if (n .gt. 128) then 705 call cftleaf(n, 1, a, nw, w) 706 else 707 call cftfx41(n, a, nw, w) 708 end if 709 call bitrv2conj(n, ip, a) 710 else if (n .eq. 32) then 711 call cftf161(a, w(nw - 8)) 712 call bitrv216neg(a) 713 else 714 call cftf081(a, w) 715 call bitrv208neg(a) 716 end if 717 else if (n .eq. 8) then 718 call cftb040(a) 719 else if (n .eq. 4) then 720 call cftx020(a) 721 end if 722 end 723! 724 subroutine bitrv2(n, ip, a) 725 integer n, ip(0 : *), j, j1, k, k1, l, m, nh, nm 726 real*8 a(0 : n - 1), xr, xi, yr, yi 727 m = 1 728 l = n / 4 729 do while (l .gt. 8) 730 m = m * 2 731 l = l / 4 732 end do 733 nh = n / 2 734 nm = 4 * m 735 if (l .eq. 8) then 736 do k = 0, m - 1 737 do j = 0, k - 1 738 j1 = 4 * j + 2 * ip(m + k) 739 k1 = 4 * k + 2 * ip(m + j) 740 xr = a(j1) 741 xi = a(j1 + 1) 742 yr = a(k1) 743 yi = a(k1 + 1) 744 a(j1) = yr 745 a(j1 + 1) = yi 746 a(k1) = xr 747 a(k1 + 1) = xi 748 j1 = j1 + nm 749 k1 = k1 + 2 * nm 750 xr = a(j1) 751 xi = a(j1 + 1) 752 yr = a(k1) 753 yi = a(k1 + 1) 754 a(j1) = yr 755 a(j1 + 1) = yi 756 a(k1) = xr 757 a(k1 + 1) = xi 758 j1 = j1 + nm 759 k1 = k1 - nm 760 xr = a(j1) 761 xi = a(j1 + 1) 762 yr = a(k1) 763 yi = a(k1 + 1) 764 a(j1) = yr 765 a(j1 + 1) = yi 766 a(k1) = xr 767 a(k1 + 1) = xi 768 j1 = j1 + nm 769 k1 = k1 + 2 * nm 770 xr = a(j1) 771 xi = a(j1 + 1) 772 yr = a(k1) 773 yi = a(k1 + 1) 774 a(j1) = yr 775 a(j1 + 1) = yi 776 a(k1) = xr 777 a(k1 + 1) = xi 778 j1 = j1 + nh 779 k1 = k1 + 2 780 xr = a(j1) 781 xi = a(j1 + 1) 782 yr = a(k1) 783 yi = a(k1 + 1) 784 a(j1) = yr 785 a(j1 + 1) = yi 786 a(k1) = xr 787 a(k1 + 1) = xi 788 j1 = j1 - nm 789 k1 = k1 - 2 * nm 790 xr = a(j1) 791 xi = a(j1 + 1) 792 yr = a(k1) 793 yi = a(k1 + 1) 794 a(j1) = yr 795 a(j1 + 1) = yi 796 a(k1) = xr 797 a(k1 + 1) = xi 798 j1 = j1 - nm 799 k1 = k1 + nm 800 xr = a(j1) 801 xi = a(j1 + 1) 802 yr = a(k1) 803 yi = a(k1 + 1) 804 a(j1) = yr 805 a(j1 + 1) = yi 806 a(k1) = xr 807 a(k1 + 1) = xi 808 j1 = j1 - nm 809 k1 = k1 - 2 * nm 810 xr = a(j1) 811 xi = a(j1 + 1) 812 yr = a(k1) 813 yi = a(k1 + 1) 814 a(j1) = yr 815 a(j1 + 1) = yi 816 a(k1) = xr 817 a(k1 + 1) = xi 818 j1 = j1 + 2 819 k1 = k1 + nh 820 xr = a(j1) 821 xi = a(j1 + 1) 822 yr = a(k1) 823 yi = a(k1 + 1) 824 a(j1) = yr 825 a(j1 + 1) = yi 826 a(k1) = xr 827 a(k1 + 1) = xi 828 j1 = j1 + nm 829 k1 = k1 + 2 * nm 830 xr = a(j1) 831 xi = a(j1 + 1) 832 yr = a(k1) 833 yi = a(k1 + 1) 834 a(j1) = yr 835 a(j1 + 1) = yi 836 a(k1) = xr 837 a(k1 + 1) = xi 838 j1 = j1 + nm 839 k1 = k1 - nm 840 xr = a(j1) 841 xi = a(j1 + 1) 842 yr = a(k1) 843 yi = a(k1 + 1) 844 a(j1) = yr 845 a(j1 + 1) = yi 846 a(k1) = xr 847 a(k1 + 1) = xi 848 j1 = j1 + nm 849 k1 = k1 + 2 * nm 850 xr = a(j1) 851 xi = a(j1 + 1) 852 yr = a(k1) 853 yi = a(k1 + 1) 854 a(j1) = yr 855 a(j1 + 1) = yi 856 a(k1) = xr 857 a(k1 + 1) = xi 858 j1 = j1 - nh 859 k1 = k1 - 2 860 xr = a(j1) 861 xi = a(j1 + 1) 862 yr = a(k1) 863 yi = a(k1 + 1) 864 a(j1) = yr 865 a(j1 + 1) = yi 866 a(k1) = xr 867 a(k1 + 1) = xi 868 j1 = j1 - nm 869 k1 = k1 - 2 * nm 870 xr = a(j1) 871 xi = a(j1 + 1) 872 yr = a(k1) 873 yi = a(k1 + 1) 874 a(j1) = yr 875 a(j1 + 1) = yi 876 a(k1) = xr 877 a(k1 + 1) = xi 878 j1 = j1 - nm 879 k1 = k1 + nm 880 xr = a(j1) 881 xi = a(j1 + 1) 882 yr = a(k1) 883 yi = a(k1 + 1) 884 a(j1) = yr 885 a(j1 + 1) = yi 886 a(k1) = xr 887 a(k1 + 1) = xi 888 j1 = j1 - nm 889 k1 = k1 - 2 * nm 890 xr = a(j1) 891 xi = a(j1 + 1) 892 yr = a(k1) 893 yi = a(k1 + 1) 894 a(j1) = yr 895 a(j1 + 1) = yi 896 a(k1) = xr 897 a(k1 + 1) = xi 898 end do 899 k1 = 4 * k + 2 * ip(m + k) 900 j1 = k1 + 2 901 k1 = k1 + nh 902 xr = a(j1) 903 xi = a(j1 + 1) 904 yr = a(k1) 905 yi = a(k1 + 1) 906 a(j1) = yr 907 a(j1 + 1) = yi 908 a(k1) = xr 909 a(k1 + 1) = xi 910 j1 = j1 + nm 911 k1 = k1 + 2 * nm 912 xr = a(j1) 913 xi = a(j1 + 1) 914 yr = a(k1) 915 yi = a(k1 + 1) 916 a(j1) = yr 917 a(j1 + 1) = yi 918 a(k1) = xr 919 a(k1 + 1) = xi 920 j1 = j1 + nm 921 k1 = k1 - nm 922 xr = a(j1) 923 xi = a(j1 + 1) 924 yr = a(k1) 925 yi = a(k1 + 1) 926 a(j1) = yr 927 a(j1 + 1) = yi 928 a(k1) = xr 929 a(k1 + 1) = xi 930 j1 = j1 - 2 931 k1 = k1 - nh 932 xr = a(j1) 933 xi = a(j1 + 1) 934 yr = a(k1) 935 yi = a(k1 + 1) 936 a(j1) = yr 937 a(j1 + 1) = yi 938 a(k1) = xr 939 a(k1 + 1) = xi 940 j1 = j1 + nh + 2 941 k1 = k1 + nh + 2 942 xr = a(j1) 943 xi = a(j1 + 1) 944 yr = a(k1) 945 yi = a(k1 + 1) 946 a(j1) = yr 947 a(j1 + 1) = yi 948 a(k1) = xr 949 a(k1 + 1) = xi 950 j1 = j1 - nh + nm 951 k1 = k1 + 2 * nm - 2 952 xr = a(j1) 953 xi = a(j1 + 1) 954 yr = a(k1) 955 yi = a(k1 + 1) 956 a(j1) = yr 957 a(j1 + 1) = yi 958 a(k1) = xr 959 a(k1 + 1) = xi 960 end do 961 else 962 do k = 0, m - 1 963 do j = 0, k - 1 964 j1 = 4 * j + ip(m + k) 965 k1 = 4 * k + ip(m + j) 966 xr = a(j1) 967 xi = a(j1 + 1) 968 yr = a(k1) 969 yi = a(k1 + 1) 970 a(j1) = yr 971 a(j1 + 1) = yi 972 a(k1) = xr 973 a(k1 + 1) = xi 974 j1 = j1 + nm 975 k1 = k1 + nm 976 xr = a(j1) 977 xi = a(j1 + 1) 978 yr = a(k1) 979 yi = a(k1 + 1) 980 a(j1) = yr 981 a(j1 + 1) = yi 982 a(k1) = xr 983 a(k1 + 1) = xi 984 j1 = j1 + nh 985 k1 = k1 + 2 986 xr = a(j1) 987 xi = a(j1 + 1) 988 yr = a(k1) 989 yi = a(k1 + 1) 990 a(j1) = yr 991 a(j1 + 1) = yi 992 a(k1) = xr 993 a(k1 + 1) = xi 994 j1 = j1 - nm 995 k1 = k1 - nm 996 xr = a(j1) 997 xi = a(j1 + 1) 998 yr = a(k1) 999 yi = a(k1 + 1) 1000 a(j1) = yr 1001 a(j1 + 1) = yi 1002 a(k1) = xr 1003 a(k1 + 1) = xi 1004 j1 = j1 + 2 1005 k1 = k1 + nh 1006 xr = a(j1) 1007 xi = a(j1 + 1) 1008 yr = a(k1) 1009 yi = a(k1 + 1) 1010 a(j1) = yr 1011 a(j1 + 1) = yi 1012 a(k1) = xr 1013 a(k1 + 1) = xi 1014 j1 = j1 + nm 1015 k1 = k1 + nm 1016 xr = a(j1) 1017 xi = a(j1 + 1) 1018 yr = a(k1) 1019 yi = a(k1 + 1) 1020 a(j1) = yr 1021 a(j1 + 1) = yi 1022 a(k1) = xr 1023 a(k1 + 1) = xi 1024 j1 = j1 - nh 1025 k1 = k1 - 2 1026 xr = a(j1) 1027 xi = a(j1 + 1) 1028 yr = a(k1) 1029 yi = a(k1 + 1) 1030 a(j1) = yr 1031 a(j1 + 1) = yi 1032 a(k1) = xr 1033 a(k1 + 1) = xi 1034 j1 = j1 - nm 1035 k1 = k1 - nm 1036 xr = a(j1) 1037 xi = a(j1 + 1) 1038 yr = a(k1) 1039 yi = a(k1 + 1) 1040 a(j1) = yr 1041 a(j1 + 1) = yi 1042 a(k1) = xr 1043 a(k1 + 1) = xi 1044 end do 1045 k1 = 4 * k + ip(m + k) 1046 j1 = k1 + 2 1047 k1 = k1 + nh 1048 xr = a(j1) 1049 xi = a(j1 + 1) 1050 yr = a(k1) 1051 yi = a(k1 + 1) 1052 a(j1) = yr 1053 a(j1 + 1) = yi 1054 a(k1) = xr 1055 a(k1 + 1) = xi 1056 j1 = j1 + nm 1057 k1 = k1 + nm 1058 xr = a(j1) 1059 xi = a(j1 + 1) 1060 yr = a(k1) 1061 yi = a(k1 + 1) 1062 a(j1) = yr 1063 a(j1 + 1) = yi 1064 a(k1) = xr 1065 a(k1 + 1) = xi 1066 end do 1067 end if 1068 end 1069! 1070 subroutine bitrv2conj(n, ip, a) 1071 integer n, ip(0 : *), j, j1, k, k1, l, m, nh, nm 1072 real*8 a(0 : n - 1), xr, xi, yr, yi 1073 m = 1 1074 l = n / 4 1075 do while (l .gt. 8) 1076 m = m * 2 1077 l = l / 4 1078 end do 1079 nh = n / 2 1080 nm = 4 * m 1081 if (l .eq. 8) then 1082 do k = 0, m - 1 1083 do j = 0, k - 1 1084 j1 = 4 * j + 2 * ip(m + k) 1085 k1 = 4 * k + 2 * ip(m + j) 1086 xr = a(j1) 1087 xi = -a(j1 + 1) 1088 yr = a(k1) 1089 yi = -a(k1 + 1) 1090 a(j1) = yr 1091 a(j1 + 1) = yi 1092 a(k1) = xr 1093 a(k1 + 1) = xi 1094 j1 = j1 + nm 1095 k1 = k1 + 2 * nm 1096 xr = a(j1) 1097 xi = -a(j1 + 1) 1098 yr = a(k1) 1099 yi = -a(k1 + 1) 1100 a(j1) = yr 1101 a(j1 + 1) = yi 1102 a(k1) = xr 1103 a(k1 + 1) = xi 1104 j1 = j1 + nm 1105 k1 = k1 - nm 1106 xr = a(j1) 1107 xi = -a(j1 + 1) 1108 yr = a(k1) 1109 yi = -a(k1 + 1) 1110 a(j1) = yr 1111 a(j1 + 1) = yi 1112 a(k1) = xr 1113 a(k1 + 1) = xi 1114 j1 = j1 + nm 1115 k1 = k1 + 2 * nm 1116 xr = a(j1) 1117 xi = -a(j1 + 1) 1118 yr = a(k1) 1119 yi = -a(k1 + 1) 1120 a(j1) = yr 1121 a(j1 + 1) = yi 1122 a(k1) = xr 1123 a(k1 + 1) = xi 1124 j1 = j1 + nh 1125 k1 = k1 + 2 1126 xr = a(j1) 1127 xi = -a(j1 + 1) 1128 yr = a(k1) 1129 yi = -a(k1 + 1) 1130 a(j1) = yr 1131 a(j1 + 1) = yi 1132 a(k1) = xr 1133 a(k1 + 1) = xi 1134 j1 = j1 - nm 1135 k1 = k1 - 2 * nm 1136 xr = a(j1) 1137 xi = -a(j1 + 1) 1138 yr = a(k1) 1139 yi = -a(k1 + 1) 1140 a(j1) = yr 1141 a(j1 + 1) = yi 1142 a(k1) = xr 1143 a(k1 + 1) = xi 1144 j1 = j1 - nm 1145 k1 = k1 + nm 1146 xr = a(j1) 1147 xi = -a(j1 + 1) 1148 yr = a(k1) 1149 yi = -a(k1 + 1) 1150 a(j1) = yr 1151 a(j1 + 1) = yi 1152 a(k1) = xr 1153 a(k1 + 1) = xi 1154 j1 = j1 - nm 1155 k1 = k1 - 2 * nm 1156 xr = a(j1) 1157 xi = -a(j1 + 1) 1158 yr = a(k1) 1159 yi = -a(k1 + 1) 1160 a(j1) = yr 1161 a(j1 + 1) = yi 1162 a(k1) = xr 1163 a(k1 + 1) = xi 1164 j1 = j1 + 2 1165 k1 = k1 + nh 1166 xr = a(j1) 1167 xi = -a(j1 + 1) 1168 yr = a(k1) 1169 yi = -a(k1 + 1) 1170 a(j1) = yr 1171 a(j1 + 1) = yi 1172 a(k1) = xr 1173 a(k1 + 1) = xi 1174 j1 = j1 + nm 1175 k1 = k1 + 2 * nm 1176 xr = a(j1) 1177 xi = -a(j1 + 1) 1178 yr = a(k1) 1179 yi = -a(k1 + 1) 1180 a(j1) = yr 1181 a(j1 + 1) = yi 1182 a(k1) = xr 1183 a(k1 + 1) = xi 1184 j1 = j1 + nm 1185 k1 = k1 - nm 1186 xr = a(j1) 1187 xi = -a(j1 + 1) 1188 yr = a(k1) 1189 yi = -a(k1 + 1) 1190 a(j1) = yr 1191 a(j1 + 1) = yi 1192 a(k1) = xr 1193 a(k1 + 1) = xi 1194 j1 = j1 + nm 1195 k1 = k1 + 2 * nm 1196 xr = a(j1) 1197 xi = -a(j1 + 1) 1198 yr = a(k1) 1199 yi = -a(k1 + 1) 1200 a(j1) = yr 1201 a(j1 + 1) = yi 1202 a(k1) = xr 1203 a(k1 + 1) = xi 1204 j1 = j1 - nh 1205 k1 = k1 - 2 1206 xr = a(j1) 1207 xi = -a(j1 + 1) 1208 yr = a(k1) 1209 yi = -a(k1 + 1) 1210 a(j1) = yr 1211 a(j1 + 1) = yi 1212 a(k1) = xr 1213 a(k1 + 1) = xi 1214 j1 = j1 - nm 1215 k1 = k1 - 2 * nm 1216 xr = a(j1) 1217 xi = -a(j1 + 1) 1218 yr = a(k1) 1219 yi = -a(k1 + 1) 1220 a(j1) = yr 1221 a(j1 + 1) = yi 1222 a(k1) = xr 1223 a(k1 + 1) = xi 1224 j1 = j1 - nm 1225 k1 = k1 + nm 1226 xr = a(j1) 1227 xi = -a(j1 + 1) 1228 yr = a(k1) 1229 yi = -a(k1 + 1) 1230 a(j1) = yr 1231 a(j1 + 1) = yi 1232 a(k1) = xr 1233 a(k1 + 1) = xi 1234 j1 = j1 - nm 1235 k1 = k1 - 2 * nm 1236 xr = a(j1) 1237 xi = -a(j1 + 1) 1238 yr = a(k1) 1239 yi = -a(k1 + 1) 1240 a(j1) = yr 1241 a(j1 + 1) = yi 1242 a(k1) = xr 1243 a(k1 + 1) = xi 1244 end do 1245 k1 = 4 * k + 2 * ip(m + k) 1246 j1 = k1 + 2 1247 k1 = k1 + nh 1248 a(j1 - 1) = -a(j1 - 1) 1249 xr = a(j1) 1250 xi = -a(j1 + 1) 1251 yr = a(k1) 1252 yi = -a(k1 + 1) 1253 a(j1) = yr 1254 a(j1 + 1) = yi 1255 a(k1) = xr 1256 a(k1 + 1) = xi 1257 a(k1 + 3) = -a(k1 + 3) 1258 j1 = j1 + nm 1259 k1 = k1 + 2 * nm 1260 xr = a(j1) 1261 xi = -a(j1 + 1) 1262 yr = a(k1) 1263 yi = -a(k1 + 1) 1264 a(j1) = yr 1265 a(j1 + 1) = yi 1266 a(k1) = xr 1267 a(k1 + 1) = xi 1268 j1 = j1 + nm 1269 k1 = k1 - nm 1270 xr = a(j1) 1271 xi = -a(j1 + 1) 1272 yr = a(k1) 1273 yi = -a(k1 + 1) 1274 a(j1) = yr 1275 a(j1 + 1) = yi 1276 a(k1) = xr 1277 a(k1 + 1) = xi 1278 j1 = j1 - 2 1279 k1 = k1 - nh 1280 xr = a(j1) 1281 xi = -a(j1 + 1) 1282 yr = a(k1) 1283 yi = -a(k1 + 1) 1284 a(j1) = yr 1285 a(j1 + 1) = yi 1286 a(k1) = xr 1287 a(k1 + 1) = xi 1288 j1 = j1 + nh + 2 1289 k1 = k1 + nh + 2 1290 xr = a(j1) 1291 xi = -a(j1 + 1) 1292 yr = a(k1) 1293 yi = -a(k1 + 1) 1294 a(j1) = yr 1295 a(j1 + 1) = yi 1296 a(k1) = xr 1297 a(k1 + 1) = xi 1298 j1 = j1 - nh + nm 1299 k1 = k1 + 2 * nm - 2 1300 a(j1 - 1) = -a(j1 - 1) 1301 xr = a(j1) 1302 xi = -a(j1 + 1) 1303 yr = a(k1) 1304 yi = -a(k1 + 1) 1305 a(j1) = yr 1306 a(j1 + 1) = yi 1307 a(k1) = xr 1308 a(k1 + 1) = xi 1309 a(k1 + 3) = -a(k1 + 3) 1310 end do 1311 else 1312 do k = 0, m - 1 1313 do j = 0, k - 1 1314 j1 = 4 * j + ip(m + k) 1315 k1 = 4 * k + ip(m + j) 1316 xr = a(j1) 1317 xi = -a(j1 + 1) 1318 yr = a(k1) 1319 yi = -a(k1 + 1) 1320 a(j1) = yr 1321 a(j1 + 1) = yi 1322 a(k1) = xr 1323 a(k1 + 1) = xi 1324 j1 = j1 + nm 1325 k1 = k1 + nm 1326 xr = a(j1) 1327 xi = -a(j1 + 1) 1328 yr = a(k1) 1329 yi = -a(k1 + 1) 1330 a(j1) = yr 1331 a(j1 + 1) = yi 1332 a(k1) = xr 1333 a(k1 + 1) = xi 1334 j1 = j1 + nh 1335 k1 = k1 + 2 1336 xr = a(j1) 1337 xi = -a(j1 + 1) 1338 yr = a(k1) 1339 yi = -a(k1 + 1) 1340 a(j1) = yr 1341 a(j1 + 1) = yi 1342 a(k1) = xr 1343 a(k1 + 1) = xi 1344 j1 = j1 - nm 1345 k1 = k1 - nm 1346 xr = a(j1) 1347 xi = -a(j1 + 1) 1348 yr = a(k1) 1349 yi = -a(k1 + 1) 1350 a(j1) = yr 1351 a(j1 + 1) = yi 1352 a(k1) = xr 1353 a(k1 + 1) = xi 1354 j1 = j1 + 2 1355 k1 = k1 + nh 1356 xr = a(j1) 1357 xi = -a(j1 + 1) 1358 yr = a(k1) 1359 yi = -a(k1 + 1) 1360 a(j1) = yr 1361 a(j1 + 1) = yi 1362 a(k1) = xr 1363 a(k1 + 1) = xi 1364 j1 = j1 + nm 1365 k1 = k1 + nm 1366 xr = a(j1) 1367 xi = -a(j1 + 1) 1368 yr = a(k1) 1369 yi = -a(k1 + 1) 1370 a(j1) = yr 1371 a(j1 + 1) = yi 1372 a(k1) = xr 1373 a(k1 + 1) = xi 1374 j1 = j1 - nh 1375 k1 = k1 - 2 1376 xr = a(j1) 1377 xi = -a(j1 + 1) 1378 yr = a(k1) 1379 yi = -a(k1 + 1) 1380 a(j1) = yr 1381 a(j1 + 1) = yi 1382 a(k1) = xr 1383 a(k1 + 1) = xi 1384 j1 = j1 - nm 1385 k1 = k1 - nm 1386 xr = a(j1) 1387 xi = -a(j1 + 1) 1388 yr = a(k1) 1389 yi = -a(k1 + 1) 1390 a(j1) = yr 1391 a(j1 + 1) = yi 1392 a(k1) = xr 1393 a(k1 + 1) = xi 1394 end do 1395 k1 = 4 * k + ip(m + k) 1396 j1 = k1 + 2 1397 k1 = k1 + nh 1398 a(j1 - 1) = -a(j1 - 1) 1399 xr = a(j1) 1400 xi = -a(j1 + 1) 1401 yr = a(k1) 1402 yi = -a(k1 + 1) 1403 a(j1) = yr 1404 a(j1 + 1) = yi 1405 a(k1) = xr 1406 a(k1 + 1) = xi 1407 a(k1 + 3) = -a(k1 + 3) 1408 j1 = j1 + nm 1409 k1 = k1 + nm 1410 a(j1 - 1) = -a(j1 - 1) 1411 xr = a(j1) 1412 xi = -a(j1 + 1) 1413 yr = a(k1) 1414 yi = -a(k1 + 1) 1415 a(j1) = yr 1416 a(j1 + 1) = yi 1417 a(k1) = xr 1418 a(k1 + 1) = xi 1419 a(k1 + 3) = -a(k1 + 3) 1420 end do 1421 end if 1422 end 1423! 1424 subroutine bitrv216(a) 1425 real*8 a(0 : 31), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i 1426 real*8 x5r, x5i, x7r, x7i, x8r, x8i, x10r, x10i 1427 real*8 x11r, x11i, x12r, x12i, x13r, x13i, x14r, x14i 1428 x1r = a(2) 1429 x1i = a(3) 1430 x2r = a(4) 1431 x2i = a(5) 1432 x3r = a(6) 1433 x3i = a(7) 1434 x4r = a(8) 1435 x4i = a(9) 1436 x5r = a(10) 1437 x5i = a(11) 1438 x7r = a(14) 1439 x7i = a(15) 1440 x8r = a(16) 1441 x8i = a(17) 1442 x10r = a(20) 1443 x10i = a(21) 1444 x11r = a(22) 1445 x11i = a(23) 1446 x12r = a(24) 1447 x12i = a(25) 1448 x13r = a(26) 1449 x13i = a(27) 1450 x14r = a(28) 1451 x14i = a(29) 1452 a(2) = x8r 1453 a(3) = x8i 1454 a(4) = x4r 1455 a(5) = x4i 1456 a(6) = x12r 1457 a(7) = x12i 1458 a(8) = x2r 1459 a(9) = x2i 1460 a(10) = x10r 1461 a(11) = x10i 1462 a(14) = x14r 1463 a(15) = x14i 1464 a(16) = x1r 1465 a(17) = x1i 1466 a(20) = x5r 1467 a(21) = x5i 1468 a(22) = x13r 1469 a(23) = x13i 1470 a(24) = x3r 1471 a(25) = x3i 1472 a(26) = x11r 1473 a(27) = x11i 1474 a(28) = x7r 1475 a(29) = x7i 1476 end 1477! 1478 subroutine bitrv216neg(a) 1479 real*8 a(0 : 31), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i 1480 real*8 x5r, x5i, x6r, x6i, x7r, x7i, x8r, x8i 1481 real*8 x9r, x9i, x10r, x10i, x11r, x11i, x12r, x12i 1482 real*8 x13r, x13i, x14r, x14i, x15r, x15i 1483 x1r = a(2) 1484 x1i = a(3) 1485 x2r = a(4) 1486 x2i = a(5) 1487 x3r = a(6) 1488 x3i = a(7) 1489 x4r = a(8) 1490 x4i = a(9) 1491 x5r = a(10) 1492 x5i = a(11) 1493 x6r = a(12) 1494 x6i = a(13) 1495 x7r = a(14) 1496 x7i = a(15) 1497 x8r = a(16) 1498 x8i = a(17) 1499 x9r = a(18) 1500 x9i = a(19) 1501 x10r = a(20) 1502 x10i = a(21) 1503 x11r = a(22) 1504 x11i = a(23) 1505 x12r = a(24) 1506 x12i = a(25) 1507 x13r = a(26) 1508 x13i = a(27) 1509 x14r = a(28) 1510 x14i = a(29) 1511 x15r = a(30) 1512 x15i = a(31) 1513 a(2) = x15r 1514 a(3) = x15i 1515 a(4) = x7r 1516 a(5) = x7i 1517 a(6) = x11r 1518 a(7) = x11i 1519 a(8) = x3r 1520 a(9) = x3i 1521 a(10) = x13r 1522 a(11) = x13i 1523 a(12) = x5r 1524 a(13) = x5i 1525 a(14) = x9r 1526 a(15) = x9i 1527 a(16) = x1r 1528 a(17) = x1i 1529 a(18) = x14r 1530 a(19) = x14i 1531 a(20) = x6r 1532 a(21) = x6i 1533 a(22) = x10r 1534 a(23) = x10i 1535 a(24) = x2r 1536 a(25) = x2i 1537 a(26) = x12r 1538 a(27) = x12i 1539 a(28) = x4r 1540 a(29) = x4i 1541 a(30) = x8r 1542 a(31) = x8i 1543 end 1544! 1545 subroutine bitrv208(a) 1546 real*8 a(0 : 15), x1r, x1i, x3r, x3i, x4r, x4i, x6r, x6i 1547 x1r = a(2) 1548 x1i = a(3) 1549 x3r = a(6) 1550 x3i = a(7) 1551 x4r = a(8) 1552 x4i = a(9) 1553 x6r = a(12) 1554 x6i = a(13) 1555 a(2) = x4r 1556 a(3) = x4i 1557 a(6) = x6r 1558 a(7) = x6i 1559 a(8) = x1r 1560 a(9) = x1i 1561 a(12) = x3r 1562 a(13) = x3i 1563 end 1564! 1565 subroutine bitrv208neg(a) 1566 real*8 a(0 : 15), x1r, x1i, x2r, x2i, x3r, x3i, x4r, x4i 1567 real*8 x5r, x5i, x6r, x6i, x7r, x7i 1568 x1r = a(2) 1569 x1i = a(3) 1570 x2r = a(4) 1571 x2i = a(5) 1572 x3r = a(6) 1573 x3i = a(7) 1574 x4r = a(8) 1575 x4i = a(9) 1576 x5r = a(10) 1577 x5i = a(11) 1578 x6r = a(12) 1579 x6i = a(13) 1580 x7r = a(14) 1581 x7i = a(15) 1582 a(2) = x7r 1583 a(3) = x7i 1584 a(4) = x3r 1585 a(5) = x3i 1586 a(6) = x5r 1587 a(7) = x5i 1588 a(8) = x1r 1589 a(9) = x1i 1590 a(10) = x6r 1591 a(11) = x6i 1592 a(12) = x2r 1593 a(13) = x2i 1594 a(14) = x4r 1595 a(15) = x4i 1596 end 1597! 1598 subroutine cftf1st(n, a, w) 1599 integer n, j, j0, j1, j2, j3, k, m, mh 1600 real*8 a(0 : n - 1), w(0 : *) 1601 real*8 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i 1602 real*8 wd1r, wd1i, wd3r, wd3i 1603 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 1604 real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i 1605 mh = n / 8 1606 m = 2 * mh 1607 j1 = m 1608 j2 = j1 + m 1609 j3 = j2 + m 1610 x0r = a(0) + a(j2) 1611 x0i = a(1) + a(j2 + 1) 1612 x1r = a(0) - a(j2) 1613 x1i = a(1) - a(j2 + 1) 1614 x2r = a(j1) + a(j3) 1615 x2i = a(j1 + 1) + a(j3 + 1) 1616 x3r = a(j1) - a(j3) 1617 x3i = a(j1 + 1) - a(j3 + 1) 1618 a(0) = x0r + x2r 1619 a(1) = x0i + x2i 1620 a(j1) = x0r - x2r 1621 a(j1 + 1) = x0i - x2i 1622 a(j2) = x1r - x3i 1623 a(j2 + 1) = x1i + x3r 1624 a(j3) = x1r + x3i 1625 a(j3 + 1) = x1i - x3r 1626 wn4r = w(1) 1627 csc1 = w(2) 1628 csc3 = w(3) 1629 wd1r = 1 1630 wd1i = 0 1631 wd3r = 1 1632 wd3i = 0 1633 k = 0 1634 do j = 2, mh - 6, 4 1635 k = k + 4 1636 wk1r = csc1 * (wd1r + w(k)) 1637 wk1i = csc1 * (wd1i + w(k + 1)) 1638 wk3r = csc3 * (wd3r + w(k + 2)) 1639 wk3i = csc3 * (wd3i + w(k + 3)) 1640 wd1r = w(k) 1641 wd1i = w(k + 1) 1642 wd3r = w(k + 2) 1643 wd3i = w(k + 3) 1644 j1 = j + m 1645 j2 = j1 + m 1646 j3 = j2 + m 1647 x0r = a(j) + a(j2) 1648 x0i = a(j + 1) + a(j2 + 1) 1649 x1r = a(j) - a(j2) 1650 x1i = a(j + 1) - a(j2 + 1) 1651 y0r = a(j + 2) + a(j2 + 2) 1652 y0i = a(j + 3) + a(j2 + 3) 1653 y1r = a(j + 2) - a(j2 + 2) 1654 y1i = a(j + 3) - a(j2 + 3) 1655 x2r = a(j1) + a(j3) 1656 x2i = a(j1 + 1) + a(j3 + 1) 1657 x3r = a(j1) - a(j3) 1658 x3i = a(j1 + 1) - a(j3 + 1) 1659 y2r = a(j1 + 2) + a(j3 + 2) 1660 y2i = a(j1 + 3) + a(j3 + 3) 1661 y3r = a(j1 + 2) - a(j3 + 2) 1662 y3i = a(j1 + 3) - a(j3 + 3) 1663 a(j) = x0r + x2r 1664 a(j + 1) = x0i + x2i 1665 a(j + 2) = y0r + y2r 1666 a(j + 3) = y0i + y2i 1667 a(j1) = x0r - x2r 1668 a(j1 + 1) = x0i - x2i 1669 a(j1 + 2) = y0r - y2r 1670 a(j1 + 3) = y0i - y2i 1671 x0r = x1r - x3i 1672 x0i = x1i + x3r 1673 a(j2) = wk1r * x0r - wk1i * x0i 1674 a(j2 + 1) = wk1r * x0i + wk1i * x0r 1675 x0r = y1r - y3i 1676 x0i = y1i + y3r 1677 a(j2 + 2) = wd1r * x0r - wd1i * x0i 1678 a(j2 + 3) = wd1r * x0i + wd1i * x0r 1679 x0r = x1r + x3i 1680 x0i = x1i - x3r 1681 a(j3) = wk3r * x0r + wk3i * x0i 1682 a(j3 + 1) = wk3r * x0i - wk3i * x0r 1683 x0r = y1r + y3i 1684 x0i = y1i - y3r 1685 a(j3 + 2) = wd3r * x0r + wd3i * x0i 1686 a(j3 + 3) = wd3r * x0i - wd3i * x0r 1687 j0 = m - j 1688 j1 = j0 + m 1689 j2 = j1 + m 1690 j3 = j2 + m 1691 x0r = a(j0) + a(j2) 1692 x0i = a(j0 + 1) + a(j2 + 1) 1693 x1r = a(j0) - a(j2) 1694 x1i = a(j0 + 1) - a(j2 + 1) 1695 y0r = a(j0 - 2) + a(j2 - 2) 1696 y0i = a(j0 - 1) + a(j2 - 1) 1697 y1r = a(j0 - 2) - a(j2 - 2) 1698 y1i = a(j0 - 1) - a(j2 - 1) 1699 x2r = a(j1) + a(j3) 1700 x2i = a(j1 + 1) + a(j3 + 1) 1701 x3r = a(j1) - a(j3) 1702 x3i = a(j1 + 1) - a(j3 + 1) 1703 y2r = a(j1 - 2) + a(j3 - 2) 1704 y2i = a(j1 - 1) + a(j3 - 1) 1705 y3r = a(j1 - 2) - a(j3 - 2) 1706 y3i = a(j1 - 1) - a(j3 - 1) 1707 a(j0) = x0r + x2r 1708 a(j0 + 1) = x0i + x2i 1709 a(j0 - 2) = y0r + y2r 1710 a(j0 - 1) = y0i + y2i 1711 a(j1) = x0r - x2r 1712 a(j1 + 1) = x0i - x2i 1713 a(j1 - 2) = y0r - y2r 1714 a(j1 - 1) = y0i - y2i 1715 x0r = x1r - x3i 1716 x0i = x1i + x3r 1717 a(j2) = wk1i * x0r - wk1r * x0i 1718 a(j2 + 1) = wk1i * x0i + wk1r * x0r 1719 x0r = y1r - y3i 1720 x0i = y1i + y3r 1721 a(j2 - 2) = wd1i * x0r - wd1r * x0i 1722 a(j2 - 1) = wd1i * x0i + wd1r * x0r 1723 x0r = x1r + x3i 1724 x0i = x1i - x3r 1725 a(j3) = wk3i * x0r + wk3r * x0i 1726 a(j3 + 1) = wk3i * x0i - wk3r * x0r 1727 x0r = y1r + y3i 1728 x0i = y1i - y3r 1729 a(j3 - 2) = wd3i * x0r + wd3r * x0i 1730 a(j3 - 1) = wd3i * x0i - wd3r * x0r 1731 end do 1732 wk1r = csc1 * (wd1r + wn4r) 1733 wk1i = csc1 * (wd1i + wn4r) 1734 wk3r = csc3 * (wd3r - wn4r) 1735 wk3i = csc3 * (wd3i - wn4r) 1736 j0 = mh 1737 j1 = j0 + m 1738 j2 = j1 + m 1739 j3 = j2 + m 1740 x0r = a(j0 - 2) + a(j2 - 2) 1741 x0i = a(j0 - 1) + a(j2 - 1) 1742 x1r = a(j0 - 2) - a(j2 - 2) 1743 x1i = a(j0 - 1) - a(j2 - 1) 1744 x2r = a(j1 - 2) + a(j3 - 2) 1745 x2i = a(j1 - 1) + a(j3 - 1) 1746 x3r = a(j1 - 2) - a(j3 - 2) 1747 x3i = a(j1 - 1) - a(j3 - 1) 1748 a(j0 - 2) = x0r + x2r 1749 a(j0 - 1) = x0i + x2i 1750 a(j1 - 2) = x0r - x2r 1751 a(j1 - 1) = x0i - x2i 1752 x0r = x1r - x3i 1753 x0i = x1i + x3r 1754 a(j2 - 2) = wk1r * x0r - wk1i * x0i 1755 a(j2 - 1) = wk1r * x0i + wk1i * x0r 1756 x0r = x1r + x3i 1757 x0i = x1i - x3r 1758 a(j3 - 2) = wk3r * x0r + wk3i * x0i 1759 a(j3 - 1) = wk3r * x0i - wk3i * x0r 1760 x0r = a(j0) + a(j2) 1761 x0i = a(j0 + 1) + a(j2 + 1) 1762 x1r = a(j0) - a(j2) 1763 x1i = a(j0 + 1) - a(j2 + 1) 1764 x2r = a(j1) + a(j3) 1765 x2i = a(j1 + 1) + a(j3 + 1) 1766 x3r = a(j1) - a(j3) 1767 x3i = a(j1 + 1) - a(j3 + 1) 1768 a(j0) = x0r + x2r 1769 a(j0 + 1) = x0i + x2i 1770 a(j1) = x0r - x2r 1771 a(j1 + 1) = x0i - x2i 1772 x0r = x1r - x3i 1773 x0i = x1i + x3r 1774 a(j2) = wn4r * (x0r - x0i) 1775 a(j2 + 1) = wn4r * (x0i + x0r) 1776 x0r = x1r + x3i 1777 x0i = x1i - x3r 1778 a(j3) = -wn4r * (x0r + x0i) 1779 a(j3 + 1) = -wn4r * (x0i - x0r) 1780 x0r = a(j0 + 2) + a(j2 + 2) 1781 x0i = a(j0 + 3) + a(j2 + 3) 1782 x1r = a(j0 + 2) - a(j2 + 2) 1783 x1i = a(j0 + 3) - a(j2 + 3) 1784 x2r = a(j1 + 2) + a(j3 + 2) 1785 x2i = a(j1 + 3) + a(j3 + 3) 1786 x3r = a(j1 + 2) - a(j3 + 2) 1787 x3i = a(j1 + 3) - a(j3 + 3) 1788 a(j0 + 2) = x0r + x2r 1789 a(j0 + 3) = x0i + x2i 1790 a(j1 + 2) = x0r - x2r 1791 a(j1 + 3) = x0i - x2i 1792 x0r = x1r - x3i 1793 x0i = x1i + x3r 1794 a(j2 + 2) = wk1i * x0r - wk1r * x0i 1795 a(j2 + 3) = wk1i * x0i + wk1r * x0r 1796 x0r = x1r + x3i 1797 x0i = x1i - x3r 1798 a(j3 + 2) = wk3i * x0r + wk3r * x0i 1799 a(j3 + 3) = wk3i * x0i - wk3r * x0r 1800 end 1801! 1802 subroutine cftb1st(n, a, w) 1803 integer n, j, j0, j1, j2, j3, k, m, mh 1804 real*8 a(0 : n - 1), w(0 : *) 1805 real*8 wn4r, csc1, csc3, wk1r, wk1i, wk3r, wk3i 1806 real*8 wd1r, wd1i, wd3r, wd3i 1807 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 1808 real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i 1809 mh = n / 8 1810 m = 2 * mh 1811 j1 = m 1812 j2 = j1 + m 1813 j3 = j2 + m 1814 x0r = a(0) + a(j2) 1815 x0i = -a(1) - a(j2 + 1) 1816 x1r = a(0) - a(j2) 1817 x1i = -a(1) + a(j2 + 1) 1818 x2r = a(j1) + a(j3) 1819 x2i = a(j1 + 1) + a(j3 + 1) 1820 x3r = a(j1) - a(j3) 1821 x3i = a(j1 + 1) - a(j3 + 1) 1822 a(0) = x0r + x2r 1823 a(1) = x0i - x2i 1824 a(j1) = x0r - x2r 1825 a(j1 + 1) = x0i + x2i 1826 a(j2) = x1r + x3i 1827 a(j2 + 1) = x1i + x3r 1828 a(j3) = x1r - x3i 1829 a(j3 + 1) = x1i - x3r 1830 wn4r = w(1) 1831 csc1 = w(2) 1832 csc3 = w(3) 1833 wd1r = 1 1834 wd1i = 0 1835 wd3r = 1 1836 wd3i = 0 1837 k = 0 1838 do j = 2, mh - 6, 4 1839 k = k + 4 1840 wk1r = csc1 * (wd1r + w(k)) 1841 wk1i = csc1 * (wd1i + w(k + 1)) 1842 wk3r = csc3 * (wd3r + w(k + 2)) 1843 wk3i = csc3 * (wd3i + w(k + 3)) 1844 wd1r = w(k) 1845 wd1i = w(k + 1) 1846 wd3r = w(k + 2) 1847 wd3i = w(k + 3) 1848 j1 = j + m 1849 j2 = j1 + m 1850 j3 = j2 + m 1851 x0r = a(j) + a(j2) 1852 x0i = -a(j + 1) - a(j2 + 1) 1853 x1r = a(j) - a(j2) 1854 x1i = -a(j + 1) + a(j2 + 1) 1855 y0r = a(j + 2) + a(j2 + 2) 1856 y0i = -a(j + 3) - a(j2 + 3) 1857 y1r = a(j + 2) - a(j2 + 2) 1858 y1i = -a(j + 3) + a(j2 + 3) 1859 x2r = a(j1) + a(j3) 1860 x2i = a(j1 + 1) + a(j3 + 1) 1861 x3r = a(j1) - a(j3) 1862 x3i = a(j1 + 1) - a(j3 + 1) 1863 y2r = a(j1 + 2) + a(j3 + 2) 1864 y2i = a(j1 + 3) + a(j3 + 3) 1865 y3r = a(j1 + 2) - a(j3 + 2) 1866 y3i = a(j1 + 3) - a(j3 + 3) 1867 a(j) = x0r + x2r 1868 a(j + 1) = x0i - x2i 1869 a(j + 2) = y0r + y2r 1870 a(j + 3) = y0i - y2i 1871 a(j1) = x0r - x2r 1872 a(j1 + 1) = x0i + x2i 1873 a(j1 + 2) = y0r - y2r 1874 a(j1 + 3) = y0i + y2i 1875 x0r = x1r + x3i 1876 x0i = x1i + x3r 1877 a(j2) = wk1r * x0r - wk1i * x0i 1878 a(j2 + 1) = wk1r * x0i + wk1i * x0r 1879 x0r = y1r + y3i 1880 x0i = y1i + y3r 1881 a(j2 + 2) = wd1r * x0r - wd1i * x0i 1882 a(j2 + 3) = wd1r * x0i + wd1i * x0r 1883 x0r = x1r - x3i 1884 x0i = x1i - x3r 1885 a(j3) = wk3r * x0r + wk3i * x0i 1886 a(j3 + 1) = wk3r * x0i - wk3i * x0r 1887 x0r = y1r - y3i 1888 x0i = y1i - y3r 1889 a(j3 + 2) = wd3r * x0r + wd3i * x0i 1890 a(j3 + 3) = wd3r * x0i - wd3i * x0r 1891 j0 = m - j 1892 j1 = j0 + m 1893 j2 = j1 + m 1894 j3 = j2 + m 1895 x0r = a(j0) + a(j2) 1896 x0i = -a(j0 + 1) - a(j2 + 1) 1897 x1r = a(j0) - a(j2) 1898 x1i = -a(j0 + 1) + a(j2 + 1) 1899 y0r = a(j0 - 2) + a(j2 - 2) 1900 y0i = -a(j0 - 1) - a(j2 - 1) 1901 y1r = a(j0 - 2) - a(j2 - 2) 1902 y1i = -a(j0 - 1) + a(j2 - 1) 1903 x2r = a(j1) + a(j3) 1904 x2i = a(j1 + 1) + a(j3 + 1) 1905 x3r = a(j1) - a(j3) 1906 x3i = a(j1 + 1) - a(j3 + 1) 1907 y2r = a(j1 - 2) + a(j3 - 2) 1908 y2i = a(j1 - 1) + a(j3 - 1) 1909 y3r = a(j1 - 2) - a(j3 - 2) 1910 y3i = a(j1 - 1) - a(j3 - 1) 1911 a(j0) = x0r + x2r 1912 a(j0 + 1) = x0i - x2i 1913 a(j0 - 2) = y0r + y2r 1914 a(j0 - 1) = y0i - y2i 1915 a(j1) = x0r - x2r 1916 a(j1 + 1) = x0i + x2i 1917 a(j1 - 2) = y0r - y2r 1918 a(j1 - 1) = y0i + y2i 1919 x0r = x1r + x3i 1920 x0i = x1i + x3r 1921 a(j2) = wk1i * x0r - wk1r * x0i 1922 a(j2 + 1) = wk1i * x0i + wk1r * x0r 1923 x0r = y1r + y3i 1924 x0i = y1i + y3r 1925 a(j2 - 2) = wd1i * x0r - wd1r * x0i 1926 a(j2 - 1) = wd1i * x0i + wd1r * x0r 1927 x0r = x1r - x3i 1928 x0i = x1i - x3r 1929 a(j3) = wk3i * x0r + wk3r * x0i 1930 a(j3 + 1) = wk3i * x0i - wk3r * x0r 1931 x0r = y1r - y3i 1932 x0i = y1i - y3r 1933 a(j3 - 2) = wd3i * x0r + wd3r * x0i 1934 a(j3 - 1) = wd3i * x0i - wd3r * x0r 1935 end do 1936 wk1r = csc1 * (wd1r + wn4r) 1937 wk1i = csc1 * (wd1i + wn4r) 1938 wk3r = csc3 * (wd3r - wn4r) 1939 wk3i = csc3 * (wd3i - wn4r) 1940 j0 = mh 1941 j1 = j0 + m 1942 j2 = j1 + m 1943 j3 = j2 + m 1944 x0r = a(j0 - 2) + a(j2 - 2) 1945 x0i = -a(j0 - 1) - a(j2 - 1) 1946 x1r = a(j0 - 2) - a(j2 - 2) 1947 x1i = -a(j0 - 1) + a(j2 - 1) 1948 x2r = a(j1 - 2) + a(j3 - 2) 1949 x2i = a(j1 - 1) + a(j3 - 1) 1950 x3r = a(j1 - 2) - a(j3 - 2) 1951 x3i = a(j1 - 1) - a(j3 - 1) 1952 a(j0 - 2) = x0r + x2r 1953 a(j0 - 1) = x0i - x2i 1954 a(j1 - 2) = x0r - x2r 1955 a(j1 - 1) = x0i + x2i 1956 x0r = x1r + x3i 1957 x0i = x1i + x3r 1958 a(j2 - 2) = wk1r * x0r - wk1i * x0i 1959 a(j2 - 1) = wk1r * x0i + wk1i * x0r 1960 x0r = x1r - x3i 1961 x0i = x1i - x3r 1962 a(j3 - 2) = wk3r * x0r + wk3i * x0i 1963 a(j3 - 1) = wk3r * x0i - wk3i * x0r 1964 x0r = a(j0) + a(j2) 1965 x0i = -a(j0 + 1) - a(j2 + 1) 1966 x1r = a(j0) - a(j2) 1967 x1i = -a(j0 + 1) + a(j2 + 1) 1968 x2r = a(j1) + a(j3) 1969 x2i = a(j1 + 1) + a(j3 + 1) 1970 x3r = a(j1) - a(j3) 1971 x3i = a(j1 + 1) - a(j3 + 1) 1972 a(j0) = x0r + x2r 1973 a(j0 + 1) = x0i - x2i 1974 a(j1) = x0r - x2r 1975 a(j1 + 1) = x0i + x2i 1976 x0r = x1r + x3i 1977 x0i = x1i + x3r 1978 a(j2) = wn4r * (x0r - x0i) 1979 a(j2 + 1) = wn4r * (x0i + x0r) 1980 x0r = x1r - x3i 1981 x0i = x1i - x3r 1982 a(j3) = -wn4r * (x0r + x0i) 1983 a(j3 + 1) = -wn4r * (x0i - x0r) 1984 x0r = a(j0 + 2) + a(j2 + 2) 1985 x0i = -a(j0 + 3) - a(j2 + 3) 1986 x1r = a(j0 + 2) - a(j2 + 2) 1987 x1i = -a(j0 + 3) + a(j2 + 3) 1988 x2r = a(j1 + 2) + a(j3 + 2) 1989 x2i = a(j1 + 3) + a(j3 + 3) 1990 x3r = a(j1 + 2) - a(j3 + 2) 1991 x3i = a(j1 + 3) - a(j3 + 3) 1992 a(j0 + 2) = x0r + x2r 1993 a(j0 + 3) = x0i - x2i 1994 a(j1 + 2) = x0r - x2r 1995 a(j1 + 3) = x0i + x2i 1996 x0r = x1r + x3i 1997 x0i = x1i + x3r 1998 a(j2 + 2) = wk1i * x0r - wk1r * x0i 1999 a(j2 + 3) = wk1i * x0i + wk1r * x0r 2000 x0r = x1r - x3i 2001 x0i = x1i - x3r 2002 a(j3 + 2) = wk3i * x0r + wk3r * x0i 2003 a(j3 + 3) = wk3i * x0i - wk3r * x0r 2004 end 2005! 2006 subroutine cftrec4(n, a, nw, w) 2007 integer n, nw, cfttree, isplt, j, k, m 2008 real*8 a(0 : n - 1), w(0 : nw - 1) 2009 m = n 2010 do while (m .gt. 512) 2011 m = m / 4 2012 call cftmdl1(m, a(n - m), w(nw - m / 2)) 2013 end do 2014 call cftleaf(m, 1, a(n - m), nw, w) 2015 k = 0 2016 do j = n - m, m, -m 2017 k = k + 1 2018 isplt = cfttree(m, j, k, a, nw, w) 2019 call cftleaf(m, isplt, a(j - m), nw, w) 2020 end do 2021 end 2022! 2023 integer function cfttree(n, j, k, a, nw, w) 2024 integer n, j, k, nw, i, isplt, m 2025 real*8 a(0 : j - 1), w(0 : nw - 1) 2026 if (mod(k, 4) .ne. 0) then 2027 isplt = mod(k, 2) 2028 if (isplt .ne. 0) then 2029 call cftmdl1(n, a(j - n), w(nw - n / 2)) 2030 else 2031 call cftmdl2(n, a(j - n), w(nw - n)) 2032 end if 2033 else 2034 m = n 2035 i = k 2036 do while (mod(i, 4) .eq. 0) 2037 m = m * 4 2038 i = i / 4 2039 end do 2040 isplt = mod(i, 2) 2041 if (isplt .ne. 0) then 2042 do while (m .gt. 128) 2043 call cftmdl1(m, a(j - m), w(nw - m / 2)) 2044 m = m / 4 2045 end do 2046 else 2047 do while (m .gt. 128) 2048 call cftmdl2(m, a(j - m), w(nw - m)) 2049 m = m / 4 2050 end do 2051 end if 2052 end if 2053 cfttree = isplt 2054 end 2055! 2056 subroutine cftleaf(n, isplt, a, nw, w) 2057 integer n, isplt, nw 2058 real*8 a(0 : n - 1), w(0 : nw - 1) 2059 if (n .eq. 512) then 2060 call cftmdl1(128, a, w(nw - 64)) 2061 call cftf161(a, w(nw - 8)) 2062 call cftf162(a(32), w(nw - 32)) 2063 call cftf161(a(64), w(nw - 8)) 2064 call cftf161(a(96), w(nw - 8)) 2065 call cftmdl2(128, a(128), w(nw - 128)) 2066 call cftf161(a(128), w(nw - 8)) 2067 call cftf162(a(160), w(nw - 32)) 2068 call cftf161(a(192), w(nw - 8)) 2069 call cftf162(a(224), w(nw - 32)) 2070 call cftmdl1(128, a(256), w(nw - 64)) 2071 call cftf161(a(256), w(nw - 8)) 2072 call cftf162(a(288), w(nw - 32)) 2073 call cftf161(a(320), w(nw - 8)) 2074 call cftf161(a(352), w(nw - 8)) 2075 if (isplt .ne. 0) then 2076 call cftmdl1(128, a(384), w(nw - 64)) 2077 call cftf161(a(480), w(nw - 8)) 2078 else 2079 call cftmdl2(128, a(384), w(nw - 128)) 2080 call cftf162(a(480), w(nw - 32)) 2081 end if 2082 call cftf161(a(384), w(nw - 8)) 2083 call cftf162(a(416), w(nw - 32)) 2084 call cftf161(a(448), w(nw - 8)) 2085 else 2086 call cftmdl1(64, a, w(nw - 32)) 2087 call cftf081(a, w(nw - 8)) 2088 call cftf082(a(16), w(nw - 8)) 2089 call cftf081(a(32), w(nw - 8)) 2090 call cftf081(a(48), w(nw - 8)) 2091 call cftmdl2(64, a(64), w(nw - 64)) 2092 call cftf081(a(64), w(nw - 8)) 2093 call cftf082(a(80), w(nw - 8)) 2094 call cftf081(a(96), w(nw - 8)) 2095 call cftf082(a(112), w(nw - 8)) 2096 call cftmdl1(64, a(128), w(nw - 32)) 2097 call cftf081(a(128), w(nw - 8)) 2098 call cftf082(a(144), w(nw - 8)) 2099 call cftf081(a(160), w(nw - 8)) 2100 call cftf081(a(176), w(nw - 8)) 2101 if (isplt .ne. 0) then 2102 call cftmdl1(64, a(192), w(nw - 32)) 2103 call cftf081(a(240), w(nw - 8)) 2104 else 2105 call cftmdl2(64, a(192), w(nw - 64)) 2106 call cftf082(a(240), w(nw - 8)) 2107 end if 2108 call cftf081(a(192), w(nw - 8)) 2109 call cftf082(a(208), w(nw - 8)) 2110 call cftf081(a(224), w(nw - 8)) 2111 end if 2112 end 2113! 2114 subroutine cftmdl1(n, a, w) 2115 integer n, j, j0, j1, j2, j3, k, m, mh 2116 real*8 a(0 : n - 1), w(0 : *) 2117 real*8 wn4r, wk1r, wk1i, wk3r, wk3i 2118 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 2119 mh = n / 8 2120 m = 2 * mh 2121 j1 = m 2122 j2 = j1 + m 2123 j3 = j2 + m 2124 x0r = a(0) + a(j2) 2125 x0i = a(1) + a(j2 + 1) 2126 x1r = a(0) - a(j2) 2127 x1i = a(1) - a(j2 + 1) 2128 x2r = a(j1) + a(j3) 2129 x2i = a(j1 + 1) + a(j3 + 1) 2130 x3r = a(j1) - a(j3) 2131 x3i = a(j1 + 1) - a(j3 + 1) 2132 a(0) = x0r + x2r 2133 a(1) = x0i + x2i 2134 a(j1) = x0r - x2r 2135 a(j1 + 1) = x0i - x2i 2136 a(j2) = x1r - x3i 2137 a(j2 + 1) = x1i + x3r 2138 a(j3) = x1r + x3i 2139 a(j3 + 1) = x1i - x3r 2140 wn4r = w(1) 2141 k = 0 2142 do j = 2, mh - 2, 2 2143 k = k + 4 2144 wk1r = w(k) 2145 wk1i = w(k + 1) 2146 wk3r = w(k + 2) 2147 wk3i = w(k + 3) 2148 j1 = j + m 2149 j2 = j1 + m 2150 j3 = j2 + m 2151 x0r = a(j) + a(j2) 2152 x0i = a(j + 1) + a(j2 + 1) 2153 x1r = a(j) - a(j2) 2154 x1i = a(j + 1) - a(j2 + 1) 2155 x2r = a(j1) + a(j3) 2156 x2i = a(j1 + 1) + a(j3 + 1) 2157 x3r = a(j1) - a(j3) 2158 x3i = a(j1 + 1) - a(j3 + 1) 2159 a(j) = x0r + x2r 2160 a(j + 1) = x0i + x2i 2161 a(j1) = x0r - x2r 2162 a(j1 + 1) = x0i - x2i 2163 x0r = x1r - x3i 2164 x0i = x1i + x3r 2165 a(j2) = wk1r * x0r - wk1i * x0i 2166 a(j2 + 1) = wk1r * x0i + wk1i * x0r 2167 x0r = x1r + x3i 2168 x0i = x1i - x3r 2169 a(j3) = wk3r * x0r + wk3i * x0i 2170 a(j3 + 1) = wk3r * x0i - wk3i * x0r 2171 j0 = m - j 2172 j1 = j0 + m 2173 j2 = j1 + m 2174 j3 = j2 + m 2175 x0r = a(j0) + a(j2) 2176 x0i = a(j0 + 1) + a(j2 + 1) 2177 x1r = a(j0) - a(j2) 2178 x1i = a(j0 + 1) - a(j2 + 1) 2179 x2r = a(j1) + a(j3) 2180 x2i = a(j1 + 1) + a(j3 + 1) 2181 x3r = a(j1) - a(j3) 2182 x3i = a(j1 + 1) - a(j3 + 1) 2183 a(j0) = x0r + x2r 2184 a(j0 + 1) = x0i + x2i 2185 a(j1) = x0r - x2r 2186 a(j1 + 1) = x0i - x2i 2187 x0r = x1r - x3i 2188 x0i = x1i + x3r 2189 a(j2) = wk1i * x0r - wk1r * x0i 2190 a(j2 + 1) = wk1i * x0i + wk1r * x0r 2191 x0r = x1r + x3i 2192 x0i = x1i - x3r 2193 a(j3) = wk3i * x0r + wk3r * x0i 2194 a(j3 + 1) = wk3i * x0i - wk3r * x0r 2195 end do 2196 j0 = mh 2197 j1 = j0 + m 2198 j2 = j1 + m 2199 j3 = j2 + m 2200 x0r = a(j0) + a(j2) 2201 x0i = a(j0 + 1) + a(j2 + 1) 2202 x1r = a(j0) - a(j2) 2203 x1i = a(j0 + 1) - a(j2 + 1) 2204 x2r = a(j1) + a(j3) 2205 x2i = a(j1 + 1) + a(j3 + 1) 2206 x3r = a(j1) - a(j3) 2207 x3i = a(j1 + 1) - a(j3 + 1) 2208 a(j0) = x0r + x2r 2209 a(j0 + 1) = x0i + x2i 2210 a(j1) = x0r - x2r 2211 a(j1 + 1) = x0i - x2i 2212 x0r = x1r - x3i 2213 x0i = x1i + x3r 2214 a(j2) = wn4r * (x0r - x0i) 2215 a(j2 + 1) = wn4r * (x0i + x0r) 2216 x0r = x1r + x3i 2217 x0i = x1i - x3r 2218 a(j3) = -wn4r * (x0r + x0i) 2219 a(j3 + 1) = -wn4r * (x0i - x0r) 2220 end 2221! 2222 subroutine cftmdl2(n, a, w) 2223 integer n, j, j0, j1, j2, j3, k, kr, m, mh 2224 real*8 a(0 : n - 1), w(0 : *) 2225 real*8 wn4r, wk1r, wk1i, wk3r, wk3i, wd1r, wd1i, wd3r, wd3i 2226 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 2227 real*8 y0r, y0i, y2r, y2i 2228 mh = n / 8 2229 m = 2 * mh 2230 wn4r = w(1) 2231 j1 = m 2232 j2 = j1 + m 2233 j3 = j2 + m 2234 x0r = a(0) - a(j2 + 1) 2235 x0i = a(1) + a(j2) 2236 x1r = a(0) + a(j2 + 1) 2237 x1i = a(1) - a(j2) 2238 x2r = a(j1) - a(j3 + 1) 2239 x2i = a(j1 + 1) + a(j3) 2240 x3r = a(j1) + a(j3 + 1) 2241 x3i = a(j1 + 1) - a(j3) 2242 y0r = wn4r * (x2r - x2i) 2243 y0i = wn4r * (x2i + x2r) 2244 a(0) = x0r + y0r 2245 a(1) = x0i + y0i 2246 a(j1) = x0r - y0r 2247 a(j1 + 1) = x0i - y0i 2248 y0r = wn4r * (x3r - x3i) 2249 y0i = wn4r * (x3i + x3r) 2250 a(j2) = x1r - y0i 2251 a(j2 + 1) = x1i + y0r 2252 a(j3) = x1r + y0i 2253 a(j3 + 1) = x1i - y0r 2254 k = 0 2255 kr = 2 * m 2256 do j = 2, mh - 2, 2 2257 k = k + 4 2258 wk1r = w(k) 2259 wk1i = w(k + 1) 2260 wk3r = w(k + 2) 2261 wk3i = w(k + 3) 2262 kr = kr - 4 2263 wd1i = w(kr) 2264 wd1r = w(kr + 1) 2265 wd3i = w(kr + 2) 2266 wd3r = w(kr + 3) 2267 j1 = j + m 2268 j2 = j1 + m 2269 j3 = j2 + m 2270 x0r = a(j) - a(j2 + 1) 2271 x0i = a(j + 1) + a(j2) 2272 x1r = a(j) + a(j2 + 1) 2273 x1i = a(j + 1) - a(j2) 2274 x2r = a(j1) - a(j3 + 1) 2275 x2i = a(j1 + 1) + a(j3) 2276 x3r = a(j1) + a(j3 + 1) 2277 x3i = a(j1 + 1) - a(j3) 2278 y0r = wk1r * x0r - wk1i * x0i 2279 y0i = wk1r * x0i + wk1i * x0r 2280 y2r = wd1r * x2r - wd1i * x2i 2281 y2i = wd1r * x2i + wd1i * x2r 2282 a(j) = y0r + y2r 2283 a(j + 1) = y0i + y2i 2284 a(j1) = y0r - y2r 2285 a(j1 + 1) = y0i - y2i 2286 y0r = wk3r * x1r + wk3i * x1i 2287 y0i = wk3r * x1i - wk3i * x1r 2288 y2r = wd3r * x3r + wd3i * x3i 2289 y2i = wd3r * x3i - wd3i * x3r 2290 a(j2) = y0r + y2r 2291 a(j2 + 1) = y0i + y2i 2292 a(j3) = y0r - y2r 2293 a(j3 + 1) = y0i - y2i 2294 j0 = m - j 2295 j1 = j0 + m 2296 j2 = j1 + m 2297 j3 = j2 + m 2298 x0r = a(j0) - a(j2 + 1) 2299 x0i = a(j0 + 1) + a(j2) 2300 x1r = a(j0) + a(j2 + 1) 2301 x1i = a(j0 + 1) - a(j2) 2302 x2r = a(j1) - a(j3 + 1) 2303 x2i = a(j1 + 1) + a(j3) 2304 x3r = a(j1) + a(j3 + 1) 2305 x3i = a(j1 + 1) - a(j3) 2306 y0r = wd1i * x0r - wd1r * x0i 2307 y0i = wd1i * x0i + wd1r * x0r 2308 y2r = wk1i * x2r - wk1r * x2i 2309 y2i = wk1i * x2i + wk1r * x2r 2310 a(j0) = y0r + y2r 2311 a(j0 + 1) = y0i + y2i 2312 a(j1) = y0r - y2r 2313 a(j1 + 1) = y0i - y2i 2314 y0r = wd3i * x1r + wd3r * x1i 2315 y0i = wd3i * x1i - wd3r * x1r 2316 y2r = wk3i * x3r + wk3r * x3i 2317 y2i = wk3i * x3i - wk3r * x3r 2318 a(j2) = y0r + y2r 2319 a(j2 + 1) = y0i + y2i 2320 a(j3) = y0r - y2r 2321 a(j3 + 1) = y0i - y2i 2322 end do 2323 wk1r = w(m) 2324 wk1i = w(m + 1) 2325 j0 = mh 2326 j1 = j0 + m 2327 j2 = j1 + m 2328 j3 = j2 + m 2329 x0r = a(j0) - a(j2 + 1) 2330 x0i = a(j0 + 1) + a(j2) 2331 x1r = a(j0) + a(j2 + 1) 2332 x1i = a(j0 + 1) - a(j2) 2333 x2r = a(j1) - a(j3 + 1) 2334 x2i = a(j1 + 1) + a(j3) 2335 x3r = a(j1) + a(j3 + 1) 2336 x3i = a(j1 + 1) - a(j3) 2337 y0r = wk1r * x0r - wk1i * x0i 2338 y0i = wk1r * x0i + wk1i * x0r 2339 y2r = wk1i * x2r - wk1r * x2i 2340 y2i = wk1i * x2i + wk1r * x2r 2341 a(j0) = y0r + y2r 2342 a(j0 + 1) = y0i + y2i 2343 a(j1) = y0r - y2r 2344 a(j1 + 1) = y0i - y2i 2345 y0r = wk1i * x1r - wk1r * x1i 2346 y0i = wk1i * x1i + wk1r * x1r 2347 y2r = wk1r * x3r - wk1i * x3i 2348 y2i = wk1r * x3i + wk1i * x3r 2349 a(j2) = y0r - y2r 2350 a(j2 + 1) = y0i - y2i 2351 a(j3) = y0r + y2r 2352 a(j3 + 1) = y0i + y2i 2353 end 2354! 2355 subroutine cftfx41(n, a, nw, w) 2356 integer n, nw 2357 real*8 a(0 : n - 1), w(0 : nw - 1) 2358 if (n .eq. 128) then 2359 call cftf161(a, w(nw - 8)) 2360 call cftf162(a(32), w(nw - 32)) 2361 call cftf161(a(64), w(nw - 8)) 2362 call cftf161(a(96), w(nw - 8)) 2363 else 2364 call cftf081(a, w(nw - 8)) 2365 call cftf082(a(16), w(nw - 8)) 2366 call cftf081(a(32), w(nw - 8)) 2367 call cftf081(a(48), w(nw - 8)) 2368 end if 2369 end 2370! 2371 subroutine cftf161(a, w) 2372 real*8 a(0 : 31), w(0 : *), wn4r, wk1r, wk1i 2373 real*8 x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 2374 real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i 2375 real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i 2376 real*8 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i 2377 real*8 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i 2378 wn4r = w(1) 2379 wk1r = w(2) 2380 wk1i = w(3) 2381 x0r = a(0) + a(16) 2382 x0i = a(1) + a(17) 2383 x1r = a(0) - a(16) 2384 x1i = a(1) - a(17) 2385 x2r = a(8) + a(24) 2386 x2i = a(9) + a(25) 2387 x3r = a(8) - a(24) 2388 x3i = a(9) - a(25) 2389 y0r = x0r + x2r 2390 y0i = x0i + x2i 2391 y4r = x0r - x2r 2392 y4i = x0i - x2i 2393 y8r = x1r - x3i 2394 y8i = x1i + x3r 2395 y12r = x1r + x3i 2396 y12i = x1i - x3r 2397 x0r = a(2) + a(18) 2398 x0i = a(3) + a(19) 2399 x1r = a(2) - a(18) 2400 x1i = a(3) - a(19) 2401 x2r = a(10) + a(26) 2402 x2i = a(11) + a(27) 2403 x3r = a(10) - a(26) 2404 x3i = a(11) - a(27) 2405 y1r = x0r + x2r 2406 y1i = x0i + x2i 2407 y5r = x0r - x2r 2408 y5i = x0i - x2i 2409 x0r = x1r - x3i 2410 x0i = x1i + x3r 2411 y9r = wk1r * x0r - wk1i * x0i 2412 y9i = wk1r * x0i + wk1i * x0r 2413 x0r = x1r + x3i 2414 x0i = x1i - x3r 2415 y13r = wk1i * x0r - wk1r * x0i 2416 y13i = wk1i * x0i + wk1r * x0r 2417 x0r = a(4) + a(20) 2418 x0i = a(5) + a(21) 2419 x1r = a(4) - a(20) 2420 x1i = a(5) - a(21) 2421 x2r = a(12) + a(28) 2422 x2i = a(13) + a(29) 2423 x3r = a(12) - a(28) 2424 x3i = a(13) - a(29) 2425 y2r = x0r + x2r 2426 y2i = x0i + x2i 2427 y6r = x0r - x2r 2428 y6i = x0i - x2i 2429 x0r = x1r - x3i 2430 x0i = x1i + x3r 2431 y10r = wn4r * (x0r - x0i) 2432 y10i = wn4r * (x0i + x0r) 2433 x0r = x1r + x3i 2434 x0i = x1i - x3r 2435 y14r = wn4r * (x0r + x0i) 2436 y14i = wn4r * (x0i - x0r) 2437 x0r = a(6) + a(22) 2438 x0i = a(7) + a(23) 2439 x1r = a(6) - a(22) 2440 x1i = a(7) - a(23) 2441 x2r = a(14) + a(30) 2442 x2i = a(15) + a(31) 2443 x3r = a(14) - a(30) 2444 x3i = a(15) - a(31) 2445 y3r = x0r + x2r 2446 y3i = x0i + x2i 2447 y7r = x0r - x2r 2448 y7i = x0i - x2i 2449 x0r = x1r - x3i 2450 x0i = x1i + x3r 2451 y11r = wk1i * x0r - wk1r * x0i 2452 y11i = wk1i * x0i + wk1r * x0r 2453 x0r = x1r + x3i 2454 x0i = x1i - x3r 2455 y15r = wk1r * x0r - wk1i * x0i 2456 y15i = wk1r * x0i + wk1i * x0r 2457 x0r = y12r - y14r 2458 x0i = y12i - y14i 2459 x1r = y12r + y14r 2460 x1i = y12i + y14i 2461 x2r = y13r - y15r 2462 x2i = y13i - y15i 2463 x3r = y13r + y15r 2464 x3i = y13i + y15i 2465 a(24) = x0r + x2r 2466 a(25) = x0i + x2i 2467 a(26) = x0r - x2r 2468 a(27) = x0i - x2i 2469 a(28) = x1r - x3i 2470 a(29) = x1i + x3r 2471 a(30) = x1r + x3i 2472 a(31) = x1i - x3r 2473 x0r = y8r + y10r 2474 x0i = y8i + y10i 2475 x1r = y8r - y10r 2476 x1i = y8i - y10i 2477 x2r = y9r + y11r 2478 x2i = y9i + y11i 2479 x3r = y9r - y11r 2480 x3i = y9i - y11i 2481 a(16) = x0r + x2r 2482 a(17) = x0i + x2i 2483 a(18) = x0r - x2r 2484 a(19) = x0i - x2i 2485 a(20) = x1r - x3i 2486 a(21) = x1i + x3r 2487 a(22) = x1r + x3i 2488 a(23) = x1i - x3r 2489 x0r = y5r - y7i 2490 x0i = y5i + y7r 2491 x2r = wn4r * (x0r - x0i) 2492 x2i = wn4r * (x0i + x0r) 2493 x0r = y5r + y7i 2494 x0i = y5i - y7r 2495 x3r = wn4r * (x0r - x0i) 2496 x3i = wn4r * (x0i + x0r) 2497 x0r = y4r - y6i 2498 x0i = y4i + y6r 2499 x1r = y4r + y6i 2500 x1i = y4i - y6r 2501 a(8) = x0r + x2r 2502 a(9) = x0i + x2i 2503 a(10) = x0r - x2r 2504 a(11) = x0i - x2i 2505 a(12) = x1r - x3i 2506 a(13) = x1i + x3r 2507 a(14) = x1r + x3i 2508 a(15) = x1i - x3r 2509 x0r = y0r + y2r 2510 x0i = y0i + y2i 2511 x1r = y0r - y2r 2512 x1i = y0i - y2i 2513 x2r = y1r + y3r 2514 x2i = y1i + y3i 2515 x3r = y1r - y3r 2516 x3i = y1i - y3i 2517 a(0) = x0r + x2r 2518 a(1) = x0i + x2i 2519 a(2) = x0r - x2r 2520 a(3) = x0i - x2i 2521 a(4) = x1r - x3i 2522 a(5) = x1i + x3r 2523 a(6) = x1r + x3i 2524 a(7) = x1i - x3r 2525 end 2526! 2527 subroutine cftf162(a, w) 2528 real*8 a(0 : 31), w(0 : *) 2529 real*8 wn4r, wk1r, wk1i, wk2r, wk2i, wk3r, wk3i 2530 real*8 x0r, x0i, x1r, x1i, x2r, x2i 2531 real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i 2532 real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i 2533 real*8 y8r, y8i, y9r, y9i, y10r, y10i, y11r, y11i 2534 real*8 y12r, y12i, y13r, y13i, y14r, y14i, y15r, y15i 2535 wn4r = w(1) 2536 wk1r = w(4) 2537 wk1i = w(5) 2538 wk3r = w(6) 2539 wk3i = -w(7) 2540 wk2r = w(8) 2541 wk2i = w(9) 2542 x1r = a(0) - a(17) 2543 x1i = a(1) + a(16) 2544 x0r = a(8) - a(25) 2545 x0i = a(9) + a(24) 2546 x2r = wn4r * (x0r - x0i) 2547 x2i = wn4r * (x0i + x0r) 2548 y0r = x1r + x2r 2549 y0i = x1i + x2i 2550 y4r = x1r - x2r 2551 y4i = x1i - x2i 2552 x1r = a(0) + a(17) 2553 x1i = a(1) - a(16) 2554 x0r = a(8) + a(25) 2555 x0i = a(9) - a(24) 2556 x2r = wn4r * (x0r - x0i) 2557 x2i = wn4r * (x0i + x0r) 2558 y8r = x1r - x2i 2559 y8i = x1i + x2r 2560 y12r = x1r + x2i 2561 y12i = x1i - x2r 2562 x0r = a(2) - a(19) 2563 x0i = a(3) + a(18) 2564 x1r = wk1r * x0r - wk1i * x0i 2565 x1i = wk1r * x0i + wk1i * x0r 2566 x0r = a(10) - a(27) 2567 x0i = a(11) + a(26) 2568 x2r = wk3i * x0r - wk3r * x0i 2569 x2i = wk3i * x0i + wk3r * x0r 2570 y1r = x1r + x2r 2571 y1i = x1i + x2i 2572 y5r = x1r - x2r 2573 y5i = x1i - x2i 2574 x0r = a(2) + a(19) 2575 x0i = a(3) - a(18) 2576 x1r = wk3r * x0r - wk3i * x0i 2577 x1i = wk3r * x0i + wk3i * x0r 2578 x0r = a(10) + a(27) 2579 x0i = a(11) - a(26) 2580 x2r = wk1r * x0r + wk1i * x0i 2581 x2i = wk1r * x0i - wk1i * x0r 2582 y9r = x1r - x2r 2583 y9i = x1i - x2i 2584 y13r = x1r + x2r 2585 y13i = x1i + x2i 2586 x0r = a(4) - a(21) 2587 x0i = a(5) + a(20) 2588 x1r = wk2r * x0r - wk2i * x0i 2589 x1i = wk2r * x0i + wk2i * x0r 2590 x0r = a(12) - a(29) 2591 x0i = a(13) + a(28) 2592 x2r = wk2i * x0r - wk2r * x0i 2593 x2i = wk2i * x0i + wk2r * x0r 2594 y2r = x1r + x2r 2595 y2i = x1i + x2i 2596 y6r = x1r - x2r 2597 y6i = x1i - x2i 2598 x0r = a(4) + a(21) 2599 x0i = a(5) - a(20) 2600 x1r = wk2i * x0r - wk2r * x0i 2601 x1i = wk2i * x0i + wk2r * x0r 2602 x0r = a(12) + a(29) 2603 x0i = a(13) - a(28) 2604 x2r = wk2r * x0r - wk2i * x0i 2605 x2i = wk2r * x0i + wk2i * x0r 2606 y10r = x1r - x2r 2607 y10i = x1i - x2i 2608 y14r = x1r + x2r 2609 y14i = x1i + x2i 2610 x0r = a(6) - a(23) 2611 x0i = a(7) + a(22) 2612 x1r = wk3r * x0r - wk3i * x0i 2613 x1i = wk3r * x0i + wk3i * x0r 2614 x0r = a(14) - a(31) 2615 x0i = a(15) + a(30) 2616 x2r = wk1i * x0r - wk1r * x0i 2617 x2i = wk1i * x0i + wk1r * x0r 2618 y3r = x1r + x2r 2619 y3i = x1i + x2i 2620 y7r = x1r - x2r 2621 y7i = x1i - x2i 2622 x0r = a(6) + a(23) 2623 x0i = a(7) - a(22) 2624 x1r = wk1i * x0r + wk1r * x0i 2625 x1i = wk1i * x0i - wk1r * x0r 2626 x0r = a(14) + a(31) 2627 x0i = a(15) - a(30) 2628 x2r = wk3i * x0r - wk3r * x0i 2629 x2i = wk3i * x0i + wk3r * x0r 2630 y11r = x1r + x2r 2631 y11i = x1i + x2i 2632 y15r = x1r - x2r 2633 y15i = x1i - x2i 2634 x1r = y0r + y2r 2635 x1i = y0i + y2i 2636 x2r = y1r + y3r 2637 x2i = y1i + y3i 2638 a(0) = x1r + x2r 2639 a(1) = x1i + x2i 2640 a(2) = x1r - x2r 2641 a(3) = x1i - x2i 2642 x1r = y0r - y2r 2643 x1i = y0i - y2i 2644 x2r = y1r - y3r 2645 x2i = y1i - y3i 2646 a(4) = x1r - x2i 2647 a(5) = x1i + x2r 2648 a(6) = x1r + x2i 2649 a(7) = x1i - x2r 2650 x1r = y4r - y6i 2651 x1i = y4i + y6r 2652 x0r = y5r - y7i 2653 x0i = y5i + y7r 2654 x2r = wn4r * (x0r - x0i) 2655 x2i = wn4r * (x0i + x0r) 2656 a(8) = x1r + x2r 2657 a(9) = x1i + x2i 2658 a(10) = x1r - x2r 2659 a(11) = x1i - x2i 2660 x1r = y4r + y6i 2661 x1i = y4i - y6r 2662 x0r = y5r + y7i 2663 x0i = y5i - y7r 2664 x2r = wn4r * (x0r - x0i) 2665 x2i = wn4r * (x0i + x0r) 2666 a(12) = x1r - x2i 2667 a(13) = x1i + x2r 2668 a(14) = x1r + x2i 2669 a(15) = x1i - x2r 2670 x1r = y8r + y10r 2671 x1i = y8i + y10i 2672 x2r = y9r - y11r 2673 x2i = y9i - y11i 2674 a(16) = x1r + x2r 2675 a(17) = x1i + x2i 2676 a(18) = x1r - x2r 2677 a(19) = x1i - x2i 2678 x1r = y8r - y10r 2679 x1i = y8i - y10i 2680 x2r = y9r + y11r 2681 x2i = y9i + y11i 2682 a(20) = x1r - x2i 2683 a(21) = x1i + x2r 2684 a(22) = x1r + x2i 2685 a(23) = x1i - x2r 2686 x1r = y12r - y14i 2687 x1i = y12i + y14r 2688 x0r = y13r + y15i 2689 x0i = y13i - y15r 2690 x2r = wn4r * (x0r - x0i) 2691 x2i = wn4r * (x0i + x0r) 2692 a(24) = x1r + x2r 2693 a(25) = x1i + x2i 2694 a(26) = x1r - x2r 2695 a(27) = x1i - x2i 2696 x1r = y12r + y14i 2697 x1i = y12i - y14r 2698 x0r = y13r - y15i 2699 x0i = y13i + y15r 2700 x2r = wn4r * (x0r - x0i) 2701 x2i = wn4r * (x0i + x0r) 2702 a(28) = x1r - x2i 2703 a(29) = x1i + x2r 2704 a(30) = x1r + x2i 2705 a(31) = x1i - x2r 2706 end 2707! 2708 subroutine cftf081(a, w) 2709 real*8 a(0 : 15), w(0 : *) 2710 real*8 wn4r, x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 2711 real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i 2712 real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i 2713 wn4r = w(1) 2714 x0r = a(0) + a(8) 2715 x0i = a(1) + a(9) 2716 x1r = a(0) - a(8) 2717 x1i = a(1) - a(9) 2718 x2r = a(4) + a(12) 2719 x2i = a(5) + a(13) 2720 x3r = a(4) - a(12) 2721 x3i = a(5) - a(13) 2722 y0r = x0r + x2r 2723 y0i = x0i + x2i 2724 y2r = x0r - x2r 2725 y2i = x0i - x2i 2726 y1r = x1r - x3i 2727 y1i = x1i + x3r 2728 y3r = x1r + x3i 2729 y3i = x1i - x3r 2730 x0r = a(2) + a(10) 2731 x0i = a(3) + a(11) 2732 x1r = a(2) - a(10) 2733 x1i = a(3) - a(11) 2734 x2r = a(6) + a(14) 2735 x2i = a(7) + a(15) 2736 x3r = a(6) - a(14) 2737 x3i = a(7) - a(15) 2738 y4r = x0r + x2r 2739 y4i = x0i + x2i 2740 y6r = x0r - x2r 2741 y6i = x0i - x2i 2742 x0r = x1r - x3i 2743 x0i = x1i + x3r 2744 x2r = x1r + x3i 2745 x2i = x1i - x3r 2746 y5r = wn4r * (x0r - x0i) 2747 y5i = wn4r * (x0r + x0i) 2748 y7r = wn4r * (x2r - x2i) 2749 y7i = wn4r * (x2r + x2i) 2750 a(8) = y1r + y5r 2751 a(9) = y1i + y5i 2752 a(10) = y1r - y5r 2753 a(11) = y1i - y5i 2754 a(12) = y3r - y7i 2755 a(13) = y3i + y7r 2756 a(14) = y3r + y7i 2757 a(15) = y3i - y7r 2758 a(0) = y0r + y4r 2759 a(1) = y0i + y4i 2760 a(2) = y0r - y4r 2761 a(3) = y0i - y4i 2762 a(4) = y2r - y6i 2763 a(5) = y2i + y6r 2764 a(6) = y2r + y6i 2765 a(7) = y2i - y6r 2766 end 2767! 2768 subroutine cftf082(a, w) 2769 real*8 a(0 : 15), w(0 : *) 2770 real*8 wn4r, wk1r, wk1i, x0r, x0i, x1r, x1i 2771 real*8 y0r, y0i, y1r, y1i, y2r, y2i, y3r, y3i 2772 real*8 y4r, y4i, y5r, y5i, y6r, y6i, y7r, y7i 2773 wn4r = w(1) 2774 wk1r = w(2) 2775 wk1i = w(3) 2776 y0r = a(0) - a(9) 2777 y0i = a(1) + a(8) 2778 y1r = a(0) + a(9) 2779 y1i = a(1) - a(8) 2780 x0r = a(4) - a(13) 2781 x0i = a(5) + a(12) 2782 y2r = wn4r * (x0r - x0i) 2783 y2i = wn4r * (x0i + x0r) 2784 x0r = a(4) + a(13) 2785 x0i = a(5) - a(12) 2786 y3r = wn4r * (x0r - x0i) 2787 y3i = wn4r * (x0i + x0r) 2788 x0r = a(2) - a(11) 2789 x0i = a(3) + a(10) 2790 y4r = wk1r * x0r - wk1i * x0i 2791 y4i = wk1r * x0i + wk1i * x0r 2792 x0r = a(2) + a(11) 2793 x0i = a(3) - a(10) 2794 y5r = wk1i * x0r - wk1r * x0i 2795 y5i = wk1i * x0i + wk1r * x0r 2796 x0r = a(6) - a(15) 2797 x0i = a(7) + a(14) 2798 y6r = wk1i * x0r - wk1r * x0i 2799 y6i = wk1i * x0i + wk1r * x0r 2800 x0r = a(6) + a(15) 2801 x0i = a(7) - a(14) 2802 y7r = wk1r * x0r - wk1i * x0i 2803 y7i = wk1r * x0i + wk1i * x0r 2804 x0r = y0r + y2r 2805 x0i = y0i + y2i 2806 x1r = y4r + y6r 2807 x1i = y4i + y6i 2808 a(0) = x0r + x1r 2809 a(1) = x0i + x1i 2810 a(2) = x0r - x1r 2811 a(3) = x0i - x1i 2812 x0r = y0r - y2r 2813 x0i = y0i - y2i 2814 x1r = y4r - y6r 2815 x1i = y4i - y6i 2816 a(4) = x0r - x1i 2817 a(5) = x0i + x1r 2818 a(6) = x0r + x1i 2819 a(7) = x0i - x1r 2820 x0r = y1r - y3i 2821 x0i = y1i + y3r 2822 x1r = y5r - y7r 2823 x1i = y5i - y7i 2824 a(8) = x0r + x1r 2825 a(9) = x0i + x1i 2826 a(10) = x0r - x1r 2827 a(11) = x0i - x1i 2828 x0r = y1r + y3i 2829 x0i = y1i - y3r 2830 x1r = y5r + y7r 2831 x1i = y5i + y7i 2832 a(12) = x0r - x1i 2833 a(13) = x0i + x1r 2834 a(14) = x0r + x1i 2835 a(15) = x0i - x1r 2836 end 2837! 2838 subroutine cftf040(a) 2839 real*8 a(0 : 7), x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 2840 x0r = a(0) + a(4) 2841 x0i = a(1) + a(5) 2842 x1r = a(0) - a(4) 2843 x1i = a(1) - a(5) 2844 x2r = a(2) + a(6) 2845 x2i = a(3) + a(7) 2846 x3r = a(2) - a(6) 2847 x3i = a(3) - a(7) 2848 a(0) = x0r + x2r 2849 a(1) = x0i + x2i 2850 a(2) = x1r - x3i 2851 a(3) = x1i + x3r 2852 a(4) = x0r - x2r 2853 a(5) = x0i - x2i 2854 a(6) = x1r + x3i 2855 a(7) = x1i - x3r 2856 end 2857! 2858 subroutine cftb040(a) 2859 real*8 a(0 : 7), x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i 2860 x0r = a(0) + a(4) 2861 x0i = a(1) + a(5) 2862 x1r = a(0) - a(4) 2863 x1i = a(1) - a(5) 2864 x2r = a(2) + a(6) 2865 x2i = a(3) + a(7) 2866 x3r = a(2) - a(6) 2867 x3i = a(3) - a(7) 2868 a(0) = x0r + x2r 2869 a(1) = x0i + x2i 2870 a(2) = x1r + x3i 2871 a(3) = x1i - x3r 2872 a(4) = x0r - x2r 2873 a(5) = x0i - x2i 2874 a(6) = x1r - x3i 2875 a(7) = x1i + x3r 2876 end 2877! 2878 subroutine cftx020(a) 2879 real*8 a(0 : 3), x0r, x0i 2880 x0r = a(0) - a(2) 2881 x0i = a(1) - a(3) 2882 a(0) = a(0) + a(2) 2883 a(1) = a(1) + a(3) 2884 a(2) = x0r 2885 a(3) = x0i 2886 end 2887! 2888 subroutine rftfsub(n, a, nc, c) 2889 integer n, nc, j, k, kk, ks, m 2890 real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr, xi, yr, yi 2891 m = n / 2 2892 ks = 2 * nc / m 2893 kk = 0 2894 do j = 2, m - 2, 2 2895 k = n - j 2896 kk = kk + ks 2897 wkr = 0.5d0 - c(nc - kk) 2898 wki = c(kk) 2899 xr = a(j) - a(k) 2900 xi = a(j + 1) + a(k + 1) 2901 yr = wkr * xr - wki * xi 2902 yi = wkr * xi + wki * xr 2903 a(j) = a(j) - yr 2904 a(j + 1) = a(j + 1) - yi 2905 a(k) = a(k) + yr 2906 a(k + 1) = a(k + 1) - yi 2907 end do 2908 end 2909! 2910 subroutine rftbsub(n, a, nc, c) 2911 integer n, nc, j, k, kk, ks, m 2912 real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr, xi, yr, yi 2913 m = n / 2 2914 ks = 2 * nc / m 2915 kk = 0 2916 do j = 2, m - 2, 2 2917 k = n - j 2918 kk = kk + ks 2919 wkr = 0.5d0 - c(nc - kk) 2920 wki = c(kk) 2921 xr = a(j) - a(k) 2922 xi = a(j + 1) + a(k + 1) 2923 yr = wkr * xr + wki * xi 2924 yi = wkr * xi - wki * xr 2925 a(j) = a(j) - yr 2926 a(j + 1) = a(j + 1) - yi 2927 a(k) = a(k) + yr 2928 a(k + 1) = a(k + 1) - yi 2929 end do 2930 end 2931! 2932 subroutine dctsub(n, a, nc, c) 2933 integer n, nc, j, k, kk, ks, m 2934 real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr 2935 m = n / 2 2936 ks = nc / n 2937 kk = 0 2938 do j = 1, m - 1 2939 k = n - j 2940 kk = kk + ks 2941 wkr = c(kk) - c(nc - kk) 2942 wki = c(kk) + c(nc - kk) 2943 xr = wki * a(j) - wkr * a(k) 2944 a(j) = wkr * a(j) + wki * a(k) 2945 a(k) = xr 2946 end do 2947 a(m) = c(0) * a(m) 2948 end 2949! 2950 subroutine dstsub(n, a, nc, c) 2951 integer n, nc, j, k, kk, ks, m 2952 real*8 a(0 : n - 1), c(0 : nc - 1), wkr, wki, xr 2953 m = n / 2 2954 ks = nc / n 2955 kk = 0 2956 do j = 1, m - 1 2957 k = n - j 2958 kk = kk + ks 2959 wkr = c(kk) - c(nc - kk) 2960 wki = c(kk) + c(nc - kk) 2961 xr = wki * a(k) - wkr * a(j) 2962 a(k) = wkr * a(k) + wki * a(j) 2963 a(j) = xr 2964 end do 2965 a(m) = c(0) * a(m) 2966 end 2967! 2968