Java tutorial
/* * 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 fvr. 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]; } } }