1! Fast Fourier/Cosine/Sine Transform 2! dimension :two 3! data length :power of 2 4! decimation :frequency 5! radix :split-radix, row-column 6! data :inplace 7! table :use 8! subroutines 9! cdft2d: Complex Discrete Fourier Transform 10! rdft2d: Real Discrete Fourier Transform 11! ddct2d: Discrete Cosine Transform 12! ddst2d: Discrete Sine Transform 13! necessary package 14! fftsg.f : 1D-FFT package 15! 16! 17! -------- Complex DFT (Discrete Fourier Transform) -------- 18! [definition] 19! <case1> 20! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) * 21! exp(2*pi*i*j1*k1/n1) * 22! exp(2*pi*i*j2*k2/n2), 23! 0<=k1<n1, 0<=k2<n2 24! <case2> 25! X(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 x(j1,j2) * 26! exp(-2*pi*i*j1*k1/n1) * 27! exp(-2*pi*i*j2*k2/n2), 28! 0<=k1<n1, 0<=k2<n2 29! (notes: sum_j=0^n-1 is a summation from j=0 to n-1) 30! [usage] 31! <case1> 32! ip(0) = 0 ! first time only 33! call cdft2d(n1max, 2*n1, n2, 1, a, t, ip, w) 34! <case2> 35! ip(0) = 0 ! first time only 36! call cdft2d(n1max, 2*n1, n2, -1, a, t, ip, w) 37! [parameters] 38! n1max :row size of the 2D array (integer) 39! 2*n1 :data length (integer) 40! n1 >= 1, n1 = power of 2 41! n2 :data length (integer) 42! n2 >= 1, n2 = power of 2 43! a(0:2*n1-1,0:n2-1) 44! :input/output data (real*8) 45! input data 46! a(2*j1,j2) = Re(x(j1,j2)), 47! a(2*j1+1,j2) = Im(x(j1,j2)), 48! 0<=j1<n1, 0<=j2<n2 49! output data 50! a(2*k1,k2) = Re(X(k1,k2)), 51! a(2*k1+1,k2) = Im(X(k1,k2)), 52! 0<=k1<n1, 0<=k2<n2 53! t(0:8*n2-1) 54! :work area (real*8) 55! length of t >= 8*n2 56! ip(0:*):work area for bit reversal (integer) 57! length of ip >= 2+sqrt(n) 58! (n = max(n1, n2)) 59! ip(0),ip(1) are pointers of the cos/sin table. 60! w(0:*) :cos/sin table (real*8) 61! length of w >= max(n1/2, n2/2) 62! w(),ip() are initialized if ip(0) = 0. 63! [remark] 64! Inverse of 65! call cdft2d(n1max, 2*n1, n2, -1, a, t, ip, w) 66! is 67! call cdft2d(n1max, 2*n1, n2, 1, a, t, ip, w) 68! do j2 = 0, n2 - 1 69! do j1 = 0, 2 * n1 - 1 70! a(j1, j2) = a(j1, j2) * (1.0d0 / n1 / n2) 71! end do 72! end do 73! . 74! 75! 76! -------- Real DFT / Inverse of Real DFT -------- 77! [definition] 78! <case1> RDFT 79! R(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 80! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), 81! 0<=k1<n1, 0<=k2<n2 82! I(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 83! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2), 84! 0<=k1<n1, 0<=k2<n2 85! <case2> IRDFT (excluding scale) 86! a(k1,k2) = (1/2) * sum_j1=0^n1-1 sum_j2=0^n2-1 87! (R(j1,j2) * 88! cos(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2) + 89! I(j1,j2) * 90! sin(2*pi*j1*k1/n1 + 2*pi*j2*k2/n2)), 91! 0<=k1<n1, 0<=k2<n2 92! (notes: R(n1-k1,n2-k2) = R(k1,k2), 93! I(n1-k1,n2-k2) = -I(k1,k2), 94! R(n1-k1,0) = R(k1,0), 95! I(n1-k1,0) = -I(k1,0), 96! R(0,n2-k2) = R(0,k2), 97! I(0,n2-k2) = -I(0,k2), 98! 0<k1<n1, 0<k2<n2) 99! [usage] 100! <case1> 101! ip(0) = 0 ! first time only 102! call rdft2d(n1max, n1, n2, 1, a, t, ip, w) 103! <case2> 104! ip(0) = 0 ! first time only 105! call rdft2d(n1max, n1, n2, -1, a, t, ip, w) 106! [parameters] 107! n1max :row size of the 2D array (integer) 108! n1 :data length (integer) 109! n1 >= 2, n1 = power of 2 110! n2 :data length (integer) 111! n2 >= 2, n2 = power of 2 112! a(0:n1-1,0:n2-1) 113! :input/output data (real*8) 114! <case1> 115! output data 116! a(2*k1,k2) = R(k1,k2) = R(n1-k1,n2-k2), 117! a(2*k1+1,k2) = I(k1,k2) = -I(n1-k1,n2-k2), 118! 0<k1<n1/2, 0<k2<n2, 119! a(2*k1,0) = R(k1,0) = R(n1-k1,0), 120! a(2*k1+1,0) = I(k1,0) = -I(n1-k1,0), 121! 0<k1<n1/2, 122! a(0,k2) = R(0,k2) = R(0,n2-k2), 123! a(1,k2) = I(0,k2) = -I(0,n2-k2), 124! a(1,n2-k2) = R(n1/2,k2) = R(n1/2,n2-k2), 125! a(0,n2-k2) = -I(n1/2,k2) = I(n1/2,n2-k2), 126! 0<k2<n2/2, 127! a(0,0) = R(0,0), 128! a(1,0) = R(n1/2,0), 129! a(0,n2/2) = R(0,n2/2), 130! a(1,n2/2) = R(n1/2,n2/2) 131! <case2> 132! input data 133! a(2*j1,j2) = R(j1,j2) = R(n1-j1,n2-j2), 134! a(2*j1+1,j2) = I(j1,j2) = -I(n1-j1,n2-j2), 135! 0<j1<n1/2, 0<j2<n2, 136! a(2*j1,0) = R(j1,0) = R(n1-j1,0), 137! a(2*j1+1,0) = I(j1,0) = -I(n1-j1,0), 138! 0<j1<n1/2, 139! a(0,j2) = R(0,j2) = R(0,n2-j2), 140! a(1,j2) = I(0,j2) = -I(0,n2-j2), 141! a(1,n2-j2) = R(n1/2,j2) = R(n1/2,n2-j2), 142! a(0,n2-j2) = -I(n1/2,j2) = I(n1/2,n2-j2), 143! 0<j2<n2/2, 144! a(0,0) = R(0,0), 145! a(1,0) = R(n1/2,0), 146! a(0,n2/2) = R(0,n2/2), 147! a(1,n2/2) = R(n1/2,n2/2) 148! ---- output ordering ---- 149! call rdft2d(n1max, n1, n2, 1, a, t, ip, w) 150! call rdft2dsort(n1max, n1, n2, 1, a) 151! ! stored data is a(0:n1-1,0:n2+1): 152! ! a(2*k1,k2) = R(k1,k2), 153! ! a(2*k1+1,k2) = I(k1,k2), 154! ! 0<=k1<=n1/2, 0<=k2<n2. 155! ! the stored data is larger than the input data! 156! ---- input ordering ---- 157! call rdft2dsort(n1max, n1, n2, -1, a) 158! call rdft2d(n1max, n1, n2, -1, a, t, ip, w) 159! t(0:8*n2-1) 160! :work area (real*8) 161! length of t >= 8*n2 162! ip(0:*):work area for bit reversal (integer) 163! length of ip >= 2+sqrt(n) 164! (n = max(n1/2, n2)) 165! ip(0),ip(1) are pointers of the cos/sin table. 166! w(0:*) :cos/sin table (real*8) 167! length of w >= max(n1/4, n2/2) + n1/4 168! w(),ip() are initialized if ip(0) = 0. 169! [remark] 170! Inverse of 171! call rdft2d(n1max, n1, n2, 1, a, t, ip, w) 172! is 173! call rdft2d(n1max, n1, n2, -1, a, t, ip, w) 174! do j2 = 0, n2 - 1 175! do j1 = 0, n1 - 1 176! a(j1, j2) = a(j1, j2) * (2.0d0 / n1 / n2) 177! end do 178! end do 179! . 180! 181! 182! -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- 183! [definition] 184! <case1> IDCT (excluding scale) 185! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 186! cos(pi*j1*(k1+1/2)/n1) * 187! cos(pi*j2*(k2+1/2)/n2), 188! 0<=k1<n1, 0<=k2<n2 189! <case2> DCT 190! C(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 191! cos(pi*(j1+1/2)*k1/n1) * 192! cos(pi*(j2+1/2)*k2/n2), 193! 0<=k1<n1, 0<=k2<n2 194! [usage] 195! <case1> 196! ip(0) = 0 ! first time only 197! call ddct2d(n1max, n1, n2, 1, a, t, ip, w) 198! <case2> 199! ip(0) = 0 ! first time only 200! call ddct2d(n1max, n1, n2, -1, a, t, ip, w) 201! [parameters] 202! n1max :row size of the 2D array (integer) 203! n1 :data length (integer) 204! n1 >= 2, n1 = power of 2 205! n2 :data length (integer) 206! n2 >= 2, n2 = power of 2 207! a(0:n1-1,0:n2-1) 208! :input/output data (real*8) 209! output data 210! a(k1,k2) = C(k1,k2), 0<=k1<n1, 0<=k2<n2 211! t(0:4*n2-1) 212! :work area (real*8) 213! length of t >= 4*n2 214! ip(0:*):work area for bit reversal (integer) 215! length of ip >= 2+sqrt(n) 216! (n = max(n1/2, n2/2)) 217! ip(0),ip(1) are pointers of the cos/sin table. 218! w(0:*) :cos/sin table (real*8) 219! length of w >= max(n1*3/2, n2*3/2) 220! w(),ip() are initialized if ip(0) = 0. 221! [remark] 222! Inverse of 223! call ddct2d(n1max, n1, n2, -1, a, t, ip, w) 224! is 225! do j1 = 0, n1 - 1 226! a(j1, 0) = a(j1, 0) * 0.5d0 227! end do 228! do j2 = 0, n2 - 1 229! a(0, j2) = a(0, j2) * 0.5d0 230! end do 231! call ddct2d(n1max, n1, n2, 1, a, t, ip, w) 232! do j2 = 0, n2 - 1 233! do j1 = 0, n1 - 1 234! a(j1, j2) = a(j1, j2) * (4.0d0 / n1 / n2) 235! end do 236! end do 237! . 238! 239! 240! -------- DST (Discrete Sine Transform) / Inverse of DST -------- 241! [definition] 242! <case1> IDST (excluding scale) 243! S(k1,k2) = sum_j1=1^n1 sum_j2=1^n2 A(j1,j2) * 244! sin(pi*j1*(k1+1/2)/n1) * 245! sin(pi*j2*(k2+1/2)/n2), 246! 0<=k1<n1, 0<=k2<n2 247! <case2> DST 248! S(k1,k2) = sum_j1=0^n1-1 sum_j2=0^n2-1 a(j1,j2) * 249! sin(pi*(j1+1/2)*k1/n1) * 250! sin(pi*(j2+1/2)*k2/n2), 251! 0<k1<=n1, 0<k2<=n2 252! [usage] 253! <case1> 254! ip(0) = 0 ! first time only 255! call ddst2d(n1max, n1, n2, 1, a, t, ip, w) 256! <case2> 257! ip(0) = 0 ! first time only 258! call ddst2d(n1max, n1, n2, -1, a, t, ip, w) 259! [parameters] 260! n1max :row size of the 2D array (integer) 261! n1 :data length (integer) 262! n1 >= 2, n1 = power of 2 263! n2 :data length (integer) 264! n2 >= 2, n2 = power of 2 265! a(0:n1-1,0:n2-1) 266! :input/output data (real*8) 267! <case1> 268! input data 269! a(j1,j2) = A(j1,j2), 0<j1<n1, 0<j2<n2, 270! a(j1,0) = A(j1,n2), 0<j1<n1, 271! a(0,j2) = A(n1,j2), 0<j2<n2, 272! a(0,0) = A(n1,n2) 273! (i.e. A(j1,j2) = a(mod(j1,n1),mod(j2,n2))) 274! output data 275! a(k1,k2) = S(k1,k2), 0<=k1<n1, 0<=k2<n2 276! <case2> 277! output data 278! a(k1,k2) = S(k1,k2), 0<k1<n1, 0<k2<n2, 279! a(k1,0) = S(k1,n2), 0<k1<n1, 280! a(0,k2) = S(n1,k2), 0<k2<n2, 281! a(0,0) = S(n1,n2) 282! (i.e. S(k1,k2) = a(mod(k1,n1),mod(k2,n2))) 283! t(0:4*n2-1) 284! :work area (real*8) 285! length of t >= 4*n2 286! ip(0:*):work area for bit reversal (integer) 287! length of ip >= 2+sqrt(n) 288! (n = max(n1/2, n2/2)) 289! ip(0),ip(1) are pointers of the cos/sin table. 290! w(0:*) :cos/sin table (real*8) 291! length of w >= max(n1*3/2, n2*3/2) 292! w(),ip() are initialized if ip(0) = 0. 293! [remark] 294! Inverse of 295! call ddst2d(n1max, n1, n2, -1, a, t, ip, w) 296! is 297! do j1 = 0, n1 - 1 298! a(j1, 0) = a(j1, 0) * 0.5d0 299! end do 300! do j2 = 0, n2 - 1 301! a(0, j2) = a(0, j2) * 0.5d0 302! end do 303! call ddst2d(n1max, n1, n2, 1, a, t, ip, w) 304! do j2 = 0, n2 - 1 305! do j1 = 0, n1 - 1 306! a(j1, j2) = a(j1, j2) * (4.0d0 / n1 / n2) 307! end do 308! end do 309! . 310! 311! 312 subroutine cdft2d(n1max, n1, n2, isgn, a, t, ip, w) 313 integer n1max, n1, n2, isgn, ip(0 : *), n, j 314 real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 * n2 - 1), 315 & w(0 : *) 316 n = max(n1, 2 * n2) 317 if (n .gt. 4 * ip(0)) then 318 call makewt(n / 4, ip, w) 319 end if 320 do j = 0, n2 - 1 321 call cdft(n1, isgn, a(0, j), ip, w) 322 end do 323 call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w) 324 end 325! 326 subroutine rdft2d(n1max, n1, n2, isgn, a, t, ip, w) 327 integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j 328 real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 * n2 - 1), 329 & w(0 : *) 330 n = max(n1, 2 * n2) 331 nw = ip(0) 332 if (n .gt. 4 * nw) then 333 nw = n / 4 334 call makewt(nw, ip, w) 335 end if 336 nc = ip(1) 337 if (n1 .gt. 4 * nc) then 338 nc = n1 / 4 339 call makect(nc, ip, w(nw)) 340 end if 341 if (isgn .lt. 0) then 342 call rdft2d_sub(n1max, n1, n2, isgn, a) 343 call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w) 344 end if 345 do j = 0, n2 - 1 346 call rdft(n1, isgn, a(0, j), ip, w) 347 end do 348 if (isgn .ge. 0) then 349 call cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w) 350 call rdft2d_sub(n1max, n1, n2, isgn, a) 351 end if 352 end 353! 354 subroutine rdft2dsort(n1max, n1, n2, isgn, a) 355 integer n1max, n1, n2, isgn, n2h, j 356 real*8 a(0 : n1max - 1, 0 : n2 - 1), x, y 357 n2h = n2 / 2 358 if (isgn .lt. 0) then 359 do j = n2h + 1, n2 - 1 360 a(0, j) = a(n1 + 1, j) 361 a(1, j) = a(n1, j) 362 end do 363 a(1, 0) = a(n1, 0) 364 a(1, n2h) = a(n1, n2h) 365 else 366 do j = n2h + 1, n2 - 1 367 y = a(0, j) 368 x = a(1, j) 369 a(n1, j) = x 370 a(n1 + 1, j) = y 371 a(n1, n2 - j) = x 372 a(n1 + 1, n2 - j) = -y 373 a(0, j) = a(0, n2 - j) 374 a(1, j) = -a(1, n2 - j) 375 end do 376 a(n1, 0) = a(1, 0) 377 a(n1 + 1, 0) = 0 378 a(1, 0) = 0 379 a(n1, n2h) = a(1, n2h) 380 a(n1 + 1, n2h) = 0 381 a(1, n2h) = 0 382 end if 383 end 384! 385 subroutine ddct2d(n1max, n1, n2, isgn, a, t, ip, w) 386 integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j 387 real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 * n2 - 1), 388 & w(0 : *) 389 n = max(n1, n2) 390 nw = ip(0) 391 if (n .gt. 4 * nw) then 392 nw = n / 4 393 call makewt(nw, ip, w) 394 end if 395 nc = ip(1) 396 if (n .gt. nc) then 397 nc = n 398 call makect(nc, ip, w(nw)) 399 end if 400 do j = 0, n2 - 1 401 call ddct(n1, isgn, a(0, j), ip, w) 402 end do 403 call ddxt2d_sub(n1max, n1, n2, 0, isgn, a, t, ip, w) 404 end 405! 406 subroutine ddst2d(n1max, n1, n2, isgn, a, t, ip, w) 407 integer n1max, n1, n2, isgn, ip(0 : *), n, nw, nc, j 408 real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 * n2 - 1), 409 & w(0 : *) 410 n = max(n1, n2) 411 nw = ip(0) 412 if (n .gt. 4 * nw) then 413 nw = n / 4 414 call makewt(nw, ip, w) 415 end if 416 nc = ip(1) 417 if (n .gt. nc) then 418 nc = n 419 call makect(nc, ip, w(nw)) 420 end if 421 do j = 0, n2 - 1 422 call ddst(n1, isgn, a(0, j), ip, w) 423 end do 424 call ddxt2d_sub(n1max, n1, n2, 1, isgn, a, t, ip, w) 425 end 426! 427! -------- child routines -------- 428! 429 subroutine cdft2d_sub(n1max, n1, n2, isgn, a, t, ip, w) 430 integer n1max, n1, n2, isgn, ip(0 : *), i, j 431 real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 8 * n2 - 1), 432 & w(0 : *) 433 if (n1 .gt. 4) then 434 do i = 0, n1 - 8, 8 435 do j = 0, n2 - 1 436 t(2 * j) = a(i, j) 437 t(2 * j + 1) = a(i + 1, j) 438 t(2 * n2 + 2 * j) = a(i + 2, j) 439 t(2 * n2 + 2 * j + 1) = a(i + 3, j) 440 t(4 * n2 + 2 * j) = a(i + 4, j) 441 t(4 * n2 + 2 * j + 1) = a(i + 5, j) 442 t(6 * n2 + 2 * j) = a(i + 6, j) 443 t(6 * n2 + 2 * j + 1) = a(i + 7, j) 444 end do 445 call cdft(2 * n2, isgn, t, ip, w) 446 call cdft(2 * n2, isgn, t(2 * n2), ip, w) 447 call cdft(2 * n2, isgn, t(4 * n2), ip, w) 448 call cdft(2 * n2, isgn, t(6 * n2), ip, w) 449 do j = 0, n2 - 1 450 a(i, j) = t(2 * j) 451 a(i + 1, j) = t(2 * j + 1) 452 a(i + 2, j) = t(2 * n2 + 2 * j) 453 a(i + 3, j) = t(2 * n2 + 2 * j + 1) 454 a(i + 4, j) = t(4 * n2 + 2 * j) 455 a(i + 5, j) = t(4 * n2 + 2 * j + 1) 456 a(i + 6, j) = t(6 * n2 + 2 * j) 457 a(i + 7, j) = t(6 * n2 + 2 * j + 1) 458 end do 459 end do 460 else if (n1 .eq. 4) then 461 do j = 0, n2 - 1 462 t(2 * j) = a(0, j) 463 t(2 * j + 1) = a(1, j) 464 t(2 * n2 + 2 * j) = a(2, j) 465 t(2 * n2 + 2 * j + 1) = a(3, j) 466 end do 467 call cdft(2 * n2, isgn, t, ip, w) 468 call cdft(2 * n2, isgn, t(2 * n2), ip, w) 469 do j = 0, n2 - 1 470 a(0, j) = t(2 * j) 471 a(1, j) = t(2 * j + 1) 472 a(2, j) = t(2 * n2 + 2 * j) 473 a(3, j) = t(2 * n2 + 2 * j + 1) 474 end do 475 else if (n1 .eq. 2) then 476 do j = 0, n2 - 1 477 t(2 * j) = a(0, j) 478 t(2 * j + 1) = a(1, j) 479 end do 480 call cdft(2 * n2, isgn, t, ip, w) 481 do j = 0, n2 - 1 482 a(0, j) = t(2 * j) 483 a(1, j) = t(2 * j + 1) 484 end do 485 end if 486 end 487! 488 subroutine rdft2d_sub(n1max, n1, n2, isgn, a) 489 integer n1max, n1, n2, isgn, n2h, i, j 490 real*8 a(0 : n1max - 1, 0 : n2 - 1), xi 491 n2h = n2 / 2 492 if (isgn .lt. 0) then 493 do i = 1, n2h - 1 494 j = n2 - i 495 xi = a(0, i) - a(0, j) 496 a(0, i) = a(0, i) + a(0, j) 497 a(0, j) = xi 498 xi = a(1, j) - a(1, i) 499 a(1, i) = a(1, i) + a(1, j) 500 a(1, j) = xi 501 end do 502 else 503 do i = 1, n2h - 1 504 j = n2 - i 505 a(0, j) = 0.5d0 * (a(0, i) - a(0, j)) 506 a(0, i) = a(0, i) - a(0, j) 507 a(1, j) = 0.5d0 * (a(1, i) + a(1, j)) 508 a(1, i) = a(1, i) - a(1, j) 509 end do 510 end if 511 end 512! 513 subroutine ddxt2d_sub(n1max, n1, n2, ics, isgn, a, t, 514 & ip, w) 515 integer n1max, n1, n2, ics, isgn, ip(0 : *), i, j 516 real*8 a(0 : n1max - 1, 0 : n2 - 1), t(0 : 4 * n2 - 1), 517 & w(0 : *) 518 if (n1 .gt. 2) then 519 do i = 0, n1 - 4, 4 520 do j = 0, n2 - 1 521 t(j) = a(i, j) 522 t(n2 + j) = a(i + 1, j) 523 t(2 * n2 + j) = a(i + 2, j) 524 t(3 * n2 + j) = a(i + 3, j) 525 end do 526 if (ics .eq. 0) then 527 call ddct(n2, isgn, t, ip, w) 528 call ddct(n2, isgn, t(n2), ip, w) 529 call ddct(n2, isgn, t(2 * n2), ip, w) 530 call ddct(n2, isgn, t(3 * n2), ip, w) 531 else 532 call ddst(n2, isgn, t, ip, w) 533 call ddst(n2, isgn, t(n2), ip, w) 534 call ddst(n2, isgn, t(2 * n2), ip, w) 535 call ddst(n2, isgn, t(3 * n2), ip, w) 536 end if 537 do j = 0, n2 - 1 538 a(i, j) = t(j) 539 a(i + 1, j) = t(n2 + j) 540 a(i + 2, j) = t(2 * n2 + j) 541 a(i + 3, j) = t(3 * n2 + j) 542 end do 543 end do 544 else if (n1 .eq. 2) then 545 do j = 0, n2 - 1 546 t(j) = a(0, j) 547 t(n2 + j) = a(1, j) 548 end do 549 if (ics .eq. 0) then 550 call ddct(n2, isgn, t, ip, w) 551 call ddct(n2, isgn, t(n2), ip, w) 552 else 553 call ddst(n2, isgn, t, ip, w) 554 call ddst(n2, isgn, t(n2), ip, w) 555 end if 556 do j = 0, n2 - 1 557 a(0, j) = t(j) 558 a(1, j) = t(n2 + j) 559 end do 560 end if 561 end 562! 563