/*
* 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;
/**
* 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: 790243 $ $Date: 2009-07-01 12:03:28 -0400 (Wed, 01 Jul 2009) $
* @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 / Math.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 / Math.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 / Math.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 / Math.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 / Math.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 / Math.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 {
int i, j, k, m, N = data.length;
Complex A, B, C, D, E, F, z, 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
j = 0;
for (i = 0; i < N; i++) {
f[i] = data[j];
k = N >> 1;
while (j >= k && k > 0) {
j -= k; k >>= 1;
}
j += k;
}
// the bottom base-4 round
for (i = 0; i < N; i += 4) {
A = f[i].add(f[i+1]);
B = f[i+2].add(f[i+3]);
C = f[i].subtract(f[i+1]);
D = f[i+2].subtract(f[i+3]);
E = C.add(D.multiply(Complex.I));
F = 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() ? F : E;
f[i+3] = roots.isForward() ? E : F;
}
// iterations from bottom to top take O(N*logN) time
for (i = 4; i < N; i <<= 1) {
m = N / (i<<1);
for (j = 0; j < N; j += i<<1) {
for (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]);
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(
"number of sample is not positive: {0}",
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(
"{0} is not a power of 2, consider padding for fix",
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(
"{0} is not a power of 2, consider padding for fix",
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(
"endpoints do not specify an interval: [{0}, {1}]",
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(
"some dimensions don't match: {0} != {1}",
0, dimensionSize.length);
}
return null;
}
if (vector.length != dimensionSize.length) {
throw MathRuntimeException.createIllegalArgumentException(
"some dimensions don't match: {0} != {1}",
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(
"some dimensions don't match: {0} != {1}",
0, dimensionSize.length);
}
return null;
}
if (vector.length != dimensionSize.length) {
throw MathRuntimeException.createIllegalArgumentException(
"some dimensions don't match: {0} != {1}",
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(
"roots of unity have not been 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(
"cannot compute 0-th root of unity, indefinite result");
}
isForward = (n > 0);
// avoid repetitive calculations
final int absN = Math.abs(n);
if (absN == omegaCount) {
return;
}
// calculate everything from scratch, for both forward and inverse versions
final double t = 2.0 * Math.PI / absN;
final double cosT = Math.cos(t);
final double sinT = Math.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(
"roots of unity have not been computed yet");
}
if ((k < 0) || (k >= omegaCount)) {
throw MathRuntimeException.createIllegalArgumentException(
"out of range root of unity index {0} (must be in [{1};{2}])",
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(
"roots of unity have not been computed yet");
}
if ((k < 0) || (k >= omegaCount)) {
throw MathRuntimeException.createIllegalArgumentException(
"out of range root of unity index {0} (must be in [{1};{2}])",
k, 0, omegaCount - 1);
}
return (isForward) ?
omegaImaginaryForward[k] : omegaImaginaryInverse[k];
}
}
}