FastSineTransform.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;

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

  20. import org.apache.commons.numbers.complex.Complex;
  21. import org.apache.commons.numbers.core.ArithmeticUtils;

  22. /**
  23.  * Implements the Fast Sine Transform for transformation of one-dimensional real
  24.  * data sets. For reference, see James S. Walker, <em>Fast Fourier
  25.  * Transforms</em>, chapter 3 (ISBN 0849371635).
  26.  * <p>
  27.  * There are several variants of the discrete sine transform. The present
  28.  * implementation corresponds to DST-I, with various normalization conventions,
  29.  * which are specified by the parameter {@link Norm}.
  30.  * <strong>It should be noted that regardless to the convention, the first
  31.  * element of the dataset to be transformed must be zero.</strong>
  32.  * <p>
  33.  * DST-I is equivalent to DFT of an <em>odd extension</em> of the data series.
  34.  * More precisely, if x<sub>0</sub>, &hellip;, x<sub>N-1</sub> is the data set
  35.  * to be sine transformed, the extended data set x<sub>0</sub><sup>&#35;</sup>,
  36.  * &hellip;, x<sub>2N-1</sub><sup>&#35;</sup> is defined as follows
  37.  * <ul>
  38.  * <li>x<sub>0</sub><sup>&#35;</sup> = x<sub>0</sub> = 0,</li>
  39.  * <li>x<sub>k</sub><sup>&#35;</sup> = x<sub>k</sub> if 1 &le; k &lt; N,</li>
  40.  * <li>x<sub>N</sub><sup>&#35;</sup> = 0,</li>
  41.  * <li>x<sub>k</sub><sup>&#35;</sup> = -x<sub>2N-k</sub> if N + 1 &le; k &lt;
  42.  * 2N.</li>
  43.  * </ul>
  44.  * <p>
  45.  * Then, the standard DST-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 i (the pure imaginary number) times the N first elements of the DFT of the
  48.  * extended data set x<sub>0</sub><sup>&#35;</sup>, &hellip;,
  49.  * x<sub>2N-1</sub><sup>&#35;</sup> <br>
  50.  * y<sub>n</sub> = (i / 2) &sum;<sub>k=0</sub><sup>2N-1</sup>
  51.  * x<sub>k</sub><sup>&#35;</sup> exp[-2&pi;i nk / (2N)]
  52.  * &nbsp;&nbsp;&nbsp;&nbsp;k = 0, &hellip;, N-1.
  53.  * <p>
  54.  * The present implementation of the discrete sine transform as a fast sine
  55.  * transform requires the length of the data to be a power of two. Besides,
  56.  * it implicitly assumes that the sampled function is odd. In particular, the
  57.  * first element of the data set must be 0, which is enforced in
  58.  * {@link #apply(DoubleUnaryOperator, double, double, int)},
  59.  * after sampling.
  60.  */
  61. public class FastSineTransform implements RealTransform {
  62.     /** Operation to be performed. */
  63.     private final UnaryOperator<double[]> op;

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

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

  79.     /**
  80.      * {@inheritDoc}
  81.      *
  82.      * The first element of the specified data set is required to be {@code 0}.
  83.      *
  84.      * @throws IllegalArgumentException if the length of the data array is
  85.      * not a power of two, or the first element of the data array is not zero.
  86.      */
  87.     @Override
  88.     public double[] apply(final double[] f) {
  89.         return op.apply(f);
  90.     }

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

  109.     /**
  110.      * Perform the FST algorithm (including inverse).
  111.      * The first element of the data set is required to be {@code 0}.
  112.      *
  113.      * @param f Data array to be transformed.
  114.      * @return the transformed array.
  115.      * @throws IllegalArgumentException if the length of the data array is
  116.      * not a power of two, or the first element of the data array is not zero.
  117.      */
  118.     private double[] fst(double[] f) {
  119.         if (!ArithmeticUtils.isPowerOfTwo(f.length)) {
  120.             throw new TransformException(TransformException.NOT_POWER_OF_TWO,
  121.                                          f.length);
  122.         }
  123.         if (f[0] != 0) {
  124.             throw new TransformException(TransformException.FIRST_ELEMENT_NOT_ZERO,
  125.                                          f[0]);
  126.         }

  127.         final double[] transformed = new double[f.length];
  128.         final int n = f.length;
  129.         if (n == 1) {
  130.             transformed[0] = 0;
  131.             return transformed;
  132.         }

  133.         // construct a new array and perform FFT on it
  134.         final double[] x = new double[n];
  135.         x[0] = 0;
  136.         final int nShifted = n >> 1;
  137.         x[nShifted] = 2 * f[nShifted];
  138.         final double piOverN = Math.PI / n;
  139.         for (int i = 1; i < nShifted; i++) {
  140.             final int nMi = n - i;
  141.             final double fi = f[i];
  142.             final double fnMi = f[nMi];
  143.             final double a = Math.sin(i * piOverN) * (fi + fnMi);
  144.             final double b = 0.5 * (fi - fnMi);
  145.             x[i] = a + b;
  146.             x[nMi] = a - b;
  147.         }

  148.         final FastFourierTransform transform = new FastFourierTransform(FastFourierTransform.Norm.STD);
  149.         final Complex[] y = transform.apply(x);

  150.         // reconstruct the FST result for the original array
  151.         transformed[0] = 0;
  152.         transformed[1] = 0.5 * y[0].getReal();
  153.         for (int i = 1; i < nShifted; i++) {
  154.             final int i2 = 2 * i;
  155.             transformed[i2] = -y[i].getImaginary();
  156.             transformed[i2 + 1] = y[i].getReal() + transformed[i2 - 1];
  157.         }

  158.         return transformed;
  159.     }

  160.     /**
  161.      * Factory method.
  162.      *
  163.      * @param normalization Normalization to be applied to the
  164.      * transformed data.
  165.      * @param inverse Whether to perform the inverse transform.
  166.      * @return the transform operator.
  167.      */
  168.     private UnaryOperator<double[]> create(final Norm normalization,
  169.                                            final boolean inverse) {
  170.         if (inverse) {
  171.             return normalization == Norm.ORTHO ?
  172.                 f -> TransformUtils.scaleInPlace(fst(f), Math.sqrt(2d / f.length)) :
  173.                 f -> TransformUtils.scaleInPlace(fst(f), 2d / f.length);
  174.         } else {
  175.             return normalization == Norm.ORTHO ?
  176.                 f -> TransformUtils.scaleInPlace(fst(f), Math.sqrt(2d / f.length)) :
  177.                 f -> fst(f);
  178.         }
  179.     }

  180.     /**
  181.      * Normalization types.
  182.      */
  183.     public enum Norm {
  184.         /**
  185.          * Should be passed to the constructor of {@link FastSineTransform} to
  186.          * use the <em>standard</em> normalization convention. The standard DST-I
  187.          * normalization convention is defined as follows
  188.          * <ul>
  189.          * <li>forward transform: y<sub>n</sub> = &sum;<sub>k=0</sub><sup>N-1</sup>
  190.          * x<sub>k</sub> sin(&pi; nk / N),</li>
  191.          * <li>inverse transform: x<sub>k</sub> = (2 / N)
  192.          * &sum;<sub>n=0</sub><sup>N-1</sup> y<sub>n</sub> sin(&pi; nk / N),</li>
  193.          * </ul>
  194.          * where N is the size of the data sample, and x<sub>0</sub> = 0.
  195.          */
  196.         STD,

  197.         /**
  198.          * Should be passed to the constructor of {@link FastSineTransform} to
  199.          * use the <em>orthogonal</em> normalization convention. The orthogonal
  200.          * DCT-I normalization convention is defined as follows
  201.          * <ul>
  202.          * <li>Forward transform: y<sub>n</sub> = &radic;(2 / N)
  203.          * &sum;<sub>k=0</sub><sup>N-1</sup> x<sub>k</sub> sin(&pi; nk / N),</li>
  204.          * <li>Inverse transform: x<sub>k</sub> = &radic;(2 / N)
  205.          * &sum;<sub>n=0</sub><sup>N-1</sup> y<sub>n</sub> sin(&pi; nk / N),</li>
  206.          * </ul>
  207.          * which makes the transform orthogonal. N is the size of the data sample,
  208.          * and x<sub>0</sub> = 0.
  209.          */
  210.         ORTHO
  211.     }
  212. }