FastCosineTransform.java
/*
* 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.math4.transform;
import java.util.function.UnaryOperator;
import java.util.function.DoubleUnaryOperator;
import org.apache.commons.numbers.complex.Complex;
import org.apache.commons.numbers.core.ArithmeticUtils;
/**
* Implements the Fast Cosine Transform for transformation of one-dimensional
* real data sets. For reference, see James S. Walker, <em>Fast Fourier
* Transforms</em>, chapter 3 (ISBN 0849371635).
* <p>
* There are several variants of the discrete cosine transform. The present
* implementation corresponds to DCT-I, with various normalization conventions,
* which are specified by the parameter {@link Norm}.
* <p>
* DCT-I is equivalent to DFT of an <em>even extension</em> of the data series.
* More precisely, if x<sub>0</sub>, …, x<sub>N-1</sub> is the data set
* to be cosine transformed, the extended data set
* x<sub>0</sub><sup>#</sup>, …, x<sub>2N-3</sub><sup>#</sup>
* is defined as follows
* <ul>
* <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 0 ≤ k < N,</li>
* <li>x<sub>k</sub><sup>#</sup> = x<sub>2N-2-k</sub>
* if N ≤ k < 2N - 2.</li>
* </ul>
* <p>
* Then, the standard DCT-I y<sub>0</sub>, …, y<sub>N-1</sub> of the real
* data set x<sub>0</sub>, …, x<sub>N-1</sub> is equal to <em>half</em>
* of the N first elements of the DFT of the extended data set
* x<sub>0</sub><sup>#</sup>, …, x<sub>2N-3</sub><sup>#</sup>
* <br>
* y<sub>n</sub> = (1 / 2) ∑<sub>k=0</sub><sup>2N-3</sup>
* x<sub>k</sub><sup>#</sup> exp[-2πi nk / (2N - 2)]
* k = 0, …, N-1.
* <p>
* The present implementation of the discrete cosine transform as a fast cosine
* transform requires the length of the data set to be a power of two plus one
* (N = 2<sup>n</sup> + 1). Besides, it implicitly assumes
* that the sampled function is even.
*/
public class FastCosineTransform implements RealTransform {
/** Operation to be performed. */
private final UnaryOperator<double[]> op;
/**
* @param normalization Normalization to be applied to the
* transformed data.
* @param inverse Whether to perform the inverse transform.
*/
public FastCosineTransform(final Norm normalization,
final boolean inverse) {
op = create(normalization, inverse);
}
/**
* @param normalization Normalization to be applied to the
* transformed data.
*/
public FastCosineTransform(final Norm normalization) {
this(normalization, false);
}
/**
* {@inheritDoc}
*
* @throws IllegalArgumentException if the length of the data array is
* not a power of two plus one.
*/
@Override
public double[] apply(final double[] f) {
return op.apply(f);
}
/**
* {@inheritDoc}
*
* @throws IllegalArgumentException if the number of sample points is
* not a power of two plus one, if the lower bound is greater than or
* equal to the upper bound, if the number of sample points is negative.
*/
@Override
public double[] apply(final DoubleUnaryOperator f,
final double min,
final double max,
final int n) {
return apply(TransformUtils.sample(f, min, max, n));
}
/**
* Perform the FCT algorithm (including inverse).
*
* @param f Data to be transformed.
* @return the transformed array.
* @throws IllegalArgumentException if the length of the data array is
* not a power of two plus one.
*/
private double[] fct(double[] f) {
final int n = f.length - 1;
if (!ArithmeticUtils.isPowerOfTwo(n)) {
throw new TransformException(TransformException.NOT_POWER_OF_TWO_PLUS_ONE,
Integer.valueOf(f.length));
}
final double[] transformed = new double[f.length];
if (n == 1) { // trivial case
transformed[0] = 0.5 * (f[0] + f[1]);
transformed[1] = 0.5 * (f[0] - f[1]);
return transformed;
}
// construct a new array and perform FFT on it
final double[] x = new double[n];
x[0] = 0.5 * (f[0] + f[n]);
final int nShifted = n >> 1;
x[nShifted] = f[nShifted];
// temporary variable for transformed[1]
double t1 = 0.5 * (f[0] - f[n]);
final double piOverN = Math.PI / n;
for (int i = 1; i < nShifted; i++) {
final int nMi = n - i;
final double fi = f[i];
final double fnMi = f[nMi];
final double a = 0.5 * (fi + fnMi);
final double arg = i * piOverN;
final double b = Math.sin(arg) * (fi - fnMi);
final double c = Math.cos(arg) * (fi - fnMi);
x[i] = a - b;
x[nMi] = a + b;
t1 += c;
}
final FastFourierTransform transformer = new FastFourierTransform(FastFourierTransform.Norm.STD,
false);
final Complex[] y = transformer.apply(x);
// reconstruct the FCT result for the original array
transformed[0] = y[0].getReal();
transformed[1] = t1;
for (int i = 1; i < nShifted; i++) {
final int i2 = 2 * i;
transformed[i2] = y[i].getReal();
transformed[i2 + 1] = transformed[i2 - 1] - y[i].getImaginary();
}
transformed[n] = y[nShifted].getReal();
return transformed;
}
/**
* Factory method.
*
* @param normalization Normalization to be applied to the
* transformed data.
* @param inverse Whether to perform the inverse transform.
* @return the transform operator.
*/
private UnaryOperator<double[]> create(final Norm normalization,
final boolean inverse) {
if (inverse) {
return normalization == Norm.ORTHO ?
f -> TransformUtils.scaleInPlace(fct(f), Math.sqrt(2d / (f.length - 1))) :
f -> TransformUtils.scaleInPlace(fct(f), 2d / (f.length - 1));
} else {
return normalization == Norm.ORTHO ?
f -> TransformUtils.scaleInPlace(fct(f), Math.sqrt(2d / (f.length - 1))) :
f -> fct(f);
}
}
/**
* Normalization types.
*/
public enum Norm {
/**
* Should be passed to the constructor of {@link FastCosineTransform}
* to use the <em>standard</em> normalization convention. The standard
* DCT-I normalization convention is defined as follows
* <ul>
* <li>forward transform:
* y<sub>n</sub> = (1/2) [x<sub>0</sub> + (-1)<sup>n</sup>x<sub>N-1</sub>]
* + ∑<sub>k=1</sub><sup>N-2</sup>
* x<sub>k</sub> cos[π nk / (N - 1)],</li>
* <li>inverse transform:
* x<sub>k</sub> = [1 / (N - 1)] [y<sub>0</sub>
* + (-1)<sup>k</sup>y<sub>N-1</sub>]
* + [2 / (N - 1)] ∑<sub>n=1</sub><sup>N-2</sup>
* y<sub>n</sub> cos[π nk / (N - 1)],</li>
* </ul>
* where N is the size of the data sample.
*/
STD,
/**
* Should be passed to the constructor of {@link FastCosineTransform}
* to use the <em>orthogonal</em> normalization convention. The orthogonal
* DCT-I normalization convention is defined as follows
* <ul>
* <li>forward transform:
* y<sub>n</sub> = [2(N - 1)]<sup>-1/2</sup> [x<sub>0</sub>
* + (-1)<sup>n</sup>x<sub>N-1</sub>]
* + [2 / (N - 1)]<sup>1/2</sup> ∑<sub>k=1</sub><sup>N-2</sup>
* x<sub>k</sub> cos[π nk / (N - 1)],</li>
* <li>inverse transform:
* x<sub>k</sub> = [2(N - 1)]<sup>-1/2</sup> [y<sub>0</sub>
* + (-1)<sup>k</sup>y<sub>N-1</sub>]
* + [2 / (N - 1)]<sup>1/2</sup> ∑<sub>n=1</sub><sup>N-2</sup>
* y<sub>n</sub> cos[π nk / (N - 1)],</li>
* </ul>
* which makes the transform orthogonal. N is the size of the data sample.
*/
ORTHO;
}
}