| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.commons.math.transform; |
| |
| import java.io.Serializable; |
| import java.lang.reflect.Array; |
| |
| import org.apache.commons.math.FunctionEvaluationException; |
| import org.apache.commons.math.MathRuntimeException; |
| import org.apache.commons.math.analysis.UnivariateRealFunction; |
| import org.apache.commons.math.complex.Complex; |
| import org.apache.commons.math.exception.util.LocalizedFormats; |
| import org.apache.commons.math.util.FastMath; |
| |
| /** |
| * Implements the <a href="http://mathworld.wolfram.com/FastFourierTransform.html"> |
| * Fast Fourier Transform</a> for transformation of one-dimensional data sets. |
| * For reference, see <b>Applied Numerical Linear Algebra</b>, ISBN 0898713897, |
| * chapter 6. |
| * <p> |
| * There are several conventions for the definition of FFT and inverse FFT, |
| * mainly on different coefficient and exponent. Here the equations are listed |
| * in the comments of the corresponding methods.</p> |
| * <p> |
| * We require the length of data set to be power of 2, this greatly simplifies |
| * and speeds up the code. Users can pad the data with zeros to meet this |
| * requirement. There are other flavors of FFT, for reference, see S. Winograd, |
| * <i>On computing the discrete Fourier transform</i>, Mathematics of Computation, |
| * 32 (1978), 175 - 199.</p> |
| * |
| * @version $Revision: 1070725 $ $Date: 2011-02-15 02:31:12 +0100 (mar. 15 févr. 2011) $ |
| * @since 1.2 |
| */ |
| public class FastFourierTransformer implements Serializable { |
| |
| /** Serializable version identifier. */ |
| static final long serialVersionUID = 5138259215438106000L; |
| |
| /** roots of unity */ |
| private RootsOfUnity roots = new RootsOfUnity(); |
| |
| /** |
| * Construct a default transformer. |
| */ |
| public FastFourierTransformer() { |
| super(); |
| } |
| |
| /** |
| * Transform the given real data set. |
| * <p> |
| * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $ |
| * </p> |
| * |
| * @param f the real data array to be transformed |
| * @return the complex transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] transform(double f[]) |
| throws IllegalArgumentException { |
| return fft(f, false); |
| } |
| |
| /** |
| * Transform the given real function, sampled on the given interval. |
| * <p> |
| * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $ |
| * </p> |
| * |
| * @param f the function to be sampled and transformed |
| * @param min the lower bound for the interval |
| * @param max the upper bound for the interval |
| * @param n the number of sample points |
| * @return the complex transformed array |
| * @throws FunctionEvaluationException if function cannot be evaluated |
| * at some point |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] transform(UnivariateRealFunction f, |
| double min, double max, int n) |
| throws FunctionEvaluationException, IllegalArgumentException { |
| double data[] = sample(f, min, max, n); |
| return fft(data, false); |
| } |
| |
| /** |
| * Transform the given complex data set. |
| * <p> |
| * The formula is $ y_n = \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k $ |
| * </p> |
| * |
| * @param f the complex data array to be transformed |
| * @return the complex transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] transform(Complex f[]) |
| throws IllegalArgumentException { |
| roots.computeOmega(f.length); |
| return fft(f); |
| } |
| |
| /** |
| * Transform the given real data set. |
| * <p> |
| * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$ |
| * </p> |
| * |
| * @param f the real data array to be transformed |
| * @return the complex transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] transform2(double f[]) |
| throws IllegalArgumentException { |
| |
| double scaling_coefficient = 1.0 / FastMath.sqrt(f.length); |
| return scaleArray(fft(f, false), scaling_coefficient); |
| } |
| |
| /** |
| * Transform the given real function, sampled on the given interval. |
| * <p> |
| * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$ |
| * </p> |
| * |
| * @param f the function to be sampled and transformed |
| * @param min the lower bound for the interval |
| * @param max the upper bound for the interval |
| * @param n the number of sample points |
| * @return the complex transformed array |
| * @throws FunctionEvaluationException if function cannot be evaluated |
| * at some point |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] transform2(UnivariateRealFunction f, |
| double min, double max, int n) |
| throws FunctionEvaluationException, IllegalArgumentException { |
| |
| double data[] = sample(f, min, max, n); |
| double scaling_coefficient = 1.0 / FastMath.sqrt(n); |
| return scaleArray(fft(data, false), scaling_coefficient); |
| } |
| |
| /** |
| * Transform the given complex data set. |
| * <p> |
| * The formula is $y_n = (1/\sqrt{N}) \Sigma_{k=0}^{N-1} e^{-2 \pi i nk/N} x_k$ |
| * </p> |
| * |
| * @param f the complex data array to be transformed |
| * @return the complex transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] transform2(Complex f[]) |
| throws IllegalArgumentException { |
| |
| roots.computeOmega(f.length); |
| double scaling_coefficient = 1.0 / FastMath.sqrt(f.length); |
| return scaleArray(fft(f), scaling_coefficient); |
| } |
| |
| /** |
| * Inversely transform the given real data set. |
| * <p> |
| * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $ |
| * </p> |
| * |
| * @param f the real data array to be inversely transformed |
| * @return the complex inversely transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] inversetransform(double f[]) |
| throws IllegalArgumentException { |
| |
| double scaling_coefficient = 1.0 / f.length; |
| return scaleArray(fft(f, true), scaling_coefficient); |
| } |
| |
| /** |
| * Inversely transform the given real function, sampled on the given interval. |
| * <p> |
| * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $ |
| * </p> |
| * |
| * @param f the function to be sampled and inversely transformed |
| * @param min the lower bound for the interval |
| * @param max the upper bound for the interval |
| * @param n the number of sample points |
| * @return the complex inversely transformed array |
| * @throws FunctionEvaluationException if function cannot be evaluated |
| * at some point |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] inversetransform(UnivariateRealFunction f, |
| double min, double max, int n) |
| throws FunctionEvaluationException, IllegalArgumentException { |
| |
| double data[] = sample(f, min, max, n); |
| double scaling_coefficient = 1.0 / n; |
| return scaleArray(fft(data, true), scaling_coefficient); |
| } |
| |
| /** |
| * Inversely transform the given complex data set. |
| * <p> |
| * The formula is $ x_k = (1/N) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n $ |
| * </p> |
| * |
| * @param f the complex data array to be inversely transformed |
| * @return the complex inversely transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] inversetransform(Complex f[]) |
| throws IllegalArgumentException { |
| |
| roots.computeOmega(-f.length); // pass negative argument |
| double scaling_coefficient = 1.0 / f.length; |
| return scaleArray(fft(f), scaling_coefficient); |
| } |
| |
| /** |
| * Inversely transform the given real data set. |
| * <p> |
| * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$ |
| * </p> |
| * |
| * @param f the real data array to be inversely transformed |
| * @return the complex inversely transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] inversetransform2(double f[]) |
| throws IllegalArgumentException { |
| |
| double scaling_coefficient = 1.0 / FastMath.sqrt(f.length); |
| return scaleArray(fft(f, true), scaling_coefficient); |
| } |
| |
| /** |
| * Inversely transform the given real function, sampled on the given interval. |
| * <p> |
| * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$ |
| * </p> |
| * |
| * @param f the function to be sampled and inversely transformed |
| * @param min the lower bound for the interval |
| * @param max the upper bound for the interval |
| * @param n the number of sample points |
| * @return the complex inversely transformed array |
| * @throws FunctionEvaluationException if function cannot be evaluated |
| * at some point |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] inversetransform2(UnivariateRealFunction f, |
| double min, double max, int n) |
| throws FunctionEvaluationException, IllegalArgumentException { |
| |
| double data[] = sample(f, min, max, n); |
| double scaling_coefficient = 1.0 / FastMath.sqrt(n); |
| return scaleArray(fft(data, true), scaling_coefficient); |
| } |
| |
| /** |
| * Inversely transform the given complex data set. |
| * <p> |
| * The formula is $x_k = (1/\sqrt{N}) \Sigma_{n=0}^{N-1} e^{2 \pi i nk/N} y_n$ |
| * </p> |
| * |
| * @param f the complex data array to be inversely transformed |
| * @return the complex inversely transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public Complex[] inversetransform2(Complex f[]) |
| throws IllegalArgumentException { |
| |
| roots.computeOmega(-f.length); // pass negative argument |
| double scaling_coefficient = 1.0 / FastMath.sqrt(f.length); |
| return scaleArray(fft(f), scaling_coefficient); |
| } |
| |
| /** |
| * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse). |
| * |
| * @param f the real data array to be transformed |
| * @param isInverse the indicator of forward or inverse transform |
| * @return the complex transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| protected Complex[] fft(double f[], boolean isInverse) |
| throws IllegalArgumentException { |
| |
| verifyDataSet(f); |
| Complex F[] = new Complex[f.length]; |
| if (f.length == 1) { |
| F[0] = new Complex(f[0], 0.0); |
| return F; |
| } |
| |
| // Rather than the naive real to complex conversion, pack 2N |
| // real numbers into N complex numbers for better performance. |
| int N = f.length >> 1; |
| Complex c[] = new Complex[N]; |
| for (int i = 0; i < N; i++) { |
| c[i] = new Complex(f[2*i], f[2*i+1]); |
| } |
| roots.computeOmega(isInverse ? -N : N); |
| Complex z[] = fft(c); |
| |
| // reconstruct the FFT result for the original array |
| roots.computeOmega(isInverse ? -2*N : 2*N); |
| F[0] = new Complex(2 * (z[0].getReal() + z[0].getImaginary()), 0.0); |
| F[N] = new Complex(2 * (z[0].getReal() - z[0].getImaginary()), 0.0); |
| for (int i = 1; i < N; i++) { |
| Complex A = z[N-i].conjugate(); |
| Complex B = z[i].add(A); |
| Complex C = z[i].subtract(A); |
| //Complex D = roots.getOmega(i).multiply(Complex.I); |
| Complex D = new Complex(-roots.getOmegaImaginary(i), |
| roots.getOmegaReal(i)); |
| F[i] = B.subtract(C.multiply(D)); |
| F[2*N-i] = F[i].conjugate(); |
| } |
| |
| return scaleArray(F, 0.5); |
| } |
| |
| /** |
| * Perform the base-4 Cooley-Tukey FFT algorithm (including inverse). |
| * |
| * @param data the complex data array to be transformed |
| * @return the complex transformed array |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| protected Complex[] fft(Complex data[]) |
| throws IllegalArgumentException { |
| |
| final int n = data.length; |
| final Complex f[] = new Complex[n]; |
| |
| // initial simple cases |
| verifyDataSet(data); |
| if (n == 1) { |
| f[0] = data[0]; |
| return f; |
| } |
| if (n == 2) { |
| f[0] = data[0].add(data[1]); |
| f[1] = data[0].subtract(data[1]); |
| return f; |
| } |
| |
| // permute original data array in bit-reversal order |
| int ii = 0; |
| for (int i = 0; i < n; i++) { |
| f[i] = data[ii]; |
| int k = n >> 1; |
| while (ii >= k && k > 0) { |
| ii -= k; k >>= 1; |
| } |
| ii += k; |
| } |
| |
| // the bottom base-4 round |
| for (int i = 0; i < n; i += 4) { |
| final Complex a = f[i].add(f[i+1]); |
| final Complex b = f[i+2].add(f[i+3]); |
| final Complex c = f[i].subtract(f[i+1]); |
| final Complex d = f[i+2].subtract(f[i+3]); |
| final Complex e1 = c.add(d.multiply(Complex.I)); |
| final Complex e2 = c.subtract(d.multiply(Complex.I)); |
| f[i] = a.add(b); |
| f[i+2] = a.subtract(b); |
| // omegaCount indicates forward or inverse transform |
| f[i+1] = roots.isForward() ? e2 : e1; |
| f[i+3] = roots.isForward() ? e1 : e2; |
| } |
| |
| // iterations from bottom to top take O(N*logN) time |
| for (int i = 4; i < n; i <<= 1) { |
| final int m = n / (i<<1); |
| for (int j = 0; j < n; j += i<<1) { |
| for (int k = 0; k < i; k++) { |
| //z = f[i+j+k].multiply(roots.getOmega(k*m)); |
| final int k_times_m = k*m; |
| final double omega_k_times_m_real = roots.getOmegaReal(k_times_m); |
| final double omega_k_times_m_imaginary = roots.getOmegaImaginary(k_times_m); |
| //z = f[i+j+k].multiply(omega[k*m]); |
| final Complex z = new Complex( |
| f[i+j+k].getReal() * omega_k_times_m_real - |
| f[i+j+k].getImaginary() * omega_k_times_m_imaginary, |
| f[i+j+k].getReal() * omega_k_times_m_imaginary + |
| f[i+j+k].getImaginary() * omega_k_times_m_real); |
| |
| f[i+j+k] = f[j+k].subtract(z); |
| f[j+k] = f[j+k].add(z); |
| } |
| } |
| } |
| return f; |
| } |
| |
| /** |
| * Sample the given univariate real function on the given interval. |
| * <p> |
| * The interval is divided equally into N sections and sample points |
| * are taken from min to max-(max-min)/N. Usually f(x) is periodic |
| * such that f(min) = f(max) (note max is not sampled), but we don't |
| * require that.</p> |
| * |
| * @param f the function to be sampled |
| * @param min the lower bound for the interval |
| * @param max the upper bound for the interval |
| * @param n the number of sample points |
| * @return the samples array |
| * @throws FunctionEvaluationException if function cannot be evaluated at some point |
| * @throws IllegalArgumentException if any parameters are invalid |
| */ |
| public static double[] sample(UnivariateRealFunction f, double min, double max, int n) |
| throws FunctionEvaluationException, IllegalArgumentException { |
| |
| if (n <= 0) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.NOT_POSITIVE_NUMBER_OF_SAMPLES, |
| n); |
| } |
| verifyInterval(min, max); |
| |
| double s[] = new double[n]; |
| double h = (max - min) / n; |
| for (int i = 0; i < n; i++) { |
| s[i] = f.value(min + i * h); |
| } |
| return s; |
| } |
| |
| /** |
| * Multiply every component in the given real array by the |
| * given real number. The change is made in place. |
| * |
| * @param f the real array to be scaled |
| * @param d the real scaling coefficient |
| * @return a reference to the scaled array |
| */ |
| public static double[] scaleArray(double f[], double d) { |
| for (int i = 0; i < f.length; i++) { |
| f[i] *= d; |
| } |
| return f; |
| } |
| |
| /** |
| * Multiply every component in the given complex array by the |
| * given real number. The change is made in place. |
| * |
| * @param f the complex array to be scaled |
| * @param d the real scaling coefficient |
| * @return a reference to the scaled array |
| */ |
| public static Complex[] scaleArray(Complex f[], double d) { |
| for (int i = 0; i < f.length; i++) { |
| f[i] = new Complex(d * f[i].getReal(), d * f[i].getImaginary()); |
| } |
| return f; |
| } |
| |
| /** |
| * Returns true if the argument is power of 2. |
| * |
| * @param n the number to test |
| * @return true if the argument is power of 2 |
| */ |
| public static boolean isPowerOf2(long n) { |
| return (n > 0) && ((n & (n - 1)) == 0); |
| } |
| |
| /** |
| * Verifies that the data set has length of power of 2. |
| * |
| * @param d the data array |
| * @throws IllegalArgumentException if array length is not power of 2 |
| */ |
| public static void verifyDataSet(double d[]) throws IllegalArgumentException { |
| if (!isPowerOf2(d.length)) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, d.length); |
| } |
| } |
| |
| /** |
| * Verifies that the data set has length of power of 2. |
| * |
| * @param o the data array |
| * @throws IllegalArgumentException if array length is not power of 2 |
| */ |
| public static void verifyDataSet(Object o[]) throws IllegalArgumentException { |
| if (!isPowerOf2(o.length)) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, o.length); |
| } |
| } |
| |
| /** |
| * Verifies that the endpoints specify an interval. |
| * |
| * @param lower lower endpoint |
| * @param upper upper endpoint |
| * @throws IllegalArgumentException if not interval |
| */ |
| public static void verifyInterval(double lower, double upper) |
| throws IllegalArgumentException { |
| |
| if (lower >= upper) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.ENDPOINTS_NOT_AN_INTERVAL, |
| lower, upper); |
| } |
| } |
| |
| /** |
| * Performs a multi-dimensional Fourier transform on a given array. |
| * Use {@link #inversetransform2(Complex[])} and |
| * {@link #transform2(Complex[])} in a row-column implementation |
| * in any number of dimensions with O(N×log(N)) complexity with |
| * N=n<sub>1</sub>×n<sub>2</sub>×n<sub>3</sub>×...×n<sub>d</sub>, |
| * n<sub>x</sub>=number of elements in dimension x, |
| * and d=total number of dimensions. |
| * |
| * @param mdca Multi-Dimensional Complex Array id est Complex[][][][] |
| * @param forward inverseTransform2 is preformed if this is false |
| * @return transform of mdca as a Multi-Dimensional Complex Array id est Complex[][][][] |
| * @throws IllegalArgumentException if any dimension is not a power of two |
| */ |
| public Object mdfft(Object mdca, boolean forward) |
| throws IllegalArgumentException { |
| MultiDimensionalComplexMatrix mdcm = (MultiDimensionalComplexMatrix) |
| new MultiDimensionalComplexMatrix(mdca).clone(); |
| int[] dimensionSize = mdcm.getDimensionSizes(); |
| //cycle through each dimension |
| for (int i = 0; i < dimensionSize.length; i++) { |
| mdfft(mdcm, forward, i, new int[0]); |
| } |
| return mdcm.getArray(); |
| } |
| |
| /** |
| * Performs one dimension of a multi-dimensional Fourier transform. |
| * |
| * @param mdcm input matrix |
| * @param forward inverseTransform2 is preformed if this is false |
| * @param d index of the dimension to process |
| * @param subVector recursion subvector |
| * @throws IllegalArgumentException if any dimension is not a power of two |
| */ |
| private void mdfft(MultiDimensionalComplexMatrix mdcm, boolean forward, |
| int d, int[] subVector) |
| throws IllegalArgumentException { |
| int[] dimensionSize = mdcm.getDimensionSizes(); |
| //if done |
| if (subVector.length == dimensionSize.length) { |
| Complex[] temp = new Complex[dimensionSize[d]]; |
| for (int i = 0; i < dimensionSize[d]; i++) { |
| //fft along dimension d |
| subVector[d] = i; |
| temp[i] = mdcm.get(subVector); |
| } |
| |
| if (forward) |
| temp = transform2(temp); |
| else |
| temp = inversetransform2(temp); |
| |
| for (int i = 0; i < dimensionSize[d]; i++) { |
| subVector[d] = i; |
| mdcm.set(temp[i], subVector); |
| } |
| } else { |
| int[] vector = new int[subVector.length + 1]; |
| System.arraycopy(subVector, 0, vector, 0, subVector.length); |
| if (subVector.length == d) { |
| //value is not important once the recursion is done. |
| //then an fft will be applied along the dimension d. |
| vector[d] = 0; |
| mdfft(mdcm, forward, d, vector); |
| } else { |
| for (int i = 0; i < dimensionSize[subVector.length]; i++) { |
| vector[subVector.length] = i; |
| //further split along the next dimension |
| mdfft(mdcm, forward, d, vector); |
| } |
| } |
| } |
| return; |
| } |
| |
| /** |
| * Complex matrix implementation. |
| * Not designed for synchronized access |
| * may eventually be replaced by jsr-83 of the java community process |
| * http://jcp.org/en/jsr/detail?id=83 |
| * may require additional exception throws for other basic requirements. |
| */ |
| private static class MultiDimensionalComplexMatrix |
| implements Cloneable { |
| |
| /** Size in all dimensions. */ |
| protected int[] dimensionSize; |
| |
| /** Storage array. */ |
| protected Object multiDimensionalComplexArray; |
| |
| /** Simple constructor. |
| * @param multiDimensionalComplexArray array containing the matrix elements |
| */ |
| public MultiDimensionalComplexMatrix(Object multiDimensionalComplexArray) { |
| |
| this.multiDimensionalComplexArray = multiDimensionalComplexArray; |
| |
| // count dimensions |
| int numOfDimensions = 0; |
| for (Object lastDimension = multiDimensionalComplexArray; |
| lastDimension instanceof Object[];) { |
| final Object[] array = (Object[]) lastDimension; |
| numOfDimensions++; |
| lastDimension = array[0]; |
| } |
| |
| // allocate array with exact count |
| dimensionSize = new int[numOfDimensions]; |
| |
| // fill array |
| numOfDimensions = 0; |
| for (Object lastDimension = multiDimensionalComplexArray; |
| lastDimension instanceof Object[];) { |
| final Object[] array = (Object[]) lastDimension; |
| dimensionSize[numOfDimensions++] = array.length; |
| lastDimension = array[0]; |
| } |
| |
| } |
| |
| /** |
| * Get a matrix element. |
| * @param vector indices of the element |
| * @return matrix element |
| * @exception IllegalArgumentException if dimensions do not match |
| */ |
| public Complex get(int... vector) |
| throws IllegalArgumentException { |
| if (vector == null) { |
| if (dimensionSize.length > 0) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, 0, dimensionSize.length); |
| } |
| return null; |
| } |
| if (vector.length != dimensionSize.length) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, vector.length, dimensionSize.length); |
| } |
| |
| Object lastDimension = multiDimensionalComplexArray; |
| |
| for (int i = 0; i < dimensionSize.length; i++) { |
| lastDimension = ((Object[]) lastDimension)[vector[i]]; |
| } |
| return (Complex) lastDimension; |
| } |
| |
| /** |
| * Set a matrix element. |
| * @param magnitude magnitude of the element |
| * @param vector indices of the element |
| * @return the previous value |
| * @exception IllegalArgumentException if dimensions do not match |
| */ |
| public Complex set(Complex magnitude, int... vector) |
| throws IllegalArgumentException { |
| if (vector == null) { |
| if (dimensionSize.length > 0) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, 0, dimensionSize.length); |
| } |
| return null; |
| } |
| if (vector.length != dimensionSize.length) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.DIMENSIONS_MISMATCH_SIMPLE, vector.length,dimensionSize.length); |
| } |
| |
| Object[] lastDimension = (Object[]) multiDimensionalComplexArray; |
| for (int i = 0; i < dimensionSize.length - 1; i++) { |
| lastDimension = (Object[]) lastDimension[vector[i]]; |
| } |
| |
| Complex lastValue = (Complex) lastDimension[vector[dimensionSize.length - 1]]; |
| lastDimension[vector[dimensionSize.length - 1]] = magnitude; |
| |
| return lastValue; |
| } |
| |
| /** |
| * Get the size in all dimensions. |
| * @return size in all dimensions |
| */ |
| public int[] getDimensionSizes() { |
| return dimensionSize.clone(); |
| } |
| |
| /** |
| * Get the underlying storage array |
| * @return underlying storage array |
| */ |
| public Object getArray() { |
| return multiDimensionalComplexArray; |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public Object clone() { |
| MultiDimensionalComplexMatrix mdcm = |
| new MultiDimensionalComplexMatrix(Array.newInstance( |
| Complex.class, dimensionSize)); |
| clone(mdcm); |
| return mdcm; |
| } |
| |
| /** |
| * Copy contents of current array into mdcm. |
| * @param mdcm array where to copy data |
| */ |
| private void clone(MultiDimensionalComplexMatrix mdcm) { |
| int[] vector = new int[dimensionSize.length]; |
| int size = 1; |
| for (int i = 0; i < dimensionSize.length; i++) { |
| size *= dimensionSize[i]; |
| } |
| int[][] vectorList = new int[size][dimensionSize.length]; |
| for (int[] nextVector: vectorList) { |
| System.arraycopy(vector, 0, nextVector, 0, |
| dimensionSize.length); |
| for (int i = 0; i < dimensionSize.length; i++) { |
| vector[i]++; |
| if (vector[i] < dimensionSize[i]) { |
| break; |
| } else { |
| vector[i] = 0; |
| } |
| } |
| } |
| |
| for (int[] nextVector: vectorList) { |
| mdcm.set(get(nextVector), nextVector); |
| } |
| } |
| } |
| |
| |
| /** Computes the n<sup>th</sup> roots of unity. |
| * A cache of already computed values is maintained. |
| */ |
| private static class RootsOfUnity implements Serializable { |
| |
| /** Serializable version id. */ |
| private static final long serialVersionUID = 6404784357747329667L; |
| |
| /** Number of roots of unity. */ |
| private int omegaCount; |
| |
| /** Real part of the roots. */ |
| private double[] omegaReal; |
| |
| /** Imaginary part of the roots for forward transform. */ |
| private double[] omegaImaginaryForward; |
| |
| /** Imaginary part of the roots for reverse transform. */ |
| private double[] omegaImaginaryInverse; |
| |
| /** Forward/reverse indicator. */ |
| private boolean isForward; |
| |
| /** |
| * Build an engine for computing then <sup>th</sup> roots of unity |
| */ |
| public RootsOfUnity() { |
| |
| omegaCount = 0; |
| omegaReal = null; |
| omegaImaginaryForward = null; |
| omegaImaginaryInverse = null; |
| isForward = true; |
| |
| } |
| |
| /** |
| * Check if computation has been done for forward or reverse transform. |
| * @return true if computation has been done for forward transform |
| * @throws IllegalStateException if no roots of unity have been computed yet |
| */ |
| public synchronized boolean isForward() throws IllegalStateException { |
| |
| if (omegaCount == 0) { |
| throw MathRuntimeException.createIllegalStateException(LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET); |
| } |
| return isForward; |
| |
| } |
| |
| /** Computes the n<sup>th</sup> roots of unity. |
| * <p>The computed omega[] = { 1, w, w<sup>2</sup>, ... w<sup>(n-1)</sup> } where |
| * w = exp(-2 π i / n), i = &sqrt;(-1).</p> |
| * <p>Note that n is positive for |
| * forward transform and negative for inverse transform.</p> |
| * @param n number of roots of unity to compute, |
| * positive for forward transform, negative for inverse transform |
| * @throws IllegalArgumentException if n = 0 |
| */ |
| public synchronized void computeOmega(int n) throws IllegalArgumentException { |
| |
| if (n == 0) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.CANNOT_COMPUTE_0TH_ROOT_OF_UNITY); |
| } |
| |
| isForward = n > 0; |
| |
| // avoid repetitive calculations |
| final int absN = FastMath.abs(n); |
| |
| if (absN == omegaCount) { |
| return; |
| } |
| |
| // calculate everything from scratch, for both forward and inverse versions |
| final double t = 2.0 * FastMath.PI / absN; |
| final double cosT = FastMath.cos(t); |
| final double sinT = FastMath.sin(t); |
| omegaReal = new double[absN]; |
| omegaImaginaryForward = new double[absN]; |
| omegaImaginaryInverse = new double[absN]; |
| omegaReal[0] = 1.0; |
| omegaImaginaryForward[0] = 0.0; |
| omegaImaginaryInverse[0] = 0.0; |
| for (int i = 1; i < absN; i++) { |
| omegaReal[i] = |
| omegaReal[i-1] * cosT + omegaImaginaryForward[i-1] * sinT; |
| omegaImaginaryForward[i] = |
| omegaImaginaryForward[i-1] * cosT - omegaReal[i-1] * sinT; |
| omegaImaginaryInverse[i] = -omegaImaginaryForward[i]; |
| } |
| omegaCount = absN; |
| |
| } |
| |
| /** |
| * Get the real part of the k<sup>th</sup> n<sup>th</sup> root of unity |
| * @param k index of the n<sup>th</sup> root of unity |
| * @return real part of the k<sup>th</sup> n<sup>th</sup> root of unity |
| * @throws IllegalStateException if no roots of unity have been computed yet |
| * @throws IllegalArgumentException if k is out of range |
| */ |
| public synchronized double getOmegaReal(int k) |
| throws IllegalStateException, IllegalArgumentException { |
| |
| if (omegaCount == 0) { |
| throw MathRuntimeException.createIllegalStateException(LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET); |
| } |
| if ((k < 0) || (k >= omegaCount)) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX, k, 0, omegaCount - 1); |
| } |
| |
| return omegaReal[k]; |
| |
| } |
| |
| /** |
| * Get the imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity |
| * @param k index of the n<sup>th</sup> root of unity |
| * @return imaginary part of the k<sup>th</sup> n<sup>th</sup> root of unity |
| * @throws IllegalStateException if no roots of unity have been computed yet |
| * @throws IllegalArgumentException if k is out of range |
| */ |
| public synchronized double getOmegaImaginary(int k) |
| throws IllegalStateException, IllegalArgumentException { |
| |
| if (omegaCount == 0) { |
| throw MathRuntimeException.createIllegalStateException(LocalizedFormats.ROOTS_OF_UNITY_NOT_COMPUTED_YET); |
| } |
| if ((k < 0) || (k >= omegaCount)) { |
| throw MathRuntimeException.createIllegalArgumentException( |
| LocalizedFormats.OUT_OF_RANGE_ROOT_OF_UNITY_INDEX, k, 0, omegaCount - 1); |
| } |
| |
| return isForward ? omegaImaginaryForward[k] : omegaImaginaryInverse[k]; |
| |
| } |
| |
| } |
| |
| } |
