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