1 /* 2 * http://www.kurims.kyoto-u.ac.jp/~ooura/fft.html 3 * Copyright Takuya OOURA, 1996-2001 4 * 5 * You may use, copy, modify and distribute this code for any purpose (include 6 * commercial use) and without fee. Please refer to this package when you modify 7 * this code. 8 * 9 * Changes: 10 * Trivial type modifications by the WebRTC authors. 11 */ 12 13 /* 14 Fast Fourier/Cosine/Sine Transform 15 dimension :one 16 data length :power of 2 17 decimation :frequency 18 radix :4, 2 19 data :inplace 20 table :use 21 functions 22 cdft: Complex Discrete Fourier Transform 23 rdft: Real Discrete Fourier Transform 24 ddct: Discrete Cosine Transform 25 ddst: Discrete Sine Transform 26 dfct: Cosine Transform of RDFT (Real Symmetric DFT) 27 dfst: Sine Transform of RDFT (Real Anti-symmetric DFT) 28 function prototypes 29 void cdft(int, int, float *, int *, float *); 30 void rdft(int, int, float *, int *, float *); 31 void ddct(int, int, float *, int *, float *); 32 void ddst(int, int, float *, int *, float *); 33 void dfct(int, float *, float *, int *, float *); 34 void dfst(int, float *, float *, int *, float *); 35 36 37 -------- Complex DFT (Discrete Fourier Transform) -------- 38 [definition] 39 <case1> 40 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n 41 <case2> 42 X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n 43 (notes: sum_j=0^n-1 is a summation from j=0 to n-1) 44 [usage] 45 <case1> 46 ip[0] = 0; // first time only 47 cdft(2*n, 1, a, ip, w); 48 <case2> 49 ip[0] = 0; // first time only 50 cdft(2*n, -1, a, ip, w); 51 [parameters] 52 2*n :data length (int) 53 n >= 1, n = power of 2 54 a[0...2*n-1] :input/output data (float *) 55 input data 56 a[2*j] = Re(x[j]), 57 a[2*j+1] = Im(x[j]), 0<=j<n 58 output data 59 a[2*k] = Re(X[k]), 60 a[2*k+1] = Im(X[k]), 0<=k<n 61 ip[0...*] :work area for bit reversal (int *) 62 length of ip >= 2+sqrt(n) 63 strictly, 64 length of ip >= 65 2+(1<<(int)(log(n+0.5)/log(2))/2). 66 ip[0],ip[1] are pointers of the cos/sin table. 67 w[0...n/2-1] :cos/sin table (float *) 68 w[],ip[] are initialized if ip[0] == 0. 69 [remark] 70 Inverse of 71 cdft(2*n, -1, a, ip, w); 72 is 73 cdft(2*n, 1, a, ip, w); 74 for (j = 0; j <= 2 * n - 1; j++) { 75 a[j] *= 1.0 / n; 76 } 77 . 78 79 80 -------- Real DFT / Inverse of Real DFT -------- 81 [definition] 82 <case1> RDFT 83 R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2 84 I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2 85 <case2> IRDFT (excluding scale) 86 a[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 87 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 88 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n 89 [usage] 90 <case1> 91 ip[0] = 0; // first time only 92 rdft(n, 1, a, ip, w); 93 <case2> 94 ip[0] = 0; // first time only 95 rdft(n, -1, a, ip, w); 96 [parameters] 97 n :data length (int) 98 n >= 2, n = power of 2 99 a[0...n-1] :input/output data (float *) 100 <case1> 101 output data 102 a[2*k] = R[k], 0<=k<n/2 103 a[2*k+1] = I[k], 0<k<n/2 104 a[1] = R[n/2] 105 <case2> 106 input data 107 a[2*j] = R[j], 0<=j<n/2 108 a[2*j+1] = I[j], 0<j<n/2 109 a[1] = R[n/2] 110 ip[0...*] :work area for bit reversal (int *) 111 length of ip >= 2+sqrt(n/2) 112 strictly, 113 length of ip >= 114 2+(1<<(int)(log(n/2+0.5)/log(2))/2). 115 ip[0],ip[1] are pointers of the cos/sin table. 116 w[0...n/2-1] :cos/sin table (float *) 117 w[],ip[] are initialized if ip[0] == 0. 118 [remark] 119 Inverse of 120 rdft(n, 1, a, ip, w); 121 is 122 rdft(n, -1, a, ip, w); 123 for (j = 0; j <= n - 1; j++) { 124 a[j] *= 2.0 / n; 125 } 126 . 127 128 129 -------- DCT (Discrete Cosine Transform) / Inverse of DCT -------- 130 [definition] 131 <case1> IDCT (excluding scale) 132 C[k] = sum_j=0^n-1 a[j]*cos(pi*j*(k+1/2)/n), 0<=k<n 133 <case2> DCT 134 C[k] = sum_j=0^n-1 a[j]*cos(pi*(j+1/2)*k/n), 0<=k<n 135 [usage] 136 <case1> 137 ip[0] = 0; // first time only 138 ddct(n, 1, a, ip, w); 139 <case2> 140 ip[0] = 0; // first time only 141 ddct(n, -1, a, ip, w); 142 [parameters] 143 n :data length (int) 144 n >= 2, n = power of 2 145 a[0...n-1] :input/output data (float *) 146 output data 147 a[k] = C[k], 0<=k<n 148 ip[0...*] :work area for bit reversal (int *) 149 length of ip >= 2+sqrt(n/2) 150 strictly, 151 length of ip >= 152 2+(1<<(int)(log(n/2+0.5)/log(2))/2). 153 ip[0],ip[1] are pointers of the cos/sin table. 154 w[0...n*5/4-1] :cos/sin table (float *) 155 w[],ip[] are initialized if ip[0] == 0. 156 [remark] 157 Inverse of 158 ddct(n, -1, a, ip, w); 159 is 160 a[0] *= 0.5; 161 ddct(n, 1, a, ip, w); 162 for (j = 0; j <= n - 1; j++) { 163 a[j] *= 2.0 / n; 164 } 165 . 166 167 168 -------- DST (Discrete Sine Transform) / Inverse of DST -------- 169 [definition] 170 <case1> IDST (excluding scale) 171 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n 172 <case2> DST 173 S[k] = sum_j=0^n-1 a[j]*sin(pi*(j+1/2)*k/n), 0<k<=n 174 [usage] 175 <case1> 176 ip[0] = 0; // first time only 177 ddst(n, 1, a, ip, w); 178 <case2> 179 ip[0] = 0; // first time only 180 ddst(n, -1, a, ip, w); 181 [parameters] 182 n :data length (int) 183 n >= 2, n = power of 2 184 a[0...n-1] :input/output data (float *) 185 <case1> 186 input data 187 a[j] = A[j], 0<j<n 188 a[0] = A[n] 189 output data 190 a[k] = S[k], 0<=k<n 191 <case2> 192 output data 193 a[k] = S[k], 0<k<n 194 a[0] = S[n] 195 ip[0...*] :work area for bit reversal (int *) 196 length of ip >= 2+sqrt(n/2) 197 strictly, 198 length of ip >= 199 2+(1<<(int)(log(n/2+0.5)/log(2))/2). 200 ip[0],ip[1] are pointers of the cos/sin table. 201 w[0...n*5/4-1] :cos/sin table (float *) 202 w[],ip[] are initialized if ip[0] == 0. 203 [remark] 204 Inverse of 205 ddst(n, -1, a, ip, w); 206 is 207 a[0] *= 0.5; 208 ddst(n, 1, a, ip, w); 209 for (j = 0; j <= n - 1; j++) { 210 a[j] *= 2.0 / n; 211 } 212 . 213 214 215 -------- Cosine Transform of RDFT (Real Symmetric DFT) -------- 216 [definition] 217 C[k] = sum_j=0^n a[j]*cos(pi*j*k/n), 0<=k<=n 218 [usage] 219 ip[0] = 0; // first time only 220 dfct(n, a, t, ip, w); 221 [parameters] 222 n :data length - 1 (int) 223 n >= 2, n = power of 2 224 a[0...n] :input/output data (float *) 225 output data 226 a[k] = C[k], 0<=k<=n 227 t[0...n/2] :work area (float *) 228 ip[0...*] :work area for bit reversal (int *) 229 length of ip >= 2+sqrt(n/4) 230 strictly, 231 length of ip >= 232 2+(1<<(int)(log(n/4+0.5)/log(2))/2). 233 ip[0],ip[1] are pointers of the cos/sin table. 234 w[0...n*5/8-1] :cos/sin table (float *) 235 w[],ip[] are initialized if ip[0] == 0. 236 [remark] 237 Inverse of 238 a[0] *= 0.5; 239 a[n] *= 0.5; 240 dfct(n, a, t, ip, w); 241 is 242 a[0] *= 0.5; 243 a[n] *= 0.5; 244 dfct(n, a, t, ip, w); 245 for (j = 0; j <= n; j++) { 246 a[j] *= 2.0 / n; 247 } 248 . 249 250 251 -------- Sine Transform of RDFT (Real Anti-symmetric DFT) -------- 252 [definition] 253 S[k] = sum_j=1^n-1 a[j]*sin(pi*j*k/n), 0<k<n 254 [usage] 255 ip[0] = 0; // first time only 256 dfst(n, a, t, ip, w); 257 [parameters] 258 n :data length + 1 (int) 259 n >= 2, n = power of 2 260 a[0...n-1] :input/output data (float *) 261 output data 262 a[k] = S[k], 0<k<n 263 (a[0] is used for work area) 264 t[0...n/2-1] :work area (float *) 265 ip[0...*] :work area for bit reversal (int *) 266 length of ip >= 2+sqrt(n/4) 267 strictly, 268 length of ip >= 269 2+(1<<(int)(log(n/4+0.5)/log(2))/2). 270 ip[0],ip[1] are pointers of the cos/sin table. 271 w[0...n*5/8-1] :cos/sin table (float *) 272 w[],ip[] are initialized if ip[0] == 0. 273 [remark] 274 Inverse of 275 dfst(n, a, t, ip, w); 276 is 277 dfst(n, a, t, ip, w); 278 for (j = 1; j <= n - 1; j++) { 279 a[j] *= 2.0 / n; 280 } 281 . 282 283 284 Appendix : 285 The cos/sin table is recalculated when the larger table required. 286 w[] and ip[] are compatible with all routines. 287 */ 288 289 void cdft(int n, int isgn, float *a, int *ip, float *w) 290 { 291 void makewt(int nw, int *ip, float *w); 292 void bitrv2(int n, int *ip, float *a); 293 void bitrv2conj(int n, int *ip, float *a); 294 void cftfsub(int n, float *a, float *w); 295 void cftbsub(int n, float *a, float *w); 296 297 if (n > (ip[0] << 2)) { 298 makewt(n >> 2, ip, w); 299 } 300 if (n > 4) { 301 if (isgn >= 0) { 302 bitrv2(n, ip + 2, a); 303 cftfsub(n, a, w); 304 } else { 305 bitrv2conj(n, ip + 2, a); 306 cftbsub(n, a, w); 307 } 308 } else if (n == 4) { 309 cftfsub(n, a, w); 310 } 311 } 312 313 314 void rdft(int n, int isgn, float *a, int *ip, float *w) 315 { 316 void makewt(int nw, int *ip, float *w); 317 void makect(int nc, int *ip, float *c); 318 void bitrv2(int n, int *ip, float *a); 319 void cftfsub(int n, float *a, float *w); 320 void cftbsub(int n, float *a, float *w); 321 void rftfsub(int n, float *a, int nc, float *c); 322 void rftbsub(int n, float *a, int nc, float *c); 323 int nw, nc; 324 float xi; 325 326 nw = ip[0]; 327 if (n > (nw << 2)) { 328 nw = n >> 2; 329 makewt(nw, ip, w); 330 } 331 nc = ip[1]; 332 if (n > (nc << 2)) { 333 nc = n >> 2; 334 makect(nc, ip, w + nw); 335 } 336 if (isgn >= 0) { 337 if (n > 4) { 338 bitrv2(n, ip + 2, a); 339 cftfsub(n, a, w); 340 rftfsub(n, a, nc, w + nw); 341 } else if (n == 4) { 342 cftfsub(n, a, w); 343 } 344 xi = a[0] - a[1]; 345 a[0] += a[1]; 346 a[1] = xi; 347 } else { 348 a[1] = 0.5f * (a[0] - a[1]); 349 a[0] -= a[1]; 350 if (n > 4) { 351 rftbsub(n, a, nc, w + nw); 352 bitrv2(n, ip + 2, a); 353 cftbsub(n, a, w); 354 } else if (n == 4) { 355 cftfsub(n, a, w); 356 } 357 } 358 } 359 360 361 void ddct(int n, int isgn, float *a, int *ip, float *w) 362 { 363 void makewt(int nw, int *ip, float *w); 364 void makect(int nc, int *ip, float *c); 365 void bitrv2(int n, int *ip, float *a); 366 void cftfsub(int n, float *a, float *w); 367 void cftbsub(int n, float *a, float *w); 368 void rftfsub(int n, float *a, int nc, float *c); 369 void rftbsub(int n, float *a, int nc, float *c); 370 void dctsub(int n, float *a, int nc, float *c); 371 int j, nw, nc; 372 float xr; 373 374 nw = ip[0]; 375 if (n > (nw << 2)) { 376 nw = n >> 2; 377 makewt(nw, ip, w); 378 } 379 nc = ip[1]; 380 if (n > nc) { 381 nc = n; 382 makect(nc, ip, w + nw); 383 } 384 if (isgn < 0) { 385 xr = a[n - 1]; 386 for (j = n - 2; j >= 2; j -= 2) { 387 a[j + 1] = a[j] - a[j - 1]; 388 a[j] += a[j - 1]; 389 } 390 a[1] = a[0] - xr; 391 a[0] += xr; 392 if (n > 4) { 393 rftbsub(n, a, nc, w + nw); 394 bitrv2(n, ip + 2, a); 395 cftbsub(n, a, w); 396 } else if (n == 4) { 397 cftfsub(n, a, w); 398 } 399 } 400 dctsub(n, a, nc, w + nw); 401 if (isgn >= 0) { 402 if (n > 4) { 403 bitrv2(n, ip + 2, a); 404 cftfsub(n, a, w); 405 rftfsub(n, a, nc, w + nw); 406 } else if (n == 4) { 407 cftfsub(n, a, w); 408 } 409 xr = a[0] - a[1]; 410 a[0] += a[1]; 411 for (j = 2; j < n; j += 2) { 412 a[j - 1] = a[j] - a[j + 1]; 413 a[j] += a[j + 1]; 414 } 415 a[n - 1] = xr; 416 } 417 } 418 419 420 void ddst(int n, int isgn, float *a, int *ip, float *w) 421 { 422 void makewt(int nw, int *ip, float *w); 423 void makect(int nc, int *ip, float *c); 424 void bitrv2(int n, int *ip, float *a); 425 void cftfsub(int n, float *a, float *w); 426 void cftbsub(int n, float *a, float *w); 427 void rftfsub(int n, float *a, int nc, float *c); 428 void rftbsub(int n, float *a, int nc, float *c); 429 void dstsub(int n, float *a, int nc, float *c); 430 int j, nw, nc; 431 float xr; 432 433 nw = ip[0]; 434 if (n > (nw << 2)) { 435 nw = n >> 2; 436 makewt(nw, ip, w); 437 } 438 nc = ip[1]; 439 if (n > nc) { 440 nc = n; 441 makect(nc, ip, w + nw); 442 } 443 if (isgn < 0) { 444 xr = a[n - 1]; 445 for (j = n - 2; j >= 2; j -= 2) { 446 a[j + 1] = -a[j] - a[j - 1]; 447 a[j] -= a[j - 1]; 448 } 449 a[1] = a[0] + xr; 450 a[0] -= xr; 451 if (n > 4) { 452 rftbsub(n, a, nc, w + nw); 453 bitrv2(n, ip + 2, a); 454 cftbsub(n, a, w); 455 } else if (n == 4) { 456 cftfsub(n, a, w); 457 } 458 } 459 dstsub(n, a, nc, w + nw); 460 if (isgn >= 0) { 461 if (n > 4) { 462 bitrv2(n, ip + 2, a); 463 cftfsub(n, a, w); 464 rftfsub(n, a, nc, w + nw); 465 } else if (n == 4) { 466 cftfsub(n, a, w); 467 } 468 xr = a[0] - a[1]; 469 a[0] += a[1]; 470 for (j = 2; j < n; j += 2) { 471 a[j - 1] = -a[j] - a[j + 1]; 472 a[j] -= a[j + 1]; 473 } 474 a[n - 1] = -xr; 475 } 476 } 477 478 479 void dfct(int n, float *a, float *t, int *ip, float *w) 480 { 481 void makewt(int nw, int *ip, float *w); 482 void makect(int nc, int *ip, float *c); 483 void bitrv2(int n, int *ip, float *a); 484 void cftfsub(int n, float *a, float *w); 485 void rftfsub(int n, float *a, int nc, float *c); 486 void dctsub(int n, float *a, int nc, float *c); 487 int j, k, l, m, mh, nw, nc; 488 float xr, xi, yr, yi; 489 490 nw = ip[0]; 491 if (n > (nw << 3)) { 492 nw = n >> 3; 493 makewt(nw, ip, w); 494 } 495 nc = ip[1]; 496 if (n > (nc << 1)) { 497 nc = n >> 1; 498 makect(nc, ip, w + nw); 499 } 500 m = n >> 1; 501 yi = a[m]; 502 xi = a[0] + a[n]; 503 a[0] -= a[n]; 504 t[0] = xi - yi; 505 t[m] = xi + yi; 506 if (n > 2) { 507 mh = m >> 1; 508 for (j = 1; j < mh; j++) { 509 k = m - j; 510 xr = a[j] - a[n - j]; 511 xi = a[j] + a[n - j]; 512 yr = a[k] - a[n - k]; 513 yi = a[k] + a[n - k]; 514 a[j] = xr; 515 a[k] = yr; 516 t[j] = xi - yi; 517 t[k] = xi + yi; 518 } 519 t[mh] = a[mh] + a[n - mh]; 520 a[mh] -= a[n - mh]; 521 dctsub(m, a, nc, w + nw); 522 if (m > 4) { 523 bitrv2(m, ip + 2, a); 524 cftfsub(m, a, w); 525 rftfsub(m, a, nc, w + nw); 526 } else if (m == 4) { 527 cftfsub(m, a, w); 528 } 529 a[n - 1] = a[0] - a[1]; 530 a[1] = a[0] + a[1]; 531 for (j = m - 2; j >= 2; j -= 2) { 532 a[2 * j + 1] = a[j] + a[j + 1]; 533 a[2 * j - 1] = a[j] - a[j + 1]; 534 } 535 l = 2; 536 m = mh; 537 while (m >= 2) { 538 dctsub(m, t, nc, w + nw); 539 if (m > 4) { 540 bitrv2(m, ip + 2, t); 541 cftfsub(m, t, w); 542 rftfsub(m, t, nc, w + nw); 543 } else if (m == 4) { 544 cftfsub(m, t, w); 545 } 546 a[n - l] = t[0] - t[1]; 547 a[l] = t[0] + t[1]; 548 k = 0; 549 for (j = 2; j < m; j += 2) { 550 k += l << 2; 551 a[k - l] = t[j] - t[j + 1]; 552 a[k + l] = t[j] + t[j + 1]; 553 } 554 l <<= 1; 555 mh = m >> 1; 556 for (j = 0; j < mh; j++) { 557 k = m - j; 558 t[j] = t[m + k] - t[m + j]; 559 t[k] = t[m + k] + t[m + j]; 560 } 561 t[mh] = t[m + mh]; 562 m = mh; 563 } 564 a[l] = t[0]; 565 a[n] = t[2] - t[1]; 566 a[0] = t[2] + t[1]; 567 } else { 568 a[1] = a[0]; 569 a[2] = t[0]; 570 a[0] = t[1]; 571 } 572 } 573 574 575 void dfst(int n, float *a, float *t, int *ip, float *w) 576 { 577 void makewt(int nw, int *ip, float *w); 578 void makect(int nc, int *ip, float *c); 579 void bitrv2(int n, int *ip, float *a); 580 void cftfsub(int n, float *a, float *w); 581 void rftfsub(int n, float *a, int nc, float *c); 582 void dstsub(int n, float *a, int nc, float *c); 583 int j, k, l, m, mh, nw, nc; 584 float xr, xi, yr, yi; 585 586 nw = ip[0]; 587 if (n > (nw << 3)) { 588 nw = n >> 3; 589 makewt(nw, ip, w); 590 } 591 nc = ip[1]; 592 if (n > (nc << 1)) { 593 nc = n >> 1; 594 makect(nc, ip, w + nw); 595 } 596 if (n > 2) { 597 m = n >> 1; 598 mh = m >> 1; 599 for (j = 1; j < mh; j++) { 600 k = m - j; 601 xr = a[j] + a[n - j]; 602 xi = a[j] - a[n - j]; 603 yr = a[k] + a[n - k]; 604 yi = a[k] - a[n - k]; 605 a[j] = xr; 606 a[k] = yr; 607 t[j] = xi + yi; 608 t[k] = xi - yi; 609 } 610 t[0] = a[mh] - a[n - mh]; 611 a[mh] += a[n - mh]; 612 a[0] = a[m]; 613 dstsub(m, a, nc, w + nw); 614 if (m > 4) { 615 bitrv2(m, ip + 2, a); 616 cftfsub(m, a, w); 617 rftfsub(m, a, nc, w + nw); 618 } else if (m == 4) { 619 cftfsub(m, a, w); 620 } 621 a[n - 1] = a[1] - a[0]; 622 a[1] = a[0] + a[1]; 623 for (j = m - 2; j >= 2; j -= 2) { 624 a[2 * j + 1] = a[j] - a[j + 1]; 625 a[2 * j - 1] = -a[j] - a[j + 1]; 626 } 627 l = 2; 628 m = mh; 629 while (m >= 2) { 630 dstsub(m, t, nc, w + nw); 631 if (m > 4) { 632 bitrv2(m, ip + 2, t); 633 cftfsub(m, t, w); 634 rftfsub(m, t, nc, w + nw); 635 } else if (m == 4) { 636 cftfsub(m, t, w); 637 } 638 a[n - l] = t[1] - t[0]; 639 a[l] = t[0] + t[1]; 640 k = 0; 641 for (j = 2; j < m; j += 2) { 642 k += l << 2; 643 a[k - l] = -t[j] - t[j + 1]; 644 a[k + l] = t[j] - t[j + 1]; 645 } 646 l <<= 1; 647 mh = m >> 1; 648 for (j = 1; j < mh; j++) { 649 k = m - j; 650 t[j] = t[m + k] + t[m + j]; 651 t[k] = t[m + k] - t[m + j]; 652 } 653 t[0] = t[m + mh]; 654 m = mh; 655 } 656 a[l] = t[0]; 657 } 658 a[0] = 0; 659 } 660 661 662 /* -------- initializing routines -------- */ 663 664 665 #include <math.h> 666 667 void makewt(int nw, int *ip, float *w) 668 { 669 void bitrv2(int n, int *ip, float *a); 670 int j, nwh; 671 float delta, x, y; 672 673 ip[0] = nw; 674 ip[1] = 1; 675 if (nw > 2) { 676 nwh = nw >> 1; 677 delta = (float)atan(1.0f) / nwh; 678 w[0] = 1; 679 w[1] = 0; 680 w[nwh] = (float)cos(delta * nwh); 681 w[nwh + 1] = w[nwh]; 682 if (nwh > 2) { 683 for (j = 2; j < nwh; j += 2) { 684 x = (float)cos(delta * j); 685 y = (float)sin(delta * j); 686 w[j] = x; 687 w[j + 1] = y; 688 w[nw - j] = y; 689 w[nw - j + 1] = x; 690 } 691 bitrv2(nw, ip + 2, w); 692 } 693 } 694 } 695 696 697 void makect(int nc, int *ip, float *c) 698 { 699 int j, nch; 700 float delta; 701 702 ip[1] = nc; 703 if (nc > 1) { 704 nch = nc >> 1; 705 delta = (float)atan(1.0f) / nch; 706 c[0] = (float)cos(delta * nch); 707 c[nch] = 0.5f * c[0]; 708 for (j = 1; j < nch; j++) { 709 c[j] = 0.5f * (float)cos(delta * j); 710 c[nc - j] = 0.5f * (float)sin(delta * j); 711 } 712 } 713 } 714 715 716 /* -------- child routines -------- */ 717 718 719 void bitrv2(int n, int *ip, float *a) 720 { 721 int j, j1, k, k1, l, m, m2; 722 float xr, xi, yr, yi; 723 724 ip[0] = 0; 725 l = n; 726 m = 1; 727 while ((m << 3) < l) { 728 l >>= 1; 729 for (j = 0; j < m; j++) { 730 ip[m + j] = ip[j] + l; 731 } 732 m <<= 1; 733 } 734 m2 = 2 * m; 735 if ((m << 3) == l) { 736 for (k = 0; k < m; k++) { 737 for (j = 0; j < k; j++) { 738 j1 = 2 * j + ip[k]; 739 k1 = 2 * k + ip[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 += m2; 749 k1 += 2 * m2; 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 += m2; 759 k1 -= m2; 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 += m2; 769 k1 += 2 * m2; 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 } 779 j1 = 2 * k + m2 + ip[k]; 780 k1 = j1 + m2; 781 xr = a[j1]; 782 xi = a[j1 + 1]; 783 yr = a[k1]; 784 yi = a[k1 + 1]; 785 a[j1] = yr; 786 a[j1 + 1] = yi; 787 a[k1] = xr; 788 a[k1 + 1] = xi; 789 } 790 } else { 791 for (k = 1; k < m; k++) { 792 for (j = 0; j < k; j++) { 793 j1 = 2 * j + ip[k]; 794 k1 = 2 * k + ip[j]; 795 xr = a[j1]; 796 xi = a[j1 + 1]; 797 yr = a[k1]; 798 yi = a[k1 + 1]; 799 a[j1] = yr; 800 a[j1 + 1] = yi; 801 a[k1] = xr; 802 a[k1 + 1] = xi; 803 j1 += m2; 804 k1 += m2; 805 xr = a[j1]; 806 xi = a[j1 + 1]; 807 yr = a[k1]; 808 yi = a[k1 + 1]; 809 a[j1] = yr; 810 a[j1 + 1] = yi; 811 a[k1] = xr; 812 a[k1 + 1] = xi; 813 } 814 } 815 } 816 } 817 818 819 void bitrv2conj(int n, int *ip, float *a) 820 { 821 int j, j1, k, k1, l, m, m2; 822 float xr, xi, yr, yi; 823 824 ip[0] = 0; 825 l = n; 826 m = 1; 827 while ((m << 3) < l) { 828 l >>= 1; 829 for (j = 0; j < m; j++) { 830 ip[m + j] = ip[j] + l; 831 } 832 m <<= 1; 833 } 834 m2 = 2 * m; 835 if ((m << 3) == l) { 836 for (k = 0; k < m; k++) { 837 for (j = 0; j < k; j++) { 838 j1 = 2 * j + ip[k]; 839 k1 = 2 * k + ip[j]; 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 += m2; 849 k1 += 2 * m2; 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 += m2; 859 k1 -= m2; 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 += m2; 869 k1 += 2 * m2; 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 } 879 k1 = 2 * k + ip[k]; 880 a[k1 + 1] = -a[k1 + 1]; 881 j1 = k1 + m2; 882 k1 = j1 + m2; 883 xr = a[j1]; 884 xi = -a[j1 + 1]; 885 yr = a[k1]; 886 yi = -a[k1 + 1]; 887 a[j1] = yr; 888 a[j1 + 1] = yi; 889 a[k1] = xr; 890 a[k1 + 1] = xi; 891 k1 += m2; 892 a[k1 + 1] = -a[k1 + 1]; 893 } 894 } else { 895 a[1] = -a[1]; 896 a[m2 + 1] = -a[m2 + 1]; 897 for (k = 1; k < m; k++) { 898 for (j = 0; j < k; j++) { 899 j1 = 2 * j + ip[k]; 900 k1 = 2 * k + ip[j]; 901 xr = a[j1]; 902 xi = -a[j1 + 1]; 903 yr = a[k1]; 904 yi = -a[k1 + 1]; 905 a[j1] = yr; 906 a[j1 + 1] = yi; 907 a[k1] = xr; 908 a[k1 + 1] = xi; 909 j1 += m2; 910 k1 += m2; 911 xr = a[j1]; 912 xi = -a[j1 + 1]; 913 yr = a[k1]; 914 yi = -a[k1 + 1]; 915 a[j1] = yr; 916 a[j1 + 1] = yi; 917 a[k1] = xr; 918 a[k1 + 1] = xi; 919 } 920 k1 = 2 * k + ip[k]; 921 a[k1 + 1] = -a[k1 + 1]; 922 a[k1 + m2 + 1] = -a[k1 + m2 + 1]; 923 } 924 } 925 } 926 927 928 void cftfsub(int n, float *a, float *w) 929 { 930 void cft1st(int n, float *a, float *w); 931 void cftmdl(int n, int l, float *a, float *w); 932 int j, j1, j2, j3, l; 933 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 934 935 l = 2; 936 if (n > 8) { 937 cft1st(n, a, w); 938 l = 8; 939 while ((l << 2) < n) { 940 cftmdl(n, l, a, w); 941 l <<= 2; 942 } 943 } 944 if ((l << 2) == n) { 945 for (j = 0; j < l; j += 2) { 946 j1 = j + l; 947 j2 = j1 + l; 948 j3 = j2 + l; 949 x0r = a[j] + a[j1]; 950 x0i = a[j + 1] + a[j1 + 1]; 951 x1r = a[j] - a[j1]; 952 x1i = a[j + 1] - a[j1 + 1]; 953 x2r = a[j2] + a[j3]; 954 x2i = a[j2 + 1] + a[j3 + 1]; 955 x3r = a[j2] - a[j3]; 956 x3i = a[j2 + 1] - a[j3 + 1]; 957 a[j] = x0r + x2r; 958 a[j + 1] = x0i + x2i; 959 a[j2] = x0r - x2r; 960 a[j2 + 1] = x0i - x2i; 961 a[j1] = x1r - x3i; 962 a[j1 + 1] = x1i + x3r; 963 a[j3] = x1r + x3i; 964 a[j3 + 1] = x1i - x3r; 965 } 966 } else { 967 for (j = 0; j < l; j += 2) { 968 j1 = j + l; 969 x0r = a[j] - a[j1]; 970 x0i = a[j + 1] - a[j1 + 1]; 971 a[j] += a[j1]; 972 a[j + 1] += a[j1 + 1]; 973 a[j1] = x0r; 974 a[j1 + 1] = x0i; 975 } 976 } 977 } 978 979 980 void cftbsub(int n, float *a, float *w) 981 { 982 void cft1st(int n, float *a, float *w); 983 void cftmdl(int n, int l, float *a, float *w); 984 int j, j1, j2, j3, l; 985 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 986 987 l = 2; 988 if (n > 8) { 989 cft1st(n, a, w); 990 l = 8; 991 while ((l << 2) < n) { 992 cftmdl(n, l, a, w); 993 l <<= 2; 994 } 995 } 996 if ((l << 2) == n) { 997 for (j = 0; j < l; j += 2) { 998 j1 = j + l; 999 j2 = j1 + l; 1000 j3 = j2 + l; 1001 x0r = a[j] + a[j1]; 1002 x0i = -a[j + 1] - a[j1 + 1]; 1003 x1r = a[j] - a[j1]; 1004 x1i = -a[j + 1] + a[j1 + 1]; 1005 x2r = a[j2] + a[j3]; 1006 x2i = a[j2 + 1] + a[j3 + 1]; 1007 x3r = a[j2] - a[j3]; 1008 x3i = a[j2 + 1] - a[j3 + 1]; 1009 a[j] = x0r + x2r; 1010 a[j + 1] = x0i - x2i; 1011 a[j2] = x0r - x2r; 1012 a[j2 + 1] = x0i + x2i; 1013 a[j1] = x1r - x3i; 1014 a[j1 + 1] = x1i - x3r; 1015 a[j3] = x1r + x3i; 1016 a[j3 + 1] = x1i + x3r; 1017 } 1018 } else { 1019 for (j = 0; j < l; j += 2) { 1020 j1 = j + l; 1021 x0r = a[j] - a[j1]; 1022 x0i = -a[j + 1] + a[j1 + 1]; 1023 a[j] += a[j1]; 1024 a[j + 1] = -a[j + 1] - a[j1 + 1]; 1025 a[j1] = x0r; 1026 a[j1 + 1] = x0i; 1027 } 1028 } 1029 } 1030 1031 1032 void cft1st(int n, float *a, float *w) 1033 { 1034 int j, k1, k2; 1035 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; 1036 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 1037 1038 x0r = a[0] + a[2]; 1039 x0i = a[1] + a[3]; 1040 x1r = a[0] - a[2]; 1041 x1i = a[1] - a[3]; 1042 x2r = a[4] + a[6]; 1043 x2i = a[5] + a[7]; 1044 x3r = a[4] - a[6]; 1045 x3i = a[5] - a[7]; 1046 a[0] = x0r + x2r; 1047 a[1] = x0i + x2i; 1048 a[4] = x0r - x2r; 1049 a[5] = x0i - x2i; 1050 a[2] = x1r - x3i; 1051 a[3] = x1i + x3r; 1052 a[6] = x1r + x3i; 1053 a[7] = x1i - x3r; 1054 wk1r = w[2]; 1055 x0r = a[8] + a[10]; 1056 x0i = a[9] + a[11]; 1057 x1r = a[8] - a[10]; 1058 x1i = a[9] - a[11]; 1059 x2r = a[12] + a[14]; 1060 x2i = a[13] + a[15]; 1061 x3r = a[12] - a[14]; 1062 x3i = a[13] - a[15]; 1063 a[8] = x0r + x2r; 1064 a[9] = x0i + x2i; 1065 a[12] = x2i - x0i; 1066 a[13] = x0r - x2r; 1067 x0r = x1r - x3i; 1068 x0i = x1i + x3r; 1069 a[10] = wk1r * (x0r - x0i); 1070 a[11] = wk1r * (x0r + x0i); 1071 x0r = x3i + x1r; 1072 x0i = x3r - x1i; 1073 a[14] = wk1r * (x0i - x0r); 1074 a[15] = wk1r * (x0i + x0r); 1075 k1 = 0; 1076 for (j = 16; j < n; j += 16) { 1077 k1 += 2; 1078 k2 = 2 * k1; 1079 wk2r = w[k1]; 1080 wk2i = w[k1 + 1]; 1081 wk1r = w[k2]; 1082 wk1i = w[k2 + 1]; 1083 wk3r = wk1r - 2 * wk2i * wk1i; 1084 wk3i = 2 * wk2i * wk1r - wk1i; 1085 x0r = a[j] + a[j + 2]; 1086 x0i = a[j + 1] + a[j + 3]; 1087 x1r = a[j] - a[j + 2]; 1088 x1i = a[j + 1] - a[j + 3]; 1089 x2r = a[j + 4] + a[j + 6]; 1090 x2i = a[j + 5] + a[j + 7]; 1091 x3r = a[j + 4] - a[j + 6]; 1092 x3i = a[j + 5] - a[j + 7]; 1093 a[j] = x0r + x2r; 1094 a[j + 1] = x0i + x2i; 1095 x0r -= x2r; 1096 x0i -= x2i; 1097 a[j + 4] = wk2r * x0r - wk2i * x0i; 1098 a[j + 5] = wk2r * x0i + wk2i * x0r; 1099 x0r = x1r - x3i; 1100 x0i = x1i + x3r; 1101 a[j + 2] = wk1r * x0r - wk1i * x0i; 1102 a[j + 3] = wk1r * x0i + wk1i * x0r; 1103 x0r = x1r + x3i; 1104 x0i = x1i - x3r; 1105 a[j + 6] = wk3r * x0r - wk3i * x0i; 1106 a[j + 7] = wk3r * x0i + wk3i * x0r; 1107 wk1r = w[k2 + 2]; 1108 wk1i = w[k2 + 3]; 1109 wk3r = wk1r - 2 * wk2r * wk1i; 1110 wk3i = 2 * wk2r * wk1r - wk1i; 1111 x0r = a[j + 8] + a[j + 10]; 1112 x0i = a[j + 9] + a[j + 11]; 1113 x1r = a[j + 8] - a[j + 10]; 1114 x1i = a[j + 9] - a[j + 11]; 1115 x2r = a[j + 12] + a[j + 14]; 1116 x2i = a[j + 13] + a[j + 15]; 1117 x3r = a[j + 12] - a[j + 14]; 1118 x3i = a[j + 13] - a[j + 15]; 1119 a[j + 8] = x0r + x2r; 1120 a[j + 9] = x0i + x2i; 1121 x0r -= x2r; 1122 x0i -= x2i; 1123 a[j + 12] = -wk2i * x0r - wk2r * x0i; 1124 a[j + 13] = -wk2i * x0i + wk2r * x0r; 1125 x0r = x1r - x3i; 1126 x0i = x1i + x3r; 1127 a[j + 10] = wk1r * x0r - wk1i * x0i; 1128 a[j + 11] = wk1r * x0i + wk1i * x0r; 1129 x0r = x1r + x3i; 1130 x0i = x1i - x3r; 1131 a[j + 14] = wk3r * x0r - wk3i * x0i; 1132 a[j + 15] = wk3r * x0i + wk3i * x0r; 1133 } 1134 } 1135 1136 1137 void cftmdl(int n, int l, float *a, float *w) 1138 { 1139 int j, j1, j2, j3, k, k1, k2, m, m2; 1140 float wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; 1141 float x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 1142 1143 m = l << 2; 1144 for (j = 0; j < l; j += 2) { 1145 j1 = j + l; 1146 j2 = j1 + l; 1147 j3 = j2 + l; 1148 x0r = a[j] + a[j1]; 1149 x0i = a[j + 1] + a[j1 + 1]; 1150 x1r = a[j] - a[j1]; 1151 x1i = a[j + 1] - a[j1 + 1]; 1152 x2r = a[j2] + a[j3]; 1153 x2i = a[j2 + 1] + a[j3 + 1]; 1154 x3r = a[j2] - a[j3]; 1155 x3i = a[j2 + 1] - a[j3 + 1]; 1156 a[j] = x0r + x2r; 1157 a[j + 1] = x0i + x2i; 1158 a[j2] = x0r - x2r; 1159 a[j2 + 1] = x0i - x2i; 1160 a[j1] = x1r - x3i; 1161 a[j1 + 1] = x1i + x3r; 1162 a[j3] = x1r + x3i; 1163 a[j3 + 1] = x1i - x3r; 1164 } 1165 wk1r = w[2]; 1166 for (j = m; j < l + m; j += 2) { 1167 j1 = j + l; 1168 j2 = j1 + l; 1169 j3 = j2 + l; 1170 x0r = a[j] + a[j1]; 1171 x0i = a[j + 1] + a[j1 + 1]; 1172 x1r = a[j] - a[j1]; 1173 x1i = a[j + 1] - a[j1 + 1]; 1174 x2r = a[j2] + a[j3]; 1175 x2i = a[j2 + 1] + a[j3 + 1]; 1176 x3r = a[j2] - a[j3]; 1177 x3i = a[j2 + 1] - a[j3 + 1]; 1178 a[j] = x0r + x2r; 1179 a[j + 1] = x0i + x2i; 1180 a[j2] = x2i - x0i; 1181 a[j2 + 1] = x0r - x2r; 1182 x0r = x1r - x3i; 1183 x0i = x1i + x3r; 1184 a[j1] = wk1r * (x0r - x0i); 1185 a[j1 + 1] = wk1r * (x0r + x0i); 1186 x0r = x3i + x1r; 1187 x0i = x3r - x1i; 1188 a[j3] = wk1r * (x0i - x0r); 1189 a[j3 + 1] = wk1r * (x0i + x0r); 1190 } 1191 k1 = 0; 1192 m2 = 2 * m; 1193 for (k = m2; k < n; k += m2) { 1194 k1 += 2; 1195 k2 = 2 * k1; 1196 wk2r = w[k1]; 1197 wk2i = w[k1 + 1]; 1198 wk1r = w[k2]; 1199 wk1i = w[k2 + 1]; 1200 wk3r = wk1r - 2 * wk2i * wk1i; 1201 wk3i = 2 * wk2i * wk1r - wk1i; 1202 for (j = k; j < l + k; j += 2) { 1203 j1 = j + l; 1204 j2 = j1 + l; 1205 j3 = j2 + l; 1206 x0r = a[j] + a[j1]; 1207 x0i = a[j + 1] + a[j1 + 1]; 1208 x1r = a[j] - a[j1]; 1209 x1i = a[j + 1] - a[j1 + 1]; 1210 x2r = a[j2] + a[j3]; 1211 x2i = a[j2 + 1] + a[j3 + 1]; 1212 x3r = a[j2] - a[j3]; 1213 x3i = a[j2 + 1] - a[j3 + 1]; 1214 a[j] = x0r + x2r; 1215 a[j + 1] = x0i + x2i; 1216 x0r -= x2r; 1217 x0i -= x2i; 1218 a[j2] = wk2r * x0r - wk2i * x0i; 1219 a[j2 + 1] = wk2r * x0i + wk2i * x0r; 1220 x0r = x1r - x3i; 1221 x0i = x1i + x3r; 1222 a[j1] = wk1r * x0r - wk1i * x0i; 1223 a[j1 + 1] = wk1r * x0i + wk1i * x0r; 1224 x0r = x1r + x3i; 1225 x0i = x1i - x3r; 1226 a[j3] = wk3r * x0r - wk3i * x0i; 1227 a[j3 + 1] = wk3r * x0i + wk3i * x0r; 1228 } 1229 wk1r = w[k2 + 2]; 1230 wk1i = w[k2 + 3]; 1231 wk3r = wk1r - 2 * wk2r * wk1i; 1232 wk3i = 2 * wk2r * wk1r - wk1i; 1233 for (j = k + m; j < l + (k + m); j += 2) { 1234 j1 = j + l; 1235 j2 = j1 + l; 1236 j3 = j2 + l; 1237 x0r = a[j] + a[j1]; 1238 x0i = a[j + 1] + a[j1 + 1]; 1239 x1r = a[j] - a[j1]; 1240 x1i = a[j + 1] - a[j1 + 1]; 1241 x2r = a[j2] + a[j3]; 1242 x2i = a[j2 + 1] + a[j3 + 1]; 1243 x3r = a[j2] - a[j3]; 1244 x3i = a[j2 + 1] - a[j3 + 1]; 1245 a[j] = x0r + x2r; 1246 a[j + 1] = x0i + x2i; 1247 x0r -= x2r; 1248 x0i -= x2i; 1249 a[j2] = -wk2i * x0r - wk2r * x0i; 1250 a[j2 + 1] = -wk2i * x0i + wk2r * x0r; 1251 x0r = x1r - x3i; 1252 x0i = x1i + x3r; 1253 a[j1] = wk1r * x0r - wk1i * x0i; 1254 a[j1 + 1] = wk1r * x0i + wk1i * x0r; 1255 x0r = x1r + x3i; 1256 x0i = x1i - x3r; 1257 a[j3] = wk3r * x0r - wk3i * x0i; 1258 a[j3 + 1] = wk3r * x0i + wk3i * x0r; 1259 } 1260 } 1261 } 1262 1263 1264 void rftfsub(int n, float *a, int nc, float *c) 1265 { 1266 int j, k, kk, ks, m; 1267 float wkr, wki, xr, xi, yr, yi; 1268 1269 m = n >> 1; 1270 ks = 2 * nc / m; 1271 kk = 0; 1272 for (j = 2; j < m; j += 2) { 1273 k = n - j; 1274 kk += ks; 1275 wkr = 0.5f - c[nc - kk]; 1276 wki = c[kk]; 1277 xr = a[j] - a[k]; 1278 xi = a[j + 1] + a[k + 1]; 1279 yr = wkr * xr - wki * xi; 1280 yi = wkr * xi + wki * xr; 1281 a[j] -= yr; 1282 a[j + 1] -= yi; 1283 a[k] += yr; 1284 a[k + 1] -= yi; 1285 } 1286 } 1287 1288 1289 void rftbsub(int n, float *a, int nc, float *c) 1290 { 1291 int j, k, kk, ks, m; 1292 float wkr, wki, xr, xi, yr, yi; 1293 1294 a[1] = -a[1]; 1295 m = n >> 1; 1296 ks = 2 * nc / m; 1297 kk = 0; 1298 for (j = 2; j < m; j += 2) { 1299 k = n - j; 1300 kk += ks; 1301 wkr = 0.5f - c[nc - kk]; 1302 wki = c[kk]; 1303 xr = a[j] - a[k]; 1304 xi = a[j + 1] + a[k + 1]; 1305 yr = wkr * xr + wki * xi; 1306 yi = wkr * xi - wki * xr; 1307 a[j] -= yr; 1308 a[j + 1] = yi - a[j + 1]; 1309 a[k] += yr; 1310 a[k + 1] = yi - a[k + 1]; 1311 } 1312 a[m + 1] = -a[m + 1]; 1313 } 1314 1315 1316 void dctsub(int n, float *a, int nc, float *c) 1317 { 1318 int j, k, kk, ks, m; 1319 float wkr, wki, xr; 1320 1321 m = n >> 1; 1322 ks = nc / n; 1323 kk = 0; 1324 for (j = 1; j < m; j++) { 1325 k = n - j; 1326 kk += ks; 1327 wkr = c[kk] - c[nc - kk]; 1328 wki = c[kk] + c[nc - kk]; 1329 xr = wki * a[j] - wkr * a[k]; 1330 a[j] = wkr * a[j] + wki * a[k]; 1331 a[k] = xr; 1332 } 1333 a[m] *= c[0]; 1334 } 1335 1336 1337 void dstsub(int n, float *a, int nc, float *c) 1338 { 1339 int j, k, kk, ks, m; 1340 float wkr, wki, xr; 1341 1342 m = n >> 1; 1343 ks = nc / n; 1344 kk = 0; 1345 for (j = 1; j < m; j++) { 1346 k = n - j; 1347 kk += ks; 1348 wkr = c[kk] - c[nc - kk]; 1349 wki = c[kk] + c[nc - kk]; 1350 xr = wki * a[k] - wkr * a[j]; 1351 a[k] = wkr * a[k] + wki * a[j]; 1352 a[j] = xr; 1353 } 1354 a[m] *= c[0]; 1355 } 1356 1357
