Support complex filters in Convolve1D adjoints

examples/plot_convolve.py

@@ -59,6 +59,25 @@
 ax.set_xlim(1.9, 2.1)
 plt.tight_layout()
 
+###############################################################################
+# We repeat the same example with a complex-valued filter
+Cop = pylops.signalprocessing.Convolve1D(
+    nt, h=h + 1j * h, offset=hcenter, dtype="complex128"
+)
+y = Cop * x
+xinv = Cop / y
+
+fig, axs = plt.subplots(1, 2, figsize=(10, 3))
+for ax, component, title in zip(axs, (np.real, np.imag), ("Real", "Imaginary"), strict=True):
+    ax.plot(t, component(x), "k", lw=2, label=r"$x$")
+    ax.plot(t, component(y), "r", lw=2, label=r"$y=Ax$")
+    ax.plot(t, component(xinv), "--g", lw=2, label=r"$x_{ext}$")
+    ax.set_title(title, fontsize=14, fontweight="bold")
+    ax.legend()
+    ax.set_xlim(1.9, 2.1)
+fig.suptitle("Convolve 1d data with complex filter")
+plt.tight_layout()
+
 ###############################################################################
 # We show now that also a filter with mixed phase (i.e., not centered
 # around zero) can be applied and inverted for using the

pylops/signalprocessing/convolve1d.py

@@ -160,7 +160,7 @@ def __init__(
         # ``flip`` returns a view with negative strides, which the direct
         # convolution of some backends (notably cupy) reads incorrectly, so the
         # flipped filter is made contiguous before it reaches ``convfunc``
-        self.hstar = ncp.ascontiguousarray(ncp.flip(self.h, axis=haxis))
+        self.hstar = ncp.ascontiguousarray(ncp.flip(self.h, axis=haxis).conj())
 
         # add dimensions to filter to match dimensions of model and data
         if self.h.ndim == 1:
@@ -232,7 +232,7 @@ def __init__(
         # ``flip`` returns a view with negative strides, which the direct
         # convolution of some backends (notably cupy) reads incorrectly, so the
         # flipped filter is made contiguous before it reaches ``convfunc``
-        self.hstar = ncp.ascontiguousarray(ncp.flip(self.h, axis=haxis))
+        self.hstar = ncp.ascontiguousarray(ncp.flip(self.h, axis=haxis).conj())
 
         self.pad = np.zeros((len(dims), 2), dtype=int)
         self.pad[self.axis, 0] = max(self.offset, 0)
@@ -328,7 +328,7 @@ class Convolve1D(LinearOperator):
     nh : :obj:`int`
         Length of the filter
     hstar : :obj:`numpy.ndarray`
-        Time-reversed filter used in adjoint
+        Conjugate time-reversed filter used in adjoint
     convfunc : :obj:`callable`
         Function handler used to perform convolution
     dims : :obj:`tuple`
@@ -391,12 +391,12 @@ class Convolve1D(LinearOperator):
     In time domain:
 
     .. math::
-        x(t) = \int\limits_{-\infty}^{\infty} h(t+\tau) x(\tau) \,\mathrm{d}\tau
+        x(t) = \int\limits_{-\infty}^{\infty} h(\tau)^* y(t+\tau) \,\mathrm{d}\tau
 
     or in frequency domain:
 
     .. math::
-        y(t) = \mathscr{F}^{-1} (H(f)^* * X(f))
+        x(t) = \mathscr{F}^{-1} (H(f)^* * Y(f))
 
     """
 

pytests/test_convolve.py

@@ -219,7 +219,7 @@ def test_Convolve1D_hstar_contiguous(dims, hlen):
         (par2_1d),
     ],
 )
-@pytest.mark.parametrize("dtype", [np.float32, np.float64])
+@pytest.mark.parametrize("dtype", [np.float32, np.float64, np.complex64, np.complex128])
 @pytest.mark.parametrize("method", methods_1d)
 @pytest.mark.parametrize("hlen", ["short", "long"])
 def test_Convolve1D_1d(par, dtype, method, hlen):
@@ -238,13 +238,17 @@ def test_Convolve1D_1d(par, dtype, method, hlen):
     else:
         nx, x, h, offset = nfilt[0], h1.astype(dtype), spike, nfilt[0] // 2
 
+    if np.issubdtype(dtype, np.complexfloating):
+        h = h + 1j * h
+
     # Operator
     Cop = Convolve1D(nx, h=h, offset=offset, method=method, dtype=dtype)
     assert dottest(
         Cop,
         par["nx"],
         nx,
-        rtol=1e-4 if dtype == np.float32 else 1e-6,
+        rtol=1e-4 if dtype in (np.float32, np.complex64) else 1e-6,
+        complexflag=3 if np.issubdtype(dtype, np.complexfloating) else 0,
         backend=backend,
     )