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