FastCosineTransform.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math4.transform;

  19. import java.util.function.UnaryOperator;
  20. import java.util.function.DoubleUnaryOperator;

  22. import org.apache.commons.numbers.complex.Complex;
  23. import org.apache.commons.numbers.core.ArithmeticUtils;

  25. /**
  26.  * Implements the Fast Cosine Transform for transformation of one-dimensional
  27.  * real data sets. For reference, see James S. Walker, <em>Fast Fourier
  28.  * Transforms</em>, chapter 3 (ISBN 0849371635).
  29.  * <p>
  30.  * There are several variants of the discrete cosine transform. The present
  31.  * implementation corresponds to DCT-I, with various normalization conventions,
  32.  * which are specified by the parameter {@link Norm}.
  33.  * <p>
  34.  * DCT-I is equivalent to DFT of an <em>even extension</em> of the data series.
  35.  * More precisely, if x<sub>0</sub>, &hellip;, x<sub>N-1</sub> is the data set
  36.  * to be cosine transformed, the extended data set
  37.  * x<sub>0</sub><sup>&#35;</sup>, &hellip;, x<sub>2N-3</sub><sup>&#35;</sup>
  38.  * is defined as follows
  39.  * <ul>
  40.  * <li>x<sub>k</sub><sup>&#35;</sup> = x<sub>k</sub> if 0 &le; k &lt; N,</li>
  41.  * <li>x<sub>k</sub><sup>&#35;</sup> = x<sub>2N-2-k</sub>
  42.  * if N &le; k &lt; 2N - 2.</li>
  43.  * </ul>
  44.  * <p>
  45.  * Then, the standard DCT-I y<sub>0</sub>, &hellip;, y<sub>N-1</sub> of the real
  46.  * data set x<sub>0</sub>, &hellip;, x<sub>N-1</sub> is equal to <em>half</em>
  47.  * of the N first elements of the DFT of the extended data set
  48.  * x<sub>0</sub><sup>&#35;</sup>, &hellip;, x<sub>2N-3</sub><sup>&#35;</sup>
  49.  * <br>
  50.  * y<sub>n</sub> = (1 / 2) &sum;<sub>k=0</sub><sup>2N-3</sup>
  51.  * x<sub>k</sub><sup>&#35;</sup> exp[-2&pi;i nk / (2N - 2)]
  52.  * &nbsp;&nbsp;&nbsp;&nbsp;k = 0, &hellip;, N-1.
  53.  * <p>
  54.  * The present implementation of the discrete cosine transform as a fast cosine
  55.  * transform requires the length of the data set to be a power of two plus one
  56.  * (N&nbsp;=&nbsp;2<sup>n</sup>&nbsp;+&nbsp;1). Besides, it implicitly assumes
  57.  * that the sampled function is even.
  58.  */
  59. public class FastCosineTransform implements RealTransform {
  60.     /** Operation to be performed. */
  61.     private final UnaryOperator<double[]> op;

  63.     /**
  64.      * @param normalization Normalization to be applied to the
  65.      * transformed data.
  66.      * @param inverse Whether to perform the inverse transform.
  67.      */
  68.     public FastCosineTransform(final Norm normalization,
  69.                                final boolean inverse) {
  70.         op = create(normalization, inverse);
  71.     }

  73.     /**
  74.      * @param normalization Normalization to be applied to the
  75.      * transformed data.
  76.      */
  77.     public FastCosineTransform(final Norm normalization) {
  78.         this(normalization, false);
  79.     }

  81.     /**
  82.      * {@inheritDoc}
  83.      *
  84.      * @throws IllegalArgumentException if the length of the data array is
  85.      * not a power of two plus one.
  86.      */
  87.     @Override
  88.     public double[] apply(final double[] f) {
  89.         return op.apply(f);
  90.     }

  92.     /**
  93.      * {@inheritDoc}
  94.      *
  95.      * @throws IllegalArgumentException if the number of sample points is
  96.      * not a power of two plus one, if the lower bound is greater than or
  97.      * equal to the upper bound, if the number of sample points is negative.
  98.      */
  99.     @Override
  100.     public double[] apply(final DoubleUnaryOperator f,
  101.                           final double min,
  102.                           final double max,
  103.                           final int n) {
  104.         return apply(TransformUtils.sample(f, min, max, n));
  105.     }

  107.     /**
  108.      * Perform the FCT algorithm (including inverse).
  109.      *
  110.      * @param f Data to be transformed.
  111.      * @return the transformed array.
  112.      * @throws IllegalArgumentException if the length of the data array is
  113.      * not a power of two plus one.
  114.      */
  115.     private double[] fct(double[] f) {
  116.         final int n = f.length - 1;
  117.         if (!ArithmeticUtils.isPowerOfTwo(n)) {
  118.             throw new TransformException(TransformException.NOT_POWER_OF_TWO_PLUS_ONE,
  119.                                          Integer.valueOf(f.length));
  120.         }

  122.         final double[] transformed = new double[f.length];

  124.         if (n == 1) {       // trivial case
  125.             transformed[0] = 0.5 * (f[0] + f[1]);
  126.             transformed[1] = 0.5 * (f[0] - f[1]);
  127.             return transformed;
  128.         }

  130.         // construct a new array and perform FFT on it
  131.         final double[] x = new double[n];
  132.         x[0] = 0.5 * (f[0] + f[n]);
  133.         final int nShifted = n >> 1;
  134.         x[nShifted] = f[nShifted];
  135.         // temporary variable for transformed[1]
  136.         double t1 = 0.5 * (f[0] - f[n]);
  137.         final double piOverN = Math.PI / n;
  138.         for (int i = 1; i < nShifted; i++) {
  139.             final int nMi = n - i;
  140.             final double fi = f[i];
  141.             final double fnMi = f[nMi];
  142.             final double a = 0.5 * (fi + fnMi);
  143.             final double arg = i * piOverN;
  144.             final double b = Math.sin(arg) * (fi - fnMi);
  145.             final double c = Math.cos(arg) * (fi - fnMi);
  146.             x[i] = a - b;
  147.             x[nMi] = a + b;
  148.             t1 += c;
  149.         }
  150.         final FastFourierTransform transformer = new FastFourierTransform(FastFourierTransform.Norm.STD,
  151.                                                                           false);
  152.         final Complex[] y = transformer.apply(x);

  154.         // reconstruct the FCT result for the original array
  155.         transformed[0] = y[0].getReal();
  156.         transformed[1] = t1;
  157.         for (int i = 1; i < nShifted; i++) {
  158.             final int i2 = 2 * i;
  159.             transformed[i2] = y[i].getReal();
  160.             transformed[i2 + 1] = transformed[i2 - 1] - y[i].getImaginary();
  161.         }
  162.         transformed[n] = y[nShifted].getReal();

  164.         return transformed;
  165.     }

  167.     /**
  168.      * Factory method.
  169.      *
  170.      * @param normalization Normalization to be applied to the
  171.      * transformed data.
  172.      * @param inverse Whether to perform the inverse transform.
  173.      * @return the transform operator.
  174.      */
  175.     private UnaryOperator<double[]> create(final Norm normalization,
  176.                                            final boolean inverse) {
  177.         if (inverse) {
  178.             return normalization == Norm.ORTHO ?
  179.                 f -> TransformUtils.scaleInPlace(fct(f), Math.sqrt(2d / (f.length - 1))) :
  180.                 f -> TransformUtils.scaleInPlace(fct(f), 2d / (f.length - 1));
  181.         } else {
  182.             return normalization == Norm.ORTHO ?
  183.                 f -> TransformUtils.scaleInPlace(fct(f), Math.sqrt(2d / (f.length - 1))) :
  184.                 f -> fct(f);
  185.         }
  186.     }

  188.     /**
  189.      * Normalization types.
  190.      */
  191.     public enum Norm {
  192.         /**
  193.          * Should be passed to the constructor of {@link FastCosineTransform}
  194.          * to use the <em>standard</em> normalization convention.  The standard
  195.          * DCT-I normalization convention is defined as follows
  196.          * <ul>
  197.          * <li>forward transform:
  198.          * y<sub>n</sub> = (1/2) [x<sub>0</sub> + (-1)<sup>n</sup>x<sub>N-1</sub>]
  199.          * + &sum;<sub>k=1</sub><sup>N-2</sup>
  200.          * x<sub>k</sub> cos[&pi; nk / (N - 1)],</li>
  201.          * <li>inverse transform:
  202.          * x<sub>k</sub> = [1 / (N - 1)] [y<sub>0</sub>
  203.          * + (-1)<sup>k</sup>y<sub>N-1</sub>]
  204.          * + [2 / (N - 1)] &sum;<sub>n=1</sub><sup>N-2</sup>
  205.          * y<sub>n</sub> cos[&pi; nk / (N - 1)],</li>
  206.          * </ul>
  207.          * where N is the size of the data sample.
  208.          */
  209.         STD,

  211.         /**
  212.          * Should be passed to the constructor of {@link FastCosineTransform}
  213.          * to use the <em>orthogonal</em> normalization convention. The orthogonal
  214.          * DCT-I normalization convention is defined as follows
  215.          * <ul>
  216.          * <li>forward transform:
  217.          * y<sub>n</sub> = [2(N - 1)]<sup>-1/2</sup> [x<sub>0</sub>
  218.          * + (-1)<sup>n</sup>x<sub>N-1</sub>]
  219.          * + [2 / (N - 1)]<sup>1/2</sup> &sum;<sub>k=1</sub><sup>N-2</sup>
  220.          * x<sub>k</sub> cos[&pi; nk / (N - 1)],</li>
  221.          * <li>inverse transform:
  222.          * x<sub>k</sub> = [2(N - 1)]<sup>-1/2</sup> [y<sub>0</sub>
  223.          * + (-1)<sup>k</sup>y<sub>N-1</sub>]
  224.          * + [2 / (N - 1)]<sup>1/2</sup> &sum;<sub>n=1</sub><sup>N-2</sup>
  225.          * y<sub>n</sub> cos[&pi; nk / (N - 1)],</li>
  226.          * </ul>
  227.          * which makes the transform orthogonal. N is the size of the data sample.
  228.          */
  229.         ORTHO;
  230.     }
  231. }
Created with JaCoCo 0.8.10.202304240956