All Downloads are FREE. Search and download functionalities are using the official Maven repository.

org.apache.commons.math3.transform.FastCosineTransformer Maven / Gradle / Ivy

There is a newer version: 2.12.15
Show newest version
/*
 * 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.math3.transform;

import java.io.Serializable;

import org.apache.commons.math3.analysis.FunctionUtils;
import org.apache.commons.math3.analysis.UnivariateFunction;
import org.apache.commons.math3.complex.Complex;
import org.apache.commons.math3.exception.MathIllegalArgumentException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.util.ArithmeticUtils;
import org.apache.commons.math3.util.FastMath;

/**
 * Implements the Fast Cosine Transform for transformation of one-dimensional
 * real data sets. For reference, see James S. Walker, Fast Fourier
 * Transforms, chapter 3 (ISBN 0849371635).
 * 

* 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 DctNormalization}. *

* DCT-I is equivalent to DFT of an even extension of the data series. * More precisely, if x0, …, xN-1 is the data set * to be cosine transformed, the extended data set * x0#, …, x2N-3# * is defined as follows *

    *
  • xk# = xk if 0 ≤ k < N,
  • *
  • xk# = x2N-2-k * if N ≤ k < 2N - 2.
  • *
*

* Then, the standard DCT-I y0, …, yN-1 of the real * data set x0, …, xN-1 is equal to half * of the N first elements of the DFT of the extended data set * x0#, …, x2N-3# *
* yn = (1 / 2) ∑k=02N-3 * xk# exp[-2πi nk / (2N - 2)] *     k = 0, …, N-1. *

* 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 = 2n + 1). Besides, it implicitly assumes * that the sampled function is even. * * @since 1.2 */ public class FastCosineTransformer implements RealTransformer, Serializable { /** Serializable version identifier. */ static final long serialVersionUID = 20120212L; /** The type of DCT to be performed. */ private final DctNormalization normalization; /** * Creates a new instance of this class, with various normalization * conventions. * * @param normalization the type of normalization to be applied to the * transformed data */ public FastCosineTransformer(final DctNormalization normalization) { this.normalization = normalization; } /** * {@inheritDoc} * * @throws MathIllegalArgumentException if the length of the data array is * not a power of two plus one */ public double[] transform(final double[] f, final TransformType type) throws MathIllegalArgumentException { if (type == TransformType.FORWARD) { if (normalization == DctNormalization.ORTHOGONAL_DCT_I) { final double s = FastMath.sqrt(2.0 / (f.length - 1)); return TransformUtils.scaleArray(fct(f), s); } return fct(f); } final double s2 = 2.0 / (f.length - 1); final double s1; if (normalization == DctNormalization.ORTHOGONAL_DCT_I) { s1 = FastMath.sqrt(s2); } else { s1 = s2; } return TransformUtils.scaleArray(fct(f), s1); } /** * {@inheritDoc} * * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException * if the lower bound is greater than, or equal to the upper bound * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException * if the number of sample points is negative * @throws MathIllegalArgumentException if the number of sample points is * not a power of two plus one */ public double[] transform(final UnivariateFunction f, final double min, final double max, final int n, final TransformType type) throws MathIllegalArgumentException { final double[] data = FunctionUtils.sample(f, min, max, n); return transform(data, type); } /** * Perform the FCT algorithm (including inverse). * * @param f the real data array to be transformed * @return the real transformed array * @throws MathIllegalArgumentException if the length of the data array is * not a power of two plus one */ protected double[] fct(double[] f) throws MathIllegalArgumentException { final double[] transformed = new double[f.length]; final int n = f.length - 1; if (!ArithmeticUtils.isPowerOfTwo(n)) { throw new MathIllegalArgumentException( LocalizedFormats.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]); x[n >> 1] = f[n >> 1]; // temporary variable for transformed[1] double t1 = 0.5 * (f[0] - f[n]); for (int i = 1; i < (n >> 1); i++) { final double a = 0.5 * (f[i] + f[n - i]); final double b = FastMath.sin(i * FastMath.PI / n) * (f[i] - f[n - i]); final double c = FastMath.cos(i * FastMath.PI / n) * (f[i] - f[n - i]); x[i] = a - b; x[n - i] = a + b; t1 += c; } FastFourierTransformer transformer; transformer = new FastFourierTransformer(DftNormalization.STANDARD); Complex[] y = transformer.transform(x, TransformType.FORWARD); // reconstruct the FCT result for the original array transformed[0] = y[0].getReal(); transformed[1] = t1; for (int i = 1; i < (n >> 1); i++) { transformed[2 * i] = y[i].getReal(); transformed[2 * i + 1] = transformed[2 * i - 1] - y[i].getImaginary(); } transformed[n] = y[n >> 1].getReal(); return transformed; } }





© 2015 - 2024 Weber Informatics LLC | Privacy Policy