diff --git a/examples/plot_convolve.py b/examples/plot_convolve.py index ae4e6c78..65c65bfb 100644 --- a/examples/plot_convolve.py +++ b/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 diff --git a/pylops/signalprocessing/convolve1d.py b/pylops/signalprocessing/convolve1d.py index ba4084d8..965d8712 100644 --- a/pylops/signalprocessing/convolve1d.py +++ b/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)) """ diff --git a/pytests/test_convolve.py b/pytests/test_convolve.py index b68e05d3..81947735 100644 --- a/pytests/test_convolve.py +++ b/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, )