diff options
Diffstat (limited to 'src/main/java/org/apache/commons/math3/analysis/function/Sinc.java')
| -rw-r--r-- | src/main/java/org/apache/commons/math3/analysis/function/Sinc.java | 205 |
1 files changed, 205 insertions, 0 deletions
diff --git a/src/main/java/org/apache/commons/math3/analysis/function/Sinc.java b/src/main/java/org/apache/commons/math3/analysis/function/Sinc.java new file mode 100644 index 0000000..553cfff --- /dev/null +++ b/src/main/java/org/apache/commons/math3/analysis/function/Sinc.java @@ -0,0 +1,205 @@ +/* + * 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.analysis.function; + +import org.apache.commons.math3.analysis.DifferentiableUnivariateFunction; +import org.apache.commons.math3.analysis.FunctionUtils; +import org.apache.commons.math3.analysis.UnivariateFunction; +import org.apache.commons.math3.analysis.differentiation.DerivativeStructure; +import org.apache.commons.math3.analysis.differentiation.UnivariateDifferentiableFunction; +import org.apache.commons.math3.exception.DimensionMismatchException; +import org.apache.commons.math3.util.FastMath; + +/** + * <a href="http://en.wikipedia.org/wiki/Sinc_function">Sinc</a> function, + * defined by + * <pre><code> + * sinc(x) = 1 if x = 0, + * sin(x) / x otherwise. + * </code></pre> + * + * @since 3.0 + */ +public class Sinc implements UnivariateDifferentiableFunction, DifferentiableUnivariateFunction { + /** + * Value below which the computations are done using Taylor series. + * <p> + * The Taylor series for sinc even order derivatives are: + * <pre> + * d^(2n)sinc/dx^(2n) = Sum_(k>=0) (-1)^(n+k) / ((2k)!(2n+2k+1)) x^(2k) + * = (-1)^n [ 1/(2n+1) - x^2/(4n+6) + x^4/(48n+120) - x^6/(1440n+5040) + O(x^8) ] + * </pre> + * </p> + * <p> + * The Taylor series for sinc odd order derivatives are: + * <pre> + * d^(2n+1)sinc/dx^(2n+1) = Sum_(k>=0) (-1)^(n+k+1) / ((2k+1)!(2n+2k+3)) x^(2k+1) + * = (-1)^(n+1) [ x/(2n+3) - x^3/(12n+30) + x^5/(240n+840) - x^7/(10080n+45360) + O(x^9) ] + * </pre> + * </p> + * <p> + * So the ratio of the fourth term with respect to the first term + * is always smaller than x^6/720, for all derivative orders. + * This implies that neglecting this term and using only the first three terms induces + * a relative error bounded by x^6/720. The SHORTCUT value is chosen such that this + * relative error is below double precision accuracy when |x| <= SHORTCUT. + * </p> + */ + private static final double SHORTCUT = 6.0e-3; + /** For normalized sinc function. */ + private final boolean normalized; + + /** + * The sinc function, {@code sin(x) / x}. + */ + public Sinc() { + this(false); + } + + /** + * Instantiates the sinc function. + * + * @param normalized If {@code true}, the function is + * <code> sin(πx) / πx</code>, otherwise {@code sin(x) / x}. + */ + public Sinc(boolean normalized) { + this.normalized = normalized; + } + + /** {@inheritDoc} */ + public double value(final double x) { + final double scaledX = normalized ? FastMath.PI * x : x; + if (FastMath.abs(scaledX) <= SHORTCUT) { + // use Taylor series + final double scaledX2 = scaledX * scaledX; + return ((scaledX2 - 20) * scaledX2 + 120) / 120; + } else { + // use definition expression + return FastMath.sin(scaledX) / scaledX; + } + } + + /** {@inheritDoc} + * @deprecated as of 3.1, replaced by {@link #value(DerivativeStructure)} + */ + @Deprecated + public UnivariateFunction derivative() { + return FunctionUtils.toDifferentiableUnivariateFunction(this).derivative(); + } + + /** {@inheritDoc} + * @since 3.1 + */ + public DerivativeStructure value(final DerivativeStructure t) + throws DimensionMismatchException { + + final double scaledX = (normalized ? FastMath.PI : 1) * t.getValue(); + final double scaledX2 = scaledX * scaledX; + + double[] f = new double[t.getOrder() + 1]; + + if (FastMath.abs(scaledX) <= SHORTCUT) { + + for (int i = 0; i < f.length; ++i) { + final int k = i / 2; + if ((i & 0x1) == 0) { + // even derivation order + f[i] = (((k & 0x1) == 0) ? 1 : -1) * + (1.0 / (i + 1) - scaledX2 * (1.0 / (2 * i + 6) - scaledX2 / (24 * i + 120))); + } else { + // odd derivation order + f[i] = (((k & 0x1) == 0) ? -scaledX : scaledX) * + (1.0 / (i + 2) - scaledX2 * (1.0 / (6 * i + 24) - scaledX2 / (120 * i + 720))); + } + } + + } else { + + final double inv = 1 / scaledX; + final double cos = FastMath.cos(scaledX); + final double sin = FastMath.sin(scaledX); + + f[0] = inv * sin; + + // the nth order derivative of sinc has the form: + // dn(sinc(x)/dxn = [S_n(x) sin(x) + C_n(x) cos(x)] / x^(n+1) + // where S_n(x) is an even polynomial with degree n-1 or n (depending on parity) + // and C_n(x) is an odd polynomial with degree n-1 or n (depending on parity) + // S_0(x) = 1, S_1(x) = -1, S_2(x) = -x^2 + 2, S_3(x) = 3x^2 - 6... + // C_0(x) = 0, C_1(x) = x, C_2(x) = -2x, C_3(x) = -x^3 + 6x... + // the general recurrence relations for S_n and C_n are: + // S_n(x) = x S_(n-1)'(x) - n S_(n-1)(x) - x C_(n-1)(x) + // C_n(x) = x C_(n-1)'(x) - n C_(n-1)(x) + x S_(n-1)(x) + // as per polynomials parity, we can store both S_n and C_n in the same array + final double[] sc = new double[f.length]; + sc[0] = 1; + + double coeff = inv; + for (int n = 1; n < f.length; ++n) { + + double s = 0; + double c = 0; + + // update and evaluate polynomials S_n(x) and C_n(x) + final int kStart; + if ((n & 0x1) == 0) { + // even derivation order, S_n is degree n and C_n is degree n-1 + sc[n] = 0; + kStart = n; + } else { + // odd derivation order, S_n is degree n-1 and C_n is degree n + sc[n] = sc[n - 1]; + c = sc[n]; + kStart = n - 1; + } + + // in this loop, k is always even + for (int k = kStart; k > 1; k -= 2) { + + // sine part + sc[k] = (k - n) * sc[k] - sc[k - 1]; + s = s * scaledX2 + sc[k]; + + // cosine part + sc[k - 1] = (k - 1 - n) * sc[k - 1] + sc[k -2]; + c = c * scaledX2 + sc[k - 1]; + + } + sc[0] *= -n; + s = s * scaledX2 + sc[0]; + + coeff *= inv; + f[n] = coeff * (s * sin + c * scaledX * cos); + + } + + } + + if (normalized) { + double scale = FastMath.PI; + for (int i = 1; i < f.length; ++i) { + f[i] *= scale; + scale *= FastMath.PI; + } + } + + return t.compose(f); + + } + +} |