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