001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 package org.apache.commons.math.transform; 018 019 import java.io.Serializable; 020 import java.lang.reflect.Array; 021 022 import org.apache.commons.math.FunctionEvaluationException; 023 import org.apache.commons.math.MathRuntimeException; 024 import org.apache.commons.math.analysis.UnivariateRealFunction; 025 import org.apache.commons.math.complex.Complex; 026 027 /** 028 * Implements the <a href="http://mathworld.wolfram.com/FastFourierTransform.html"> 029 * Fast Fourier Transform</a> for transformation of one-dimensional data sets. 030 * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897, 031 * chapter 6. 032 * <p> 033 * There are several conventions for the definition of FFT and inverse FFT, 034 * mainly on different coefficient and exponent. Here the equations are listed 035 * in the comments of the corresponding methods.</p> 036 * <p> 037 * We require the length of data set to be power of 2, this greatly simplifies 038 * and speeds up the code. Users can pad the data with zeros to meet this 039 * requirement. There are other flavors of FFT, for reference, see S. Winograd, 040 * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation, 041 * 32 (1978), 175 - 199.</p> 042 * 043 * @version $Revision: 790243 $ $Date: 2009-07-01 12:03:28 -0400 (Wed, 01 Jul 2009) $ 044 * @since 1.2 045 */ 046 public class FastFourierTransformer implements Serializable { 047 048 /** Serializable version identifier. */ 049 static final long serialVersionUID = 5138259215438106000L; 050 051 /** roots of unity */ 052 private RootsOfUnity roots = new RootsOfUnity(); 053 054 /** 055 * Construct a default transformer. 056 */ 057 public FastFourierTransformer() { 058 super(); 059 } 060 061 /** 062 * Transform the given real data set. 063 * <p> 064 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $ 065 * </p> 066 * 067 * @param f the real data array to be transformed 068 * @return the complex transformed array 069 * @throws IllegalArgumentException if any parameters are invalid 070 */ 071 public Complex[] transform(double f[]) 072 throws IllegalArgumentException { 073 return fft(f, false); 074 } 075 076 /** 077 * Transform the given real function, sampled on the given interval. 078 * <p> 079 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $ 080 * </p> 081 * 082 * @param f the function to be sampled and transformed 083 * @param min the lower bound for the interval 084 * @param max the upper bound for the interval 085 * @param n the number of sample points 086 * @return the complex transformed array 087 * @throws FunctionEvaluationException if function cannot be evaluated 088 * at some point 089 * @throws IllegalArgumentException if any parameters are invalid 090 */ 091 public Complex[] transform(UnivariateRealFunction f, 092 double min, double max, int n) 093 throws FunctionEvaluationException, IllegalArgumentException { 094 double data[] = sample(f, min, max, n); 095 return fft(data, false); 096 } 097 098 /** 099 * Transform the given complex data set. 100 * <p> 101 * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $ 102 * </p> 103 * 104 * @param f the complex data array to be transformed 105 * @return the complex transformed array 106 * @throws IllegalArgumentException if any parameters are invalid 107 */ 108 public Complex[] transform(Complex f[]) 109 throws IllegalArgumentException { 110 roots.computeOmega(f.length); 111 return fft(f); 112 } 113 114 /** 115 * Transform the given real data set. 116 * <p> 117 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$ 118 * </p> 119 * 120 * @param f the real data array to be transformed 121 * @return the complex transformed array 122 * @throws IllegalArgumentException if any parameters are invalid 123 */ 124 public Complex[] transform2(double f[]) 125 throws IllegalArgumentException { 126 127 double scaling_coefficient = 1.0 / Math.sqrt(f.length); 128 return scaleArray(fft(f, false), scaling_coefficient); 129 } 130 131 /** 132 * Transform the given real function, sampled on the given interval. 133 * <p> 134 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$ 135 * </p> 136 * 137 * @param f the function to be sampled and transformed 138 * @param min the lower bound for the interval 139 * @param max the upper bound for the interval 140 * @param n the number of sample points 141 * @return the complex transformed array 142 * @throws FunctionEvaluationException if function cannot be evaluated 143 * at some point 144 * @throws IllegalArgumentException if any parameters are invalid 145 */ 146 public Complex[] transform2(UnivariateRealFunction f, 147 double min, double max, int n) 148 throws FunctionEvaluationException, IllegalArgumentException { 149 150 double data[] = sample(f, min, max, n); 151 double scaling_coefficient = 1.0 / Math.sqrt(n); 152 return scaleArray(fft(data, false), scaling_coefficient); 153 } 154 155 /** 156 * Transform the given complex data set. 157 * <p> 158 * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$ 159 * </p> 160 * 161 * @param f the complex data array to be transformed 162 * @return the complex transformed array 163 * @throws IllegalArgumentException if any parameters are invalid 164 */ 165 public Complex[] transform2(Complex f[]) 166 throws IllegalArgumentException { 167 168 roots.computeOmega(f.length); 169 double scaling_coefficient = 1.0 / Math.sqrt(f.length); 170 return scaleArray(fft(f), scaling_coefficient); 171 } 172 173 /** 174 * Inversely transform the given real data set. 175 * <p> 176 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $ 177 * </p> 178 * 179 * @param f the real data array to be inversely transformed 180 * @return the complex inversely transformed array 181 * @throws IllegalArgumentException if any parameters are invalid 182 */ 183 public Complex[] inversetransform(double f[]) 184 throws IllegalArgumentException { 185 186 double scaling_coefficient = 1.0 / f.length; 187 return scaleArray(fft(f, true), scaling_coefficient); 188 } 189 190 /** 191 * Inversely transform the given real function, sampled on the given interval. 192 * <p> 193 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $ 194 * </p> 195 * 196 * @param f the function to be sampled and inversely transformed 197 * @param min the lower bound for the interval 198 * @param max the upper bound for the interval 199 * @param n the number of sample points 200 * @return the complex inversely transformed array 201 * @throws FunctionEvaluationException if function cannot be evaluated 202 * at some point 203 * @throws IllegalArgumentException if any parameters are invalid 204 */ 205 public Complex[] inversetransform(UnivariateRealFunction f, 206 double min, double max, int n) 207 throws FunctionEvaluationException, IllegalArgumentException { 208 209 double data[] = sample(f, min, max, n); 210 double scaling_coefficient = 1.0 / n; 211 return scaleArray(fft(data, true), scaling_coefficient); 212 } 213 214 /** 215 * Inversely transform the given complex data set. 216 * <p> 217 * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $ 218 * </p> 219 * 220 * @param f the complex data array to be inversely transformed 221 * @return the complex inversely transformed array 222 * @throws IllegalArgumentException if any parameters are invalid 223 */ 224 public Complex[] inversetransform(Complex f[]) 225 throws IllegalArgumentException { 226 227 roots.computeOmega(-f.length); // pass negative argument 228 double scaling_coefficient = 1.0 / f.length; 229 return scaleArray(fft(f), scaling_coefficient); 230 } 231 232 /** 233 * Inversely transform the given real data set. 234 * <p> 235 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$ 236 * </p> 237 * 238 * @param f the real data array to be inversely transformed 239 * @return the complex inversely transformed array 240 * @throws IllegalArgumentException if any parameters are invalid 241 */ 242 public Complex[] inversetransform2(double f[]) 243 throws IllegalArgumentException { 244 245 double scaling_coefficient = 1.0 / Math.sqrt(f.length); 246 return scaleArray(fft(f, true), scaling_coefficient); 247 } 248 249 /** 250 * Inversely transform the given real function, sampled on the given interval. 251 * <p> 252 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$ 253 * </p> 254 * 255 * @param f the function to be sampled and inversely transformed 256 * @param min the lower bound for the interval 257 * @param max the upper bound for the interval 258 * @param n the number of sample points 259 * @return the complex inversely transformed array 260 * @throws FunctionEvaluationException if function cannot be evaluated 261 * at some point 262 * @throws IllegalArgumentException if any parameters are invalid 263 */ 264 public Complex[] inversetransform2(UnivariateRealFunction f, 265 double min, double max, int n) 266 throws FunctionEvaluationException, IllegalArgumentException { 267 268 double data[] = sample(f, min, max, n); 269 double scaling_coefficient = 1.0 / Math.sqrt(n); 270 return scaleArray(fft(data, true), scaling_coefficient); 271 } 272 273 /** 274 * Inversely transform the given complex data set. 275 * <p> 276 * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$ 277 * </p> 278 * 279 * @param f the complex data array to be inversely transformed 280 * @return the complex inversely transformed array 281 * @throws IllegalArgumentException if any parameters are invalid 282 */ 283 public Complex[] inversetransform2(Complex f[]) 284 throws IllegalArgumentException { 285 286 roots.computeOmega(-f.length); // pass negative argument 287 double scaling_coefficient = 1.0 / Math.sqrt(f.length); 288 return scaleArray(fft(f), scaling_coefficient); 289 } 290 291 /** 292 * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse). 293 * 294 * @param f the real data array to be transformed 295 * @param isInverse the indicator of forward or inverse transform 296 * @return the complex transformed array 297 * @throws IllegalArgumentException if any parameters are invalid 298 */ 299 protected Complex[] fft(double f[], boolean isInverse) 300 throws IllegalArgumentException { 301 302 verifyDataSet(f); 303 Complex F[] = new Complex[f.length]; 304 if (f.length == 1) { 305 F[0] = new Complex(f[0], 0.0); 306 return F; 307 } 308 309 // Rather than the naive real to complex conversion, pack 2N 310 // real numbers into N complex numbers for better performance. 311 int N = f.length >> 1; 312 Complex c[] = new Complex[N]; 313 for (int i = 0; i < N; i++) { 314 c[i] = new Complex(f[2*i], f[2*i+1]); 315 } 316 roots.computeOmega(isInverse ? -N : N); 317 Complex z[] = fft(c); 318 319 // reconstruct the FFT result for the original array 320 roots.computeOmega(isInverse ? -2*N : 2*N); 321 F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0); 322 F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0); 323 for (int i = 1; i < N; i++) { 324 Complex A = z[N-i].conjugate(); 325 Complex B = z[i].add(A); 326 Complex C = z[i].subtract(A); 327 //Complex D = roots.getOmega(i).multiply(Complex.I); 328 Complex D = new Complex(-roots.getOmegaImaginary(i), 329 roots.getOmegaReal(i)); 330 F[i] = B.subtract(C.multiply(D)); 331 F[2*N-i] = F[i].conjugate(); 332 } 333 334 return scaleArray(F, 0.5); 335 } 336 337 /** 338 * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse). 339 * 340 * @param data the complex data array to be transformed 341 * @return the complex transformed array 342 * @throws IllegalArgumentException if any parameters are invalid 343 */ 344 protected Complex[] fft(Complex data[]) 345 throws IllegalArgumentException { 346 347 int i, j, k, m, N = data.length; 348 Complex A, B, C, D, E, F, z, f[] = new Complex[N]; 349 350 // initial simple cases 351 verifyDataSet(data); 352 if (N == 1) { 353 f[0] = data[0]; 354 return f; 355 } 356 if (N == 2) { 357 f[0] = data[0].add(data[1]); 358 f[1] = data[0].subtract(data[1]); 359 return f; 360 } 361 362 // permute original data array in bit-reversal order 363 j = 0; 364 for (i = 0; i < N; i++) { 365 f[i] = data[j]; 366 k = N >> 1; 367 while (j >= k && k > 0) { 368 j -= k; k >>= 1; 369 } 370 j += k; 371 } 372 373 // the bottom base-4 round 374 for (i = 0; i < N; i += 4) { 375 A = f[i].add(f[i+1]); 376 B = f[i+2].add(f[i+3]); 377 C = f[i].subtract(f[i+1]); 378 D = f[i+2].subtract(f[i+3]); 379 E = C.add(D.multiply(Complex.I)); 380 F = C.subtract(D.multiply(Complex.I)); 381 f[i] = A.add(B); 382 f[i+2] = A.subtract(B); 383 // omegaCount indicates forward or inverse transform 384 f[i+1] = roots.isForward() ? F : E; 385 f[i+3] = roots.isForward() ? E : F; 386 } 387 388 // iterations from bottom to top take O(N*logN) time 389 for (i = 4; i < N; i <<= 1) { 390 m = N / (i<<1); 391 for (j = 0; j < N; j += i<<1) { 392 for (k = 0; k < i; k++) { 393 //z = f[i+j+k].multiply(roots.getOmega(k*m)); 394 final int k_times_m = k*m; 395 final double omega_k_times_m_real = roots.getOmegaReal(k_times_m); 396 final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m); 397 //z = f[i+j+k].multiply(omega[k*m]); 398 z = new Complex( 399 f[i+j+k].getReal() * omega_k_times_m_real - 400 f[i+j+k].getImaginary() * omega_k_times_m_imaginary, 401 f[i+j+k].getReal() * omega_k_times_m_imaginary + 402 f[i+j+k].getImaginary() * omega_k_times_m_real); 403 404 f[i+j+k] = f[j+k].subtract(z); 405 f[j+k] = f[j+k].add(z); 406 } 407 } 408 } 409 return f; 410 } 411 412 /** 413 * Sample the given univariate real function on the given interval. 414 * <p> 415 * The interval is divided equally into N sections and sample points 416 * are taken from min to max-(max-min)/N. Usually f(x) is periodic 417 * such that f(min) = f(max) (note max is not sampled), but we don't 418 * require that.</p> 419 * 420 * @param f the function to be sampled 421 * @param min the lower bound for the interval 422 * @param max the upper bound for the interval 423 * @param n the number of sample points 424 * @return the samples array 425 * @throws FunctionEvaluationException if function cannot be evaluated 426 * at some point 427 * @throws IllegalArgumentException if any parameters are invalid 428 */ 429 public static double[] sample(UnivariateRealFunction f, 430 double min, double max, int n) 431 throws FunctionEvaluationException, IllegalArgumentException { 432 433 if (n <= 0) { 434 throw MathRuntimeException.createIllegalArgumentException( 435 "number of sample is not positive: {0}", 436 n); 437 } 438 verifyInterval(min, max); 439 440 double s[] = new double[n]; 441 double h = (max - min) / n; 442 for (int i = 0; i < n; i++) { 443 s[i] = f.value(min + i * h); 444 } 445 return s; 446 } 447 448 /** 449 * Multiply every component in the given real array by the 450 * given real number. The change is made in place. 451 * 452 * @param f the real array to be scaled 453 * @param d the real scaling coefficient 454 * @return a reference to the scaled array 455 */ 456 public static double[] scaleArray(double f[], double d) { 457 for (int i = 0; i < f.length; i++) { 458 f[i] *= d; 459 } 460 return f; 461 } 462 463 /** 464 * Multiply every component in the given complex array by the 465 * given real number. The change is made in place. 466 * 467 * @param f the complex array to be scaled 468 * @param d the real scaling coefficient 469 * @return a reference to the scaled array 470 */ 471 public static Complex[] scaleArray(Complex f[], double d) { 472 for (int i = 0; i < f.length; i++) { 473 f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary()); 474 } 475 return f; 476 } 477 478 /** 479 * Returns true if the argument is power of 2. 480 * 481 * @param n the number to test 482 * @return true if the argument is power of 2 483 */ 484 public static boolean isPowerOf2(long n) { 485 return (n > 0) && ((n & (n - 1)) == 0); 486 } 487 488 /** 489 * Verifies that the data set has length of power of 2. 490 * 491 * @param d the data array 492 * @throws IllegalArgumentException if array length is not power of 2 493 */ 494 public static void verifyDataSet(double d[]) throws IllegalArgumentException { 495 if (!isPowerOf2(d.length)) { 496 throw MathRuntimeException.createIllegalArgumentException( 497 "{0} is not a power of 2, consider padding for fix", 498 d.length); 499 } 500 } 501 502 /** 503 * Verifies that the data set has length of power of 2. 504 * 505 * @param o the data array 506 * @throws IllegalArgumentException if array length is not power of 2 507 */ 508 public static void verifyDataSet(Object o[]) throws IllegalArgumentException { 509 if (!isPowerOf2(o.length)) { 510 throw MathRuntimeException.createIllegalArgumentException( 511 "{0} is not a power of 2, consider padding for fix", 512 o.length); 513 } 514 } 515 516 /** 517 * Verifies that the endpoints specify an interval. 518 * 519 * @param lower lower endpoint 520 * @param upper upper endpoint 521 * @throws IllegalArgumentException if not interval 522 */ 523 public static void verifyInterval(double lower, double upper) 524 throws IllegalArgumentException { 525 526 if (lower >= upper) { 527 throw MathRuntimeException.createIllegalArgumentException( 528 "endpoints do not specify an interval: [{0}, {1}]", 529 lower, upper); 530 } 531 } 532 533 /** 534 * Performs a multi-dimensional Fourier transform on a given array. 535 * Use {@link #inversetransform2(Complex[])} and 536 * {@link #transform2(Complex[])} in a row-column implementation 537 * in any number of dimensions with O(N×log(N)) complexity with 538 * N=n<sub>1</sub>×n<sub>2</sub>×n<sub>3</sub>×...×n<sub>d</sub>, 539 * n<sub>x</sub>=number of elements in dimension x, 540 * and d=total number of dimensions. 541 * 542 * @param mdca Multi-Dimensional Complex Array id est Complex[][][][] 543 * @param forward inverseTransform2 is preformed if this is false 544 * @return transform of mdca as a Multi-Dimensional Complex Array id est Complex[][][][] 545 * @throws IllegalArgumentException if any dimension is not a power of two 546 */ 547 public Object mdfft(Object mdca, boolean forward) 548 throws IllegalArgumentException { 549 MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix) 550 new MultiDimensionalComplexMatrix(mdca).clone(); 551 int[] dimensionSize = mdcm.getDimensionSizes(); 552 //cycle through each dimension 553 for (int i = 0; i < dimensionSize.length; i++) { 554 mdfft(mdcm, forward, i, new int[0]); 555 } 556 return mdcm.getArray(); 557 } 558 559 /** 560 * Performs one dimension of a multi-dimensional Fourier transform. 561 * 562 * @param mdcm input matrix 563 * @param forward inverseTransform2 is preformed if this is false 564 * @param d index of the dimension to process 565 * @param subVector recursion subvector 566 * @throws IllegalArgumentException if any dimension is not a power of two 567 */ 568 private void mdfft(MultiDimensionalComplexMatrix mdcm, boolean forward, 569 int d, int[] subVector) 570 throws IllegalArgumentException { 571 int[] dimensionSize = mdcm.getDimensionSizes(); 572 //if done 573 if (subVector.length == dimensionSize.length) { 574 Complex[] temp = new Complex[dimensionSize[d]]; 575 for (int i = 0; i < dimensionSize[d]; i++) { 576 //fft along dimension d 577 subVector[d] = i; 578 temp[i] = mdcm.get(subVector); 579 } 580 581 if (forward) 582 temp = transform2(temp); 583 else 584 temp = inversetransform2(temp); 585 586 for (int i = 0; i < dimensionSize[d]; i++) { 587 subVector[d] = i; 588 mdcm.set(temp[i], subVector); 589 } 590 } else { 591 int[] vector = new int[subVector.length + 1]; 592 System.arraycopy(subVector, 0, vector, 0, subVector.length); 593 if (subVector.length == d) { 594 //value is not important once the recursion is done. 595 //then an fft will be applied along the dimension d. 596 vector[d] = 0; 597 mdfft(mdcm, forward, d, vector); 598 } else { 599 for (int i = 0; i < dimensionSize[subVector.length]; i++) { 600 vector[subVector.length] = i; 601 //further split along the next dimension 602 mdfft(mdcm, forward, d, vector); 603 } 604 } 605 } 606 return; 607 } 608 609 /** 610 * Complex matrix implementation. 611 * Not designed for synchronized access 612 * may eventually be replaced by jsr-83 of the java community process 613 * http://jcp.org/en/jsr/detail?id=83 614 * may require additional exception throws for other basic requirements. 615 */ 616 private static class MultiDimensionalComplexMatrix 617 implements Cloneable { 618 619 /** Size in all dimensions. */ 620 protected int[] dimensionSize; 621 622 /** Storage array. */ 623 protected Object multiDimensionalComplexArray; 624 625 /** Simple constructor. 626 * @param multiDimensionalComplexArray array containing the matrix elements 627 */ 628 public MultiDimensionalComplexMatrix(Object multiDimensionalComplexArray) { 629 630 this.multiDimensionalComplexArray = multiDimensionalComplexArray; 631 632 // count dimensions 633 int numOfDimensions = 0; 634 for (Object lastDimension = multiDimensionalComplexArray; 635 lastDimension instanceof Object[];) { 636 final Object[] array = (Object[]) lastDimension; 637 numOfDimensions++; 638 lastDimension = array[0]; 639 } 640 641 // allocate array with exact count 642 dimensionSize = new int[numOfDimensions]; 643 644 // fill array 645 numOfDimensions = 0; 646 for (Object lastDimension = multiDimensionalComplexArray; 647 lastDimension instanceof Object[];) { 648 final Object[] array = (Object[]) lastDimension; 649 dimensionSize[numOfDimensions++] = array.length; 650 lastDimension = array[0]; 651 } 652 653 } 654 655 /** 656 * Get a matrix element. 657 * @param vector indices of the element 658 * @return matrix element 659 * @exception IllegalArgumentException if dimensions do not match 660 */ 661 public Complex get(int... vector) 662 throws IllegalArgumentException { 663 if (vector == null) { 664 if (dimensionSize.length > 0) { 665 throw MathRuntimeException.createIllegalArgumentException( 666 "some dimensions don't match: {0} != {1}", 667 0, dimensionSize.length); 668 } 669 return null; 670 } 671 if (vector.length != dimensionSize.length) { 672 throw MathRuntimeException.createIllegalArgumentException( 673 "some dimensions don't match: {0} != {1}", 674 vector.length, dimensionSize.length); 675 } 676 677 Object lastDimension = multiDimensionalComplexArray; 678 679 for (int i = 0; i < dimensionSize.length; i++) { 680 lastDimension = ((Object[]) lastDimension)[vector[i]]; 681 } 682 return (Complex) lastDimension; 683 } 684 685 /** 686 * Set a matrix element. 687 * @param magnitude magnitude of the element 688 * @param vector indices of the element 689 * @return the previous value 690 * @exception IllegalArgumentException if dimensions do not match 691 */ 692 public Complex set(Complex magnitude, int... vector) 693 throws IllegalArgumentException { 694 if (vector == null) { 695 if (dimensionSize.length > 0) { 696 throw MathRuntimeException.createIllegalArgumentException( 697 "some dimensions don't match: {0} != {1}", 698 0, dimensionSize.length); 699 } 700 return null; 701 } 702 if (vector.length != dimensionSize.length) { 703 throw MathRuntimeException.createIllegalArgumentException( 704 "some dimensions don't match: {0} != {1}", 705 vector.length,dimensionSize.length); 706 } 707 708 Object[] lastDimension = (Object[]) multiDimensionalComplexArray; 709 for (int i = 0; i < dimensionSize.length - 1; i++) { 710 lastDimension = (Object[]) lastDimension[vector[i]]; 711 } 712 713 Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]]; 714 lastDimension[vector[dimensionSize.length - 1]] = magnitude; 715 716 return lastValue; 717 } 718 719 /** 720 * Get the size in all dimensions. 721 * @return size in all dimensions 722 */ 723 public int[] getDimensionSizes() { 724 return dimensionSize.clone(); 725 } 726 727 /** 728 * Get the underlying storage array 729 * @return underlying storage array 730 */ 731 public Object getArray() { 732 return multiDimensionalComplexArray; 733 } 734 735 /** {@inheritDoc} */ 736 @Override 737 public Object clone() { 738 MultiDimensionalComplexMatrix mdcm = 739 new MultiDimensionalComplexMatrix(Array.newInstance( 740 Complex.class, dimensionSize)); 741 clone(mdcm); 742 return mdcm; 743 } 744 745 /** 746 * Copy contents of current array into mdcm. 747 * @param mdcm array where to copy data 748 */ 749 private void clone(MultiDimensionalComplexMatrix mdcm) { 750 int[] vector = new int[dimensionSize.length]; 751 int size = 1; 752 for (int i = 0; i < dimensionSize.length; i++) { 753 size *= dimensionSize[i]; 754 } 755 int[][] vectorList = new int[size][dimensionSize.length]; 756 for (int[] nextVector: vectorList) { 757 System.arraycopy(vector, 0, nextVector, 0, 758 dimensionSize.length); 759 for (int i = 0; i < dimensionSize.length; i++) { 760 vector[i]++; 761 if (vector[i] < dimensionSize[i]) { 762 break; 763 } else { 764 vector[i] = 0; 765 } 766 } 767 } 768 769 for (int[] nextVector: vectorList) { 770 mdcm.set(get(nextVector), nextVector); 771 } 772 } 773 } 774 775 776 /** Computes the n<sup>th</sup> roots of unity. 777 * A cache of already computed values is maintained. 778 */ 779 private static class RootsOfUnity implements Serializable { 780 781 /** Serializable version id. */ 782 private static final long serialVersionUID = 6404784357747329667L; 783 784 /** Number of roots of unity. */ 785 private int omegaCount; 786 787 /** Real part of the roots. */ 788 private double[] omegaReal; 789 790 /** Imaginary part of the roots for forward transform. */ 791 private double[] omegaImaginaryForward; 792 793 /** Imaginary part of the roots for reverse transform. */ 794 private double[] omegaImaginaryInverse; 795 796 /** Forward/reverse indicator. */ 797 private boolean isForward; 798 799 /** 800 * Build an engine for computing then <sup>th</sup> roots of unity 801 */ 802 public RootsOfUnity() { 803 804 omegaCount = 0; 805 omegaReal = null; 806 omegaImaginaryForward = null; 807 omegaImaginaryInverse = null; 808 isForward = true; 809 810 } 811 812 /** 813 * Check if computation has been done for forward or reverse transform. 814 * @return true if computation has been done for forward transform 815 * @throws IllegalStateException if no roots of unity have been computed yet 816 */ 817 public synchronized boolean isForward() throws IllegalStateException { 818 819 if (omegaCount == 0) { 820 throw MathRuntimeException.createIllegalStateException( 821 "roots of unity have not been computed yet"); 822 } 823 return isForward; 824 825 } 826 827 /** Computes the n<sup>th</sup> roots of unity. 828 * <p>The computed omega[] = { 1, w, w<sup>2</sup>, ... w<sup>(n-1)</sup> } where 829 * w = exp(-2 π i / n), i = &sqrt;(-1).</p> 830 * <p>Note that n is positive for 831 * forward transform and negative for inverse transform.</p> 832 * @param n number of roots of unity to compute, 833 * positive for forward transform, negative for inverse transform 834 * @throws IllegalArgumentException if n = 0 835 */ 836 public synchronized void computeOmega(int n) throws IllegalArgumentException { 837 838 if (n == 0) { 839 throw MathRuntimeException.createIllegalArgumentException( 840 "cannot compute 0-th root of unity, indefinite result"); 841 } 842 843 isForward = (n > 0); 844 845 // avoid repetitive calculations 846 final int absN = Math.abs(n); 847 848 if (absN == omegaCount) { 849 return; 850 } 851 852 // calculate everything from scratch, for both forward and inverse versions 853 final double t = 2.0 * Math.PI / absN; 854 final double cosT = Math.cos(t); 855 final double sinT = Math.sin(t); 856 omegaReal = new double[absN]; 857 omegaImaginaryForward = new double[absN]; 858 omegaImaginaryInverse = new double[absN]; 859 omegaReal[0] = 1.0; 860 omegaImaginaryForward[0] = 0.0; 861 omegaImaginaryInverse[0] = 0.0; 862 for (int i = 1; i < absN; i++) { 863 omegaReal[i] = 864 omegaReal[i-1] * cosT + omegaImaginaryForward[i-1] * sinT; 865 omegaImaginaryForward[i] = 866 omegaImaginaryForward[i-1] * cosT - omegaReal[i-1] * sinT; 867 omegaImaginaryInverse[i] = -omegaImaginaryForward[i]; 868 } 869 omegaCount = absN; 870 871 } 872 873 /** 874 * Get the real part of the k<sup>th</sup> n<sup>th</sup> root of unity 875 * @param k index of the n<sup>th</sup> root of unity 876 * @return real part of the k<sup>th</sup> n<sup>th</sup> root of unity 877 * @throws IllegalStateException if no roots of unity have been computed yet 878 * @throws IllegalArgumentException if k is out of range 879 */ 880 public synchronized double getOmegaReal(int k) 881 throws IllegalStateException, IllegalArgumentException { 882 883 if (omegaCount == 0) { 884 throw MathRuntimeException.createIllegalStateException( 885 "roots of unity have not been computed yet"); 886 } 887 if ((k < 0) || (k >= omegaCount)) { 888 throw MathRuntimeException.createIllegalArgumentException( 889 "out of range root of unity index {0} (must be in [{1};{2}])", 890 k, 0, omegaCount - 1); 891 } 892 893 return omegaReal[k]; 894 895 } 896 897 /** 898 * Get the imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity 899 * @param k index of the n<sup>th</sup> root of unity 900 * @return imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity 901 * @throws IllegalStateException if no roots of unity have been computed yet 902 * @throws IllegalArgumentException if k is out of range 903 */ 904 public synchronized double getOmegaImaginary(int k) 905 throws IllegalStateException, IllegalArgumentException { 906 907 if (omegaCount == 0) { 908 throw MathRuntimeException.createIllegalStateException( 909 "roots of unity have not been computed yet"); 910 } 911 if ((k < 0) || (k >= omegaCount)) { 912 throw MathRuntimeException.createIllegalArgumentException( 913 "out of range root of unity index {0} (must be in [{1};{2}])", 914 k, 0, omegaCount - 1); 915 } 916 917 return (isForward) ? 918 omegaImaginaryForward[k] : omegaImaginaryInverse[k]; 919 920 } 921 922 } 923 924 }