001 /*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements. See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License. You may obtain a copy of the License at
008 *
009 * http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017 package org.apache.commons.math3.transform;
018
019 import java.io.Serializable;
020
021 import org.apache.commons.math3.analysis.FunctionUtils;
022 import org.apache.commons.math3.analysis.UnivariateFunction;
023 import org.apache.commons.math3.complex.Complex;
024 import org.apache.commons.math3.exception.MathIllegalArgumentException;
025 import org.apache.commons.math3.exception.util.LocalizedFormats;
026 import org.apache.commons.math3.util.ArithmeticUtils;
027 import org.apache.commons.math3.util.FastMath;
028
029 /**
030 * Implements the Fast Sine Transform for transformation of one-dimensional real
031 * data sets. For reference, see James S. Walker, <em>Fast Fourier
032 * Transforms</em>, chapter 3 (ISBN 0849371635).
033 * <p>
034 * There are several variants of the discrete sine transform. The present
035 * implementation corresponds to DST-I, with various normalization conventions,
036 * which are specified by the parameter {@link DstNormalization}.
037 * <strong>It should be noted that regardless to the convention, the first
038 * element of the dataset to be transformed must be zero.</strong>
039 * <p>
040 * DST-I is equivalent to DFT of an <em>odd extension</em> of the data series.
041 * More precisely, if x<sub>0</sub>, …, x<sub>N-1</sub> is the data set
042 * to be sine transformed, the extended data set x<sub>0</sub><sup>#</sup>,
043 * …, x<sub>2N-1</sub><sup>#</sup> is defined as follows
044 * <ul>
045 * <li>x<sub>0</sub><sup>#</sup> = x<sub>0</sub> = 0,</li>
046 * <li>x<sub>k</sub><sup>#</sup> = x<sub>k</sub> if 1 ≤ k < N,</li>
047 * <li>x<sub>N</sub><sup>#</sup> = 0,</li>
048 * <li>x<sub>k</sub><sup>#</sup> = -x<sub>2N-k</sub> if N + 1 ≤ k <
049 * 2N.</li>
050 * </ul>
051 * <p>
052 * Then, the standard DST-I y<sub>0</sub>, …, y<sub>N-1</sub> of the real
053 * data set x<sub>0</sub>, …, x<sub>N-1</sub> is equal to <em>half</em>
054 * of i (the pure imaginary number) times the N first elements of the DFT of the
055 * extended data set x<sub>0</sub><sup>#</sup>, …,
056 * x<sub>2N-1</sub><sup>#</sup> <br />
057 * y<sub>n</sub> = (i / 2) ∑<sub>k=0</sub><sup>2N-1</sup>
058 * x<sub>k</sub><sup>#</sup> exp[-2πi nk / (2N)]
059 * k = 0, …, N-1.
060 * <p>
061 * The present implementation of the discrete sine transform as a fast sine
062 * transform requires the length of the data to be a power of two. Besides,
063 * it implicitly assumes that the sampled function is odd. In particular, the
064 * first element of the data set must be 0, which is enforced in
065 * {@link #transform(UnivariateFunction, double, double, int, TransformType)},
066 * after sampling.
067 *
068 * @version $Id: FastSineTransformer.java 1385310 2012-09-16 16:32:10Z tn $
069 * @since 1.2
070 */
071 public class FastSineTransformer implements RealTransformer, Serializable {
072
073 /** Serializable version identifier. */
074 static final long serialVersionUID = 20120211L;
075
076 /** The type of DST to be performed. */
077 private final DstNormalization normalization;
078
079 /**
080 * Creates a new instance of this class, with various normalization conventions.
081 *
082 * @param normalization the type of normalization to be applied to the transformed data
083 */
084 public FastSineTransformer(final DstNormalization normalization) {
085 this.normalization = normalization;
086 }
087
088 /**
089 * {@inheritDoc}
090 *
091 * The first element of the specified data set is required to be {@code 0}.
092 *
093 * @throws MathIllegalArgumentException if the length of the data array is
094 * not a power of two, or the first element of the data array is not zero
095 */
096 public double[] transform(final double[] f, final TransformType type) {
097 if (normalization == DstNormalization.ORTHOGONAL_DST_I) {
098 final double s = FastMath.sqrt(2.0 / f.length);
099 return TransformUtils.scaleArray(fst(f), s);
100 }
101 if (type == TransformType.FORWARD) {
102 return fst(f);
103 }
104 final double s = 2.0 / f.length;
105 return TransformUtils.scaleArray(fst(f), s);
106 }
107
108 /**
109 * {@inheritDoc}
110 *
111 * This implementation enforces {@code f(x) = 0.0} at {@code x = 0.0}.
112 *
113 * @throws org.apache.commons.math3.exception.NonMonotonicSequenceException
114 * if the lower bound is greater than, or equal to the upper bound
115 * @throws org.apache.commons.math3.exception.NotStrictlyPositiveException
116 * if the number of sample points is negative
117 * @throws MathIllegalArgumentException if the number of sample points is not a power of two
118 */
119 public double[] transform(final UnivariateFunction f,
120 final double min, final double max, final int n,
121 final TransformType type) {
122
123 final double[] data = FunctionUtils.sample(f, min, max, n);
124 data[0] = 0.0;
125 return transform(data, type);
126 }
127
128 /**
129 * Perform the FST algorithm (including inverse). The first element of the
130 * data set is required to be {@code 0}.
131 *
132 * @param f the real data array to be transformed
133 * @return the real transformed array
134 * @throws MathIllegalArgumentException if the length of the data array is
135 * not a power of two, or the first element of the data array is not zero
136 */
137 protected double[] fst(double[] f) throws MathIllegalArgumentException {
138
139 final double[] transformed = new double[f.length];
140
141 if (!ArithmeticUtils.isPowerOfTwo(f.length)) {
142 throw new MathIllegalArgumentException(
143 LocalizedFormats.NOT_POWER_OF_TWO_CONSIDER_PADDING,
144 Integer.valueOf(f.length));
145 }
146 if (f[0] != 0.0) {
147 throw new MathIllegalArgumentException(
148 LocalizedFormats.FIRST_ELEMENT_NOT_ZERO,
149 Double.valueOf(f[0]));
150 }
151 final int n = f.length;
152 if (n == 1) { // trivial case
153 transformed[0] = 0.0;
154 return transformed;
155 }
156
157 // construct a new array and perform FFT on it
158 final double[] x = new double[n];
159 x[0] = 0.0;
160 x[n >> 1] = 2.0 * f[n >> 1];
161 for (int i = 1; i < (n >> 1); i++) {
162 final double a = FastMath.sin(i * FastMath.PI / n) * (f[i] + f[n - i]);
163 final double b = 0.5 * (f[i] - f[n - i]);
164 x[i] = a + b;
165 x[n - i] = a - b;
166 }
167 FastFourierTransformer transformer;
168 transformer = new FastFourierTransformer(DftNormalization.STANDARD);
169 Complex[] y = transformer.transform(x, TransformType.FORWARD);
170
171 // reconstruct the FST result for the original array
172 transformed[0] = 0.0;
173 transformed[1] = 0.5 * y[0].getReal();
174 for (int i = 1; i < (n >> 1); i++) {
175 transformed[2 * i] = -y[i].getImaginary();
176 transformed[2 * i + 1] = y[i].getReal() + transformed[2 * i - 1];
177 }
178
179 return transformed;
180 }
181 }