8000 ENH: Added optional `scale` parameter to `sinc` to enable non-normali… · numpy/numpy@f62d50d · GitHub
[go: up one dir, main page]
Skip to content

Commit f62d50d

Browse files
madphysicistmadphysicist
committed
ENH: Added optional scale parameter to sinc to enable non-normalized function.
Currently, `sinc(x)` is defined as `sin(pi*x)/(pi*x)`. It is sometimes convenient to have the mathematical definition available: `sin(x)/(x). In fact, any scaling is now possible with a scalar input. Tests included. Also modified the docs for `where` very slightly.
1 parent aa746c6 commit f62d50d

File tree

3 files changed

+94
-14
lines changed

3 files changed

+94
-14
lines changed

‎doc/release/1.12.0-notes.rst

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -459,6 +459,15 @@ improved precision. This should be useful in packages such as Theano
459459
where the precision of float16 is adequate and its smaller footprint is
460460
desirable.
461461

462+
``sinc`` now accepts optional `scale` argument
463+
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
464+
Scale is the factor applied to `x` before passing it to the mathematical
465+
sinc function. The previous (default) scale was `np.pi`. Two string
466+
arguments are recognized in addition to any array_like that broadcasts
467+
to the shape of the input array: 'normalized' and 'unnormalized'. The
468+
current behavior is 'normalized' by default: sinc(x) = sin(pi*x)/(pi*x).
469+
The mathematical function is 'unnormalized': sinc(x) = sin(x)/x.
470+
462471

463472
Changes
464473
=======

‎numpy/lib/function_base.py

Lines changed: 45 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -3834,17 +3834,33 @@ def kaiser(M, beta):
38343834
return i0(beta * sqrt(1-((n-alpha)/alpha)**2.0))/i0(float(beta))
38353835

38363836

3837-
def sinc(x):
3838-
"""
3839-
Return the sinc function.
3837+
def sinc(x, scale = 'normalized'):
3838+
r"""
3839+
Apply the sinc function to an array, element-wise.
3840+
3841+
The normalized sinc function is :math:`\frac{\sin(\pi x)}{\pi x}`.
3842+
The unnormalized sinc function is :math:`\frac{\sin(x)}{x}`. The
3843+
normalized version is used more commonly because it has roots at
3844+
integer values of `x` and the integral is exactly one.
38403845
3841-
The sinc function is :math:`\\sin(\\pi x)/(\\pi x)`.
3846+
For arbitrary scales, the function is defined as
3847+
:math:`\frac{\sin(scale \cdot x)}{scale \cdot x}`. In this case, the
3848+
integral from negative infinity to infinity is
3849+
:math:`\frac{\pi}{scale}` and roots are at multiples of
3850+
`\frac{\pi}{scale}`.
38423851
38433852
Parameters
38443853
----------
38453854
x : ndarray
38463855
Array (possibly multi-dimensional) of values for which to to
38473856
calculate ``sinc(x)``.
3857+
scale : {'normalized', 'unnormalized', scalar or array_like}, optional
3858+
The scale to apply to `x` when computing the sinc function. The
3859+
default is the most common engineering usage, 'normalized'. The
3860+
mathematical definition is 'unnormalized'. If an arbitrary scale
3861+
is provided as an array, it must broadcast against `x`. This
3862+
parameter should be used as a keyword only parameter. Support for
3863+
it as a positional parameter may be removed in the near future.
38483864
38493865
Returns
38503866
-------
@@ -3858,16 +3874,16 @@ def sinc(x):
38583874
The name sinc is short for "sine cardinal" or "sinus cardinalis".
38593875
38603876
The sinc function is used in various signal processing applications,
3861-
including in anti-aliasing, in the construction of a Lanczos resampling
3862-
filter, and in interpolation.
3877+
including in anti-aliasing, in the construction of a Lanczos
3878+
resampling filter, and in interpolation.
38633879
38643880
For bandlimited interpolation of discrete-time signals, the ideal
38653881
interpolation kernel is proportional to the sinc function.
38663882
38673883
References
38683884
----------
3869-
.. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram Web
3870-
Resource. http://mathworld.wolfram.com/SincFunction.html
3885+
.. [1] Weisstein, Eric W. "Sinc Function." From MathWorld--A Wolfram
3886+
Web Resource. http://mathworld.wolfram.com/SincFunction.html
38713887
.. [2] Wikipedia, "Sinc function",
38723888
http://en.wikipedia.org/wiki/Sinc_function
38733889
@@ -3890,9 +3906,9 @@ def sinc(x):
38903906
-5.84680802e-02, -8.90384387e-02, -8.40918587e-02,
38913907
-4.92362781e-02, -3.89804309e-17])
38923908
3893-
>>> plt.plot(x, np.sinc(x))
3894-
[<matplotlib.lines.Line2D object at 0x...>]
3895-
>>> plt.title("Sinc Function")
3909+
>>> plt.plot(x, np.sinc(x), x, np.sinc(x, scale='unnormalized'))
3910+
[<matplotlib.lines.Line2D object at 0x...>, <matplotlib.lines.Line2D object at 0x...>]
3911+
>>> plt.title("Normalized and Un-normalized Sinc Functions")
38963912
<matplotlib.text.Text object at 0x...>
38973913
>>> plt.ylabel("Amplitude")
38983914
<matplotlib.text.Text object at 0x...>
@@ -3908,9 +3924,25 @@ def sinc(x):
39083924
<matplotlib.image.AxesImage object 67E6 at 0x...>
39093925
39103926
"""
3927+
if isinstance(scale, str):
3928+
if scale == 'normalized':
3929+
scale = pi
3930+
elif scale == 'unnormalized':
3931+
scale = 1.0
3932+
else:
3933+
raise ValueError("{} is not a valid scaling for `sinc`".format(scale))
3934+
39113935
x = np.asanyarray(x)
3912-
y = pi * where(x == 0, 1.0e-20, x)
3913-
return sin(y)/y
3936+
3937+
if isinstance(x.dtype, np.inexact):
3938+
min_val = np.finfo(x.dtype).tiny
3939+
else:
3940+
min_val = np.finfo(np.float_).tiny
3941+
3942+
y = where(x == 0, min_val, x)
3943+
if np.any(scale != 1.0):
3944+
y *= scale
3945+
return sin(y) / y
39143946

39153947

39163948
def msort(a):

‎numpy/lib/tests/test_function_base.py

Lines changed: 40 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1528,7 +1528,7 @@ def test_matrix(self):
15281528
class TestSinc(TestCase):
15291529

15301530
def test_simple(self):
1531-
assert_(sinc(0) == 1)
1531+
assert_equal(sinc(0), 1)
15321532
w = sinc(np.linspace(-1, 1, 100))
15331533
# check symmetry
15341534
assert_array_almost_equal(w, flipud(w), 7)
@@ -1541,6 +1541,45 @@ def test_array_like(self):
15411541
assert_array_equal(y1, y2)
15421542
assert_array_equal(y1, y3)
15431543

1544+
def test_scale(self):
1545+
x = np.linspace(-1, 1, 100)
1546+
1547+
# Check backwards-compatibility of default
1548+
assert_array_equal(sinc(x, scale='normalized'), sinc(x))
1549+
1550+
# Check unnormalized
1551+
assert_equal(sinc(0, scale='unnormalized'), 1)
1552+
assert_almost_equal(sinc(np.pi, scale='unnormalized'), 0)
1553+
assert_almost_equal(sinc(np.pi / 2.0, scale='unnormalized'), 2.0 / np.pi)
1554+
1555+
# Check scalar
1556+
assert_equal(sinc(0, scale=2), 1)
1557+
assert_almost_equal(sinc(np.pi / 2.0, scale=2.0), 0)
1558+
assert_almost_equal(sinc(np.pi / 4.0, scale=2.0), 2.0 / np.pi)
1559+
1560+
# Check negative
1561+
assert_equal(sinc(0, scale=-5), 1)
1562+
assert_equal(sinc(x, scale=-5.0), sinc(x, scale=5.0))
1563+
1564+
# Test huge
1565+
assert_equal(sinc(0, scale=1e20), 1)
1566+
1567+
# Test broadcast
1568+
y = np.arange(3 * 4).reshape(3, 4)
1569+
np.random.shuffle(y.ravel())
1570+
scales = np.arange(4) + 1
1571+
out = sinc(y, scale=scales)
1572+
assert_equal(out.shape, y.shape)
1573+
for i in range(scales.shape[0]):
1574+
assert_almost_equal(sinc(y[:, i], scale=scales[i]), out[:, i])
1575+
# Same size as input
1576+
assert_array_equal(sinc(y + 1, scale=1. / (y + 1)), np.sin(1))
1577+
# Incorrect broadcast
1578+
assert_raises(ValueError, sinc, y, scale=scales[0:-1])
1579+
1580+
# Check bad names
1581+
assert_raises(ValueError, sinc, x, scale='foo')
1582+
15441583

15451584
class TestHistogram(TestCase):
15461585

0 commit comments

Comments
 (0)
0