| /* |
| * 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 double[] transformed = new double[f.length]; |
| |
| final int n = f.length - 1; |
| if (!ArithmeticUtils.isPowerOfTwo(n)) { |
| throw new TransformException(TransformException.NOT_POWER_OF_TWO_PLUS_ONE, |
| Integer.valueOf(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; |
| } |
| } |
