Merge pull request #8101 from asi1024/convolve1d2o

Add `cupyx.signal.convolve1d2o`

cupyx/signal/__init__.py

@@ -1,5 +1,6 @@
 from cupyx.signal._acoustics import complex_cepstrum, real_cepstrum  # NOQA
 from cupyx.signal._acoustics import inverse_complex_cepstrum  # NOQA
 from cupyx.signal._acoustics import minimum_phase  # NOQA
+from cupyx.signal._convolution import convolve1d2o  # NOQA
 from cupyx.signal._convolution import convolve1d3o  # NOQA
 from cupyx.signal._radartools import pulse_compression, pulse_doppler, cfar_alpha  # NOQA

cupyx/signal/_convolution/__init__.py

@@ -1 +1,2 @@
+from cupyx.signal._convolution._convolve import convolve1d2o  # NOQA
 from cupyx.signal._convolution._convolve import convolve1d3o  # NOQA

cupyx/signal/_convolution/_convolve.py

@@ -26,6 +26,109 @@
 import cupy
 
 
+_convolve1d2o_kernel = cupy.ElementwiseKernel(
+    'raw T in1, raw T in2, int32 W, int32 H', 'T out',
+    """
+    T temp {};
+    for (int x = 0; x < W; x++) {
+      for (int y = 0; y < H; y++) {
+        temp += in1[i + W - x - 1] * in1[i + H - y - 1] * in2[H * x + y];
+      }
+    }
+    out = temp;
+    """,
+    "cupy_convolved2o",
+)
+
+
+def _convolve1d2o(in1, in2, mode):
+    assert mode == "valid"
+    out_dim = in1.shape[0] - max(in2.shape) + 1
+    dtype = cupy.result_type(in1, in2)
+    out = cupy.empty(out_dim, dtype=dtype)
+    _convolve1d2o_kernel(in1, in2, *in2.shape, out)
+    return out
+
+
+def convolve1d2o(in1, in2, mode='valid', method='direct'):
+    """
+    Convolve a 1-dimensional arrays with a 2nd order filter.
+    This results in a second order convolution.
+
+    Convolve `in1` and `in2`, with the output size determined by the
+    `mode` argument.
+
+    Parameters
+    ----------
+    in1 : array_like
+        First input.
+    in2 : array_like
+        Second input. Should have the same number of dimensions as `in1`.
+    mode : str {'full', 'valid', 'same'}, optional
+        A string indicating the size of the output:
+
+        ``full``
+           The output is the full discrete linear convolution
+           of the inputs. (Default)
+        ``valid``
+           The output consists only of those elements that do not
+           rely on the zero-padding. In 'valid' mode, either `in1` or `in2`
+           must be at least as large as the other in every dimension.
+        ``same``
+           The output is the same size as `in1`, centered
+           with respect to the 'full' output.
+    method : str {'auto', 'direct', 'fft'}, optional
+        A string indicating which method to use to calculate the convolution.
+
+        ``direct``
+           The convolution is determined directly from sums, the definition of
+           convolution.
+        ``fft``
+           The Fourier Transform is used to perform the convolution by calling
+           `fftconvolve`.
+        ``auto``
+           Automatically chooses direct or Fourier method based on an estimate
+           of which is faster (default).
+
+    Returns
+    -------
+    out : ndarray
+        A 1-dimensional array containing a subset of the discrete linear
+        convolution of `in1` with `in2`.
+
+    See Also
+    --------
+    convolve
+    convolve1d2o
+    convolve1d3o
+
+    Examples
+    --------
+    Convolution of a 2nd order filter on a 1d signal
+
+    >>> import cusignal as cs
+    >>> import numpy as np
+    >>> d = 50
+    >>> a = np.random.uniform(-1,1,(200))
+    >>> b = np.random.uniform(-1,1,(d,d))
+    >>> c = cs.convolve1d2o(a,b)
+
+    """
+
+    if in1.ndim != 1:
+        raise ValueError('in1 should have one dimension')
+    if in2.ndim != 2:
+        raise ValueError('in2 should have three dimension')
+
+    if mode in ["same", "full"]:
+        raise NotImplementedError("Mode == {} not implemented".format(mode))
+
+    if method == "direct":
+        return _convolve1d2o(in1, in2, mode)
+    else:
+        raise NotImplementedError("Only Direct method implemented")
+
+
 _convolve1d3o_kernel = cupy.ElementwiseKernel(
     'raw T in1, raw T in2, int32 W, int32 H, int32 D', 'T out',
     """

tests/cupyx_tests/signal_tests/convolution_tests/test_convolve.py

@@ -6,20 +6,59 @@
 from cupyx import signal
 
 
-def _convolve1d3o(in1, in2):
-    dtype = in1.dtype
-    W, H, D = in2.shape
-    size = in1.shape[0] - max(W, H, D) + 1
-    s = numpy.dtype(dtype).itemsize
-    from numpy.lib.stride_tricks import as_strided
-    X = as_strided(in1, (size, W), (s, s))[:, ::-1]
-    Y = as_strided(in1, (size, H), (s, s))[:, ::-1]
-    Z = as_strided(in1, (size, D), (s, s))[:, ::-1]
-    return numpy.einsum('ix,iy,iz,xyz->i', X, Y, Z, in2)
+class TestConvolve1d2o:
+
+    def _convolve1d2o(self, in1, in2):
+        dtype = in1.dtype
+        W, H = in2.shape
+        size = in1.shape[0] - max(W, H) + 1
+        s = numpy.dtype(dtype).itemsize
+        from numpy.lib.stride_tricks import as_strided
+        X = as_strided(in1, (size, W), (s, s))[:, ::-1]
+        Y = as_strided(in1, (size, H), (s, s))[:, ::-1]
+        return numpy.einsum('ix,iy,xy->i', X, Y, in2)
+
+    @testing.for_float_dtypes(no_float16=True)
+    @testing.numpy_cupy_allclose(rtol=2e-3)
+    @pytest.mark.parametrize('shape', [(50, 50), (40, 60)])
+    def test_convolve1d2o(self, dtype, xp, shape):
+        a = testing.shaped_random((200,), xp=xp, dtype=dtype, scale=2) - 1
+        b = testing.shaped_random(shape, xp=xp, dtype=dtype, scale=2) - 1
+        if xp is cupy:
+            return signal.convolve1d2o(a, b)
+        else:
+            assert xp is numpy
+            return self._convolve1d2o(a, b)
+
+    @testing.for_complex_dtypes()
+    @testing.numpy_cupy_allclose(rtol=2e-3)
+    @pytest.mark.parametrize('shape', [(50, 50), (40, 60)])
+    def test_convolve1d2o_complex(self, dtype, xp, shape):
+        # Just check that we can call the function
+        a = testing.shaped_random(
+            (200,), xp=xp, dtype=dtype, scale=2) - (1 + 1j)
+        b = testing.shaped_random(
+            shape, xp=xp, dtype=dtype, scale=2) - (1 + 1j)
+        if xp is cupy:
+            return signal.convolve1d2o(a, b)
+        else:
+            assert xp is numpy
+            return self._convolve1d2o(a, b)
 
 
 class TestConvolve1d3o:
 
+    def _convolve1d3o(self, in1, in2):
+        dtype = in1.dtype
+        W, H, D = in2.shape
+        size = in1.shape[0] - max(W, H, D) + 1
+        s = numpy.dtype(dtype).itemsize
+        from numpy.lib.stride_tricks import as_strided
+        X = as_strided(in1, (size, W), (s, s))[:, ::-1]
+        Y = as_strided(in1, (size, H), (s, s))[:, ::-1]
+        Z = as_strided(in1, (size, D), (s, s))[:, ::-1]
+        return numpy.einsum('ix,iy,iz,xyz->i', X, Y, Z, in2)
+
     @testing.for_float_dtypes(no_float16=True)
     @testing.numpy_cupy_allclose(rtol=2e-3)
     @pytest.mark.parametrize('shape', [(50, 50, 50), (40, 50, 60)])
@@ -30,7 +69,7 @@ def test_convolve1d3o(self, dtype, xp, shape):
             return signal.convolve1d3o(a, b)
         else:
             assert xp is numpy
-            return _convolve1d3o(a, b)
+            return self._convolve1d3o(a, b)
 
     @testing.for_complex_dtypes()
     @testing.numpy_cupy_allclose(rtol=2e-3)
@@ -45,4 +84,4 @@ def test_convolve1d3o_complex(self, dtype, xp, shape):
             return signal.convolve1d3o(a, b)
         else:
             assert xp is numpy
-            return _convolve1d3o(a, b)
+            return self._convolve1d3o(a, b)