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