diff --git a/docs/api-reference.md b/docs/api-reference.md index a459743..1dc3dfc 100644 --- a/docs/api-reference.md +++ b/docs/api-reference.md @@ -10,4 +10,5 @@ cov expand_dims kron + sinc ``` diff --git a/src/array_api_extra/__init__.py b/src/array_api_extra/__init__.py index 7a46760..3049274 100644 --- a/src/array_api_extra/__init__.py +++ b/src/array_api_extra/__init__.py @@ -1,7 +1,7 @@ from __future__ import annotations -from ._funcs import atleast_nd, cov, expand_dims, kron +from ._funcs import atleast_nd, cov, expand_dims, kron, sinc __version__ = "0.1.2.dev0" -__all__ = ["__version__", "atleast_nd", "cov", "expand_dims", "kron"] +__all__ = ["__version__", "atleast_nd", "cov", "expand_dims", "kron", "sinc"] diff --git a/src/array_api_extra/_funcs.py b/src/array_api_extra/_funcs.py index a371768..459db63 100644 --- a/src/array_api_extra/_funcs.py +++ b/src/array_api_extra/_funcs.py @@ -6,7 +6,7 @@ if TYPE_CHECKING: from ._typing import Array, ModuleType -__all__ = ["atleast_nd", "cov", "expand_dims", "kron"] +__all__ = ["atleast_nd", "cov", "expand_dims", "kron", "sinc"] def atleast_nd(x: Array, /, *, ndim: int, xp: ModuleType) -> Array: @@ -348,3 +348,88 @@ def kron(a: Array, b: Array, /, *, xp: ModuleType) -> Array: a_shape = xp.asarray(a_shape) b_shape = xp.asarray(b_shape) return xp.reshape(result, tuple(xp.multiply(a_shape, b_shape))) + + +def sinc(x: Array, /, *, xp: ModuleType) -> Array: + r""" + Return the normalized sinc function. + + The sinc function is equal to :math:`\sin(\pi x)/(\pi x)` for any argument + :math:`x\ne 0`. ``sinc(0)`` takes the limit value 1, making ``sinc`` not + only everywhere continuous but also infinitely differentiable. + + .. note:: + + Note the normalization factor of ``pi`` used in the definition. + This is the most commonly used definition in signal processing. + Use ``sinc(x / xp.pi)`` to obtain the unnormalized sinc function + :math:`\sin(x)/x` that is more common in mathematics. + + Parameters + ---------- + x : array + Array (possibly multi-dimensional) of values for which to calculate + ``sinc(x)``. Must have a real floating point dtype. + xp : array_namespace + The standard-compatible namespace for `x`. + + Returns + ------- + res : array + ``sinc(x)`` calculated elementwise, which has the same shape as the input. + + Notes + ----- + The name sinc is short for "sine cardinal" or "sinus cardinalis". + + The sinc function is used in various signal processing applications, + including in anti-aliasing, in the construction of a Lanczos resampling + filter, and in interpolation. + + For bandlimited interpolation of discrete-time signals, the ideal + interpolation kernel is proportional to the sinc function. + + References + ---------- + .. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web + Resource. https://mathworld.wolfram.com/SincFunction.html + .. [2] Wikipedia, "Sinc function", + https://en.wikipedia.org/wiki/Sinc_function + + Examples + -------- + >>> import array_api_strict as xp + >>> import array_api_extra as xpx + >>> x = xp.linspace(-4, 4, 41) + >>> xpx.sinc(x, xp=xp) + Array([-3.89817183e-17, -4.92362781e-02, + -8.40918587e-02, -8.90384387e-02, + -5.84680802e-02, 3.89817183e-17, + 6.68206631e-02, 1.16434881e-01, + 1.26137788e-01, 8.50444803e-02, + -3.89817183e-17, -1.03943254e-01, + -1.89206682e-01, -2.16236208e-01, + -1.55914881e-01, 3.89817183e-17, + 2.33872321e-01, 5.04551152e-01, + 7.56826729e-01, 9.35489284e-01, + 1.00000000e+00, 9.35489284e-01, + 7.56826729e-01, 5.04551152e-01, + 2.33872321e-01, 3.89817183e-17, + -1.55914881e-01, -2.16236208e-01, + -1.89206682e-01, -1.03943254e-01, + -3.89817183e-17, 8.50444803e-02, + 1.26137788e-01, 1.16434881e-01, + 6.68206631e-02, 3.89817183e-17, + -5.84680802e-02, -8.90384387e-02, + -8.40918587e-02, -4.92362781e-02, + -3.89817183e-17], dtype=array_api_strict.float64) + + """ + if not xp.isdtype(x.dtype, "real floating"): + err_msg = "`x` must have a real floating data type." + raise ValueError(err_msg) + # no scalars in `where` - array-api#807 + y = xp.pi * xp.where( + x, x, xp.asarray(xp.finfo(x.dtype).smallest_normal, dtype=x.dtype) + ) + return xp.sin(y) / y diff --git a/tests/test_funcs.py b/tests/test_funcs.py index 556add1..30e71eb 100644 --- a/tests/test_funcs.py +++ b/tests/test_funcs.py @@ -9,7 +9,7 @@ import pytest from numpy.testing import assert_allclose, assert_array_equal, assert_equal -from array_api_extra import atleast_nd, cov, expand_dims, kron +from array_api_extra import atleast_nd, cov, expand_dims, kron, sinc if TYPE_CHECKING: Array = Any # To be changed to a Protocol later (see array-api#589) @@ -224,3 +224,22 @@ def test_positive_negative_repeated(self): a = xp.empty((2, 3, 4, 5)) with pytest.raises(ValueError, match="Duplicate dimensions"): expand_dims(a, axis=(3, -3), xp=xp) + + +class TestSinc: + def test_simple(self): + assert_array_equal(sinc(xp.asarray(0.0), xp=xp), xp.asarray(1.0)) + w = sinc(xp.linspace(-1, 1, 100), xp=xp) + # check symmetry + assert_allclose(w, xp.flip(w, axis=0)) + + @pytest.mark.parametrize("x", [0, 1 + 3j]) + def test_dtype(self, x): + with pytest.raises(ValueError, match="real floating data type"): + sinc(xp.asarray(x), xp=xp) + + def test_3d(self): + x = xp.reshape(xp.arange(18, dtype=xp.float64), (3, 3, 2)) + expected = xp.zeros((3, 3, 2)) + expected[0, 0, 0] = 1.0 + assert_allclose(sinc(x, xp=xp), expected, atol=1e-15)