There are several variants of the discrete sine transform. The present implementation * corresponds to DST-I, with various normalization conventions, which are specified by the * parameter {@link DstNormalization}. It should be noted that regardless to the convention, * the first element of the dataset to be transformed must be zero. * *
DST-I is equivalent to DFT of an odd extension of the data series. More precisely, if * x0, …, xN-1 is the data set to be sine transformed, the extended * data set x0#, …, x2N-1# is defined as * follows * *
Then, the standard DST-I y0, …, yN-1 of the real data set
* x0, …, xN-1 is equal to half of i (the pure imaginary
* number) times the N first elements of the DFT of the extended data set
* x0#, …, x2N-1#
* yn = (i / 2) ∑k=02N-1 xk#
* exp[-2πi nk / (2N)] k = 0, …, N-1.
*
*
The present implementation of the discrete sine transform as a fast sine transform requires * the length of the data to be a power of two. Besides, it implicitly assumes that the sampled * function is odd. In particular, the first element of the data set must be 0, which is enforced in * {@link #transform(UnivariateFunction, double, double, int, TransformType)}, after sampling. * * @since 1.2 */ public class FastSineTransformer implements RealTransformer, Serializable { /** Serializable version identifier. */ static final long serialVersionUID = 20120211L; /** The type of DST to be performed. */ private final DstNormalization 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 FastSineTransformer(final DstNormalization normalization) { this.normalization = normalization; } /** * {@inheritDoc} * *
The first element of the specified data set is required to be {@code 0}. * * @throws MathIllegalArgumentException if the length of the data array is not a power of two, * or the first element of the data array is not zero */ public double[] transform(final double[] f, final TransformType type) { if (normalization == DstNormalization.ORTHOGONAL_DST_I) { final double s = FastMath.sqrt(2.0 / f.length); return TransformUtils.scaleArray(fst(f), s); } if (type == TransformType.FORWARD) { return fst(f); } final double s = 2.0 / f.length; return TransformUtils.scaleArray(fst(f), s); } /** * {@inheritDoc} * *
This implementation enforces {@code f(x) = 0.0} at {@code x = 0.0}. * * @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 */ public double[] transform( final UnivariateFunction f, final double min, final double max, final int n, final TransformType type) { final double[] data = FunctionUtils.sample(f, min, max, n); data[0] = 0.0; return transform(data, type); } /** * Perform the FST algorithm (including inverse). The first element of the data set is required * to be {@code 0}. * * @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, * or the first element of the data array is not zero */ protected double[] fst(double[] f) throws MathIllegalArgumentException { final double[] transformed = new double[f.length]; if (!ArithmeticUtils.isPowerOfTwo(f.length)) { throw new MathIllegalArgumentException( LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING, Integer.valueOf(f.length)); } if (f[0] != 0.0) { throw new MathIllegalArgumentException( LocalizedFormats.FIRST_ELEMENT_NOT_ZERO, Double.valueOf(f[0])); } final int n = f.length; if (n == 1) { // trivial case transformed[0] = 0.0; return transformed; } // construct a new array and perform FFT on it final double[] x = new double[n]; x[0] = 0.0; x[n >> 1] = 2.0 * f[n >> 1]; for (int i = 1; i < (n >> 1); i++) { final double a = FastMath.sin(i * FastMath.PI / n) * (f[i] + f[n - i]); final double b = 0.5 * (f[i] - f[n - i]); x[i] = a + b; x[n - i] = a - b; } FastFourierTransformer transformer; transformer = new FastFourierTransformer(DftNormalization.STANDARD); Complex[] y = transformer.transform(x, TransformType.FORWARD); // reconstruct the FST result for the original array transformed[0] = 0.0; transformed[1] = 0.5 * y[0].getReal(); for (int i = 1; i < (n >> 1); i++) { transformed[2 * i] = -y[i].getImaginary(); transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1]; } return transformed; } }