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