diff --git a/CHANGELOG.md b/CHANGELOG.md
index 09d6c06b..1f1c6780 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,7 @@ Changelog
 =========
 
 # 2.3.1
-* Fixed bug in :py:mod:`pylops.utils.backend` (see [Issue #606](https://github.com/PyLops/pylops/issues/606))
+* Fixed bug in `pylops.utils.backend` (see [Issue #606](https://github.com/PyLops/pylops/issues/606))
 
 # 2.3.0
 
diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst
index 1d77dc15..7b540e90 100755
--- a/docs/source/api/index.rst
+++ b/docs/source/api/index.rst
@@ -109,6 +109,7 @@ Signal processing
     Seislet
     Radon2D
     Radon3D
+    FourierRadon2D
     ChirpRadon2D
     ChirpRadon3D
     Sliding1D
diff --git a/examples/plot_fourierradon.py b/examples/plot_fourierradon.py
new file mode 100644
index 00000000..c59d21bb
--- /dev/null
+++ b/examples/plot_fourierradon.py
@@ -0,0 +1,130 @@
+r"""
+Fourier Radon Transform
+=======================
+This example shows how to use the :py:class:`pylops.signalprocessing.FourierRadon2D`
+operator to apply the linear and parabolic Radon Transform to 2-dimensional signals.
+
+This operator provides a transformation equivalent to that of
+:py:class:`pylops.signalprocessing.Radon2D`, however since the shift-and-sum step
+is performed in the frequency domain, this is analytically correct (compared to
+performing to shifting the data via nearest or linear interpolation).
+
+"""
+import matplotlib.pyplot as plt
+import numpy as np
+
+import pylops
+
+plt.close("all")
+
+###############################################################################
+# Let's start by creating a empty 2d matrix of size :math:`n_x \times n_t`
+# with a single linear event.
+
+par = {
+    "ot": 0,
+    "dt": 0.004,
+    "nt": 51,
+    "ox": -250,
+    "dx": 10,
+    "nx": 51,
+    "oy": -250,
+    "dy": 10,
+    "ny": 51,
+    "f0": 40,
+}
+theta = [10.0]
+t0 = [0.1]
+amp = [1.0]
+
+# Create axes
+t, t2, x, y = pylops.utils.seismicevents.makeaxis(par)
+dt, dx, dy = par["dt"], par["dx"], par["dy"]
+
+# Create wavelet
+wav, _, wav_c = pylops.utils.wavelets.ricker(t[:41], f0=par["f0"])
+
+# Generate data
+_, d = pylops.utils.seismicevents.linear2d(x, t, 1500.0, t0, theta, amp, wav)
+
+
+###############################################################################
+# We can now define our operators and apply the forward, adjoint and inverse
+# steps.
+nfft = int(2 ** np.ceil(np.log2(par["nt"])))
+npx, pxmax = 2 * par["nx"], 5e-4
+px = np.linspace(-pxmax, pxmax, npx)
+
+R2Op = pylops.signalprocessing.FourierRadon2D(
+    t, x, px, nfft, kind="linear", engine="numpy", dtype="float64"
+)
+dL_chirp = R2Op.H * d
+dadj_chirp = R2Op * dL_chirp
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
+axs[0].imshow(d.T, vmin=-1, vmax=1, cmap="bwr_r", extent=(x[0], x[-1], t[-1], t[0]))
+axs[0].set(xlabel=r"$x$ [m]", ylabel=r"$t$ [s]", title="Input linear")
+axs[0].axis("tight")
+axs[1].imshow(
+    dL_chirp.T,
+    cmap="bwr_r",
+    vmin=-dL_chirp.max(),
+    vmax=dL_chirp.max(),
+    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),
+)
+axs[1].scatter(1e3 * np.sin(np.deg2rad(theta[0])) / 1500.0, t0[0], s=50, color="r")
+axs[1].set(xlabel=r"$p$ [s/km]", title="Radon")
+axs[1].axis("tight")
+axs[2].imshow(
+    dadj_chirp.T,
+    cmap="bwr_r",
+    vmin=-dadj_chirp.max(),
+    vmax=dadj_chirp.max(),
+    extent=(x[0], x[-1], t[-1], t[0]),
+)
+axs[2].set(xlabel=r"$x$ [m]", title="Adj Radon")
+axs[2].axis("tight")
+plt.tight_layout()
+
+
+###############################################################################
+# We repeat now the same with a parabolic event
+
+# Generate data
+pxx = [1e-6]
+_, d = pylops.utils.seismicevents.parabolic2d(x, t, t0, 0, np.array(pxx), amp, wav)
+
+# Radon transform
+npx, pxmax = 2 * par["nx"], 5e-6
+px = np.linspace(-pxmax, pxmax, npx)
+
+R2Op = pylops.signalprocessing.FourierRadon2D(
+    t, x, px, nfft, kind="parabolic", engine="numpy", dtype="float64"
+)
+dL_chirp = R2Op.H * d
+dadj_chirp = R2Op * dL_chirp
+
+fig, axs = plt.subplots(1, 3, figsize=(12, 4), sharey=True)
+axs[0].imshow(d.T, vmin=-1, vmax=1, cmap="bwr_r", extent=(x[0], x[-1], t[-1], t[0]))
+axs[0].set(xlabel=r"$x$ [m]", ylabel=r"$t$ [s]", title="Input parabolic")
+axs[0].axis("tight")
+axs[1].imshow(
+    dL_chirp.T,
+    cmap="bwr_r",
+    vmin=-dL_chirp.max(),
+    vmax=dL_chirp.max(),
+    extent=(1e3 * px[0], 1e3 * px[-1], t[-1], t[0]),
+)
+axs[1].scatter(1e3 * pxx[0], t0[0], s=50, color="r")
+axs[1].set(xlabel=r"$p$ [s/km]", title="Radon")
+axs[1].axis("tight")
+axs[2].imshow(
+    dadj_chirp.T,
+    cmap="bwr_r",
+    vmin=-dadj_chirp.max(),
+    vmax=dadj_chirp.max(),
+    extent=(x[0], x[-1], t[-1], t[0]),
+)
+axs[2].set(xlabel=r"$x$ [m]", title="Adj Radon")
+axs[2].axis("tight")
+plt.tight_layout()
diff --git a/pylops/signalprocessing/__init__.py b/pylops/signalprocessing/__init__.py
index a8e5ed65..65c0c22e 100755
--- a/pylops/signalprocessing/__init__.py
+++ b/pylops/signalprocessing/__init__.py
@@ -28,6 +28,9 @@
     DTCWT                           Dual-Tree Complex Wavelet Transform.
     Radon2D	                        Two dimensional Radon transform.
     Radon3D	                        Three dimensional Radon transform.
+    FourierRadon2D	                Two dimensional Fourier Radon transform.
+    ChirpRadon2D	                Two dimensional Chirp Radon transform.
+    ChirpRadon3D	                Three dimensional Chirp Radon transform.
     Seislet                         Two dimensional Seislet operator.
     Sliding1D                       1D Sliding transform operator.
     Sliding2D                       2D Sliding transform operator.
@@ -52,6 +55,7 @@
 from .bilinear import *
 from .radon2d import *
 from .radon3d import *
+from .fourierradon2d import *
 from .chirpradon2d import *
 from .chirpradon3d import *
 from .sliding1d import *
diff --git a/pylops/signalprocessing/_fourierradon2d_cuda.py b/pylops/signalprocessing/_fourierradon2d_cuda.py
new file mode 100755
index 00000000..56c14963
--- /dev/null
+++ b/pylops/signalprocessing/_fourierradon2d_cuda.py
@@ -0,0 +1,86 @@
+from cmath import exp
+from math import pi
+
+from numba import cuda
+
+
+@cuda.jit()
+def _radon_inner_2d_kernel(x, y, f, px, h, flim0, flim1, npx, nh):
+    """Cuda kernels for FourierRadon2D operator
+
+    Cuda implementation of the on-the-fly kernel creation and application for the
+    FourierRadon2D operator. See :class:`pylops.signalprocessing.FourierRadon2D`
+    for details about input parameters.
+
+    """
+    ih, ifr = cuda.grid(2)
+    if ih < nh and ifr >= flim0 and ifr <= flim1:
+        for ipx in range(npx):
+            y[ih, ifr] += x[ipx, ifr] * exp(-1j * 2 * pi * f[ifr] * px[ipx] * h[ih])
+
+
+@cuda.jit()
+def _aradon_inner_2d_kernel(x, y, f, px, h, flim0, flim1, npx, nh):
+    """Cuda kernels for FourierRadon2D operator
+
+    Cuda implementation of the on-the-fly kernel creation and application for the
+    FourierRadon2D operator. See :class:`pylops.signalprocessing.FourierRadon2D`
+    for details about input parameters.
+
+    """
+    ipx, ifr = cuda.grid(2)
+    if ipx < npx and ifr >= flim0 and ifr <= flim1:
+        for ih in range(nh):
+            x[ipx, ifr] += y[ih, ifr] * exp(1j * 2 * pi * f[ifr] * px[ipx] * h[ih])
+
+
+def _radon_inner_2d_cuda(
+    x,
+    y,
+    f,
+    px,
+    h,
+    flim0,
+    flim1,
+    npx,
+    nh,
+    num_blocks=(32, 32),
+    num_threads_per_blocks=(32, 32),
+):
+    """Caller for FourierRadon2D operator
+
+    Caller for cuda implementation of matvec kernel for FourierRadon2D operator.
+    See :class:`pylops.signalprocessing.FourierRadon2D` for details about
+    input parameters.
+
+    """
+    _radon_inner_2d_kernel[num_blocks, num_threads_per_blocks](
+        x, y, f, px, h, flim0, flim1, npx, nh
+    )
+    return y
+
+
+def _aradon_inner_2d_cuda(
+    x,
+    y,
+    f,
+    px,
+    h,
+    flim0,
+    flim1,
+    npx,
+    nh,
+    num_blocks=(32, 32),
+    num_threads_per_blocks=(32, 32),
+):
+    """Caller for FourierRadon2D operator
+
+    Caller for cuda implementation of rmatvec kernel for FourierRadon2D operator.
+    See :class:`pylops.signalprocessing.FourierRadon2D` for details about
+    input parameters.
+
+    """
+    _aradon_inner_2d_kernel[num_blocks, num_threads_per_blocks](
+        x, y, f, px, h, flim0, flim1, npx, nh
+    )
+    return x
diff --git a/pylops/signalprocessing/_fourierradon2d_numba.py b/pylops/signalprocessing/_fourierradon2d_numba.py
new file mode 100755
index 00000000..1853c741
--- /dev/null
+++ b/pylops/signalprocessing/_fourierradon2d_numba.py
@@ -0,0 +1,30 @@
+import os
+
+import numpy as np
+from numba import jit, prange
+
+# detect whether to use parallel or not
+numba_threads = int(os.getenv("NUMBA_NUM_THREADS", "1"))
+parallel = True if numba_threads != 1 else False
+
+
+@jit(nopython=True, parallel=parallel, nogil=True, cache=True, fastmath=True)
+def _radon_inner_2d(X, Y, f, px, h, flim0, flim1, npx, nh):
+    for ih in prange(nh):
+        for ifr in range(flim0, flim1):
+            for ipx in range(npx):
+                Y[ih, ifr] += X[ipx, ifr] * np.exp(
+                    -1j * 2 * np.pi * f[ifr] * px[ipx] * h[ih]
+                )
+    return Y
+
+
+@jit(nopython=True, parallel=parallel, nogil=True, cache=True, fastmath=True)
+def _aradon_inner_2d(X, Y, f, px, h, flim0, flim1, npx, nh):
+    for ipx in prange(npx):
+        for ifr in range(flim0, flim1):
+            for ih in range(nh):
+                X[ipx, ifr] += Y[ih, ifr] * np.exp(
+                    1j * 2 * np.pi * f[ifr] * px[ipx] * h[ih]
+                )
+    return X
diff --git a/pylops/signalprocessing/fourierradon2d.py b/pylops/signalprocessing/fourierradon2d.py
new file mode 100755
index 00000000..7a883321
--- /dev/null
+++ b/pylops/signalprocessing/fourierradon2d.py
@@ -0,0 +1,293 @@
+__all__ = ["FourierRadon2D"]
+
+import logging
+from typing import Optional, Tuple
+
+import numpy as np
+import scipy as sp
+
+from pylops import LinearOperator
+from pylops.utils import deps
+from pylops.utils.backend import get_array_module, get_complex_dtype
+from pylops.utils.decorators import reshaped
+from pylops.utils.typing import DTypeLike, NDArray
+
+jit_message = deps.numba_import("the radon2d module")
+
+if jit_message is None:
+
+    from ._fourierradon2d_cuda import _aradon_inner_2d_cuda, _radon_inner_2d_cuda
+    from ._fourierradon2d_numba import _aradon_inner_2d, _radon_inner_2d
+
+logging.basicConfig(format="%(levelname)s: %(message)s", level=logging.WARNING)
+
+
+class FourierRadon2D(LinearOperator):
+    r"""2D Fourier Radon transform
+
+    Apply Radon forward (and adjoint) transform using Fast
+    Fourier Transform to a 2-dimensional array of size
+    :math:`[n_x \times n_t]` (both in forward and adjoint mode).
+
+    Note that forward and adjoint follow the same convention of the time-space
+    implementation in :class:`pylops.signalprocessing.Radon2D`.
+
+    Parameters
+    ----------
+    taxis : :obj:`np.ndarray`
+        Time axis
+    haxis : :obj:`np.ndarray`
+        Spatial axis
+    pxaxis : :obj:`np.ndarray`
+        Axis of scanning variable :math:`p_x` of parametric curve
+    nfft : :obj:`int`
+        Number of samples in Fourier transform
+    flims : :obj:`tuple`, optional
+        Indices of lower and upper limits of Fourier axis to be used in
+        the application of the Radon matrix (if ``None``, use entire axis)
+    kind : :obj:`str`, optional
+        Curve to be used for stacking/spreading (``linear``, ``parabolic``)
+    engine : :obj:`str`, optional
+        Engine used for computation (``numpy`` or ``numba`` or ``cuda``)
+    num_threads_per_blocks : :obj:`tuple`, optional
+        Number of threads in each block (only when ``engine=cuda``)
+    dtype : :obj:`str`, optional
+        Type of elements in input array.
+    name : :obj:`str`, optional
+        Name of operator (to be used by :func:`pylops.utils.describe.describe`)
+
+    Attributes
+    ----------
+    shape : :obj:`tuple`
+        Operator shape
+    explicit : :obj:`bool`
+        Operator contains a matrix that can be solved explicitly (``True``) or
+        not (``False``)
+
+    Notes
+    -----
+    The FourierRadon2D operator applies the Radon transform in the frequency domain.
+    After transforming a 2-dimensional array of size
+    :math:`[n_x \times n_t]` into the frequency domain, the following linear
+    transformation is applied to each frequency component in adjoint mode:
+
+    .. math::
+        \begin{bmatrix}
+            \mathbf{m}(p_{x,1}, \omega_i)  \\
+            \mathbf{m}(p_{x,2}, \omega_i)  \\
+            \vdots          \\
+            \mathbf{m}(p_{x,N_p}, \omega_i)
+        \end{bmatrix}
+        =
+        \begin{bmatrix}
+            e^{-j \omega_i p_{x,1} x^l_1}  & e^{-j \omega_i p_{x,1} x^l_2} &  \ldots & e^{-j \omega_i p_{x,2} x^l_{N_x}}  \\
+            e^{-j \omega_i p_{x,2} x^l_1}  & e^{-j \omega_i p_{x,2} x^l_2} &  \ldots & e^{-j \omega_i p_{x,2} x^l_{N_x}}  \\
+            \vdots            & \vdots           &  \ddots & \vdots            \\
+            e^{-j \omega_i p_{x,N_p} x^l_1}  & e^{-j \omega_i p_{x,N_p} x^l_2} &  \ldots & e^{-j \omega_i p_{x,N_p} x^l_{N_x}}  \\
+        \end{bmatrix}
+        \begin{bmatrix}
+            \mathbf{d}(x_1, \omega_i)  \\
+            \mathbf{d}(x_2, \omega_i)  \\
+            \vdots          \\
+            \mathbf{d}(x_{N_x}, \omega_i)
+        \end{bmatrix}
+
+    where :math:`l=1,2`. Similarly the forward mode is implemented by applying the
+    transpose and complex conjugate of the above matrix to the model transformed to
+    the Fourier domain.
+
+    Refer to [1]_ for more theoretical and implementation details.
+
+    .. [1] Sacchi, M. "Statistical and Transform Methods for
+        Geophysical Signal Processing", 2007.
+
+    """
+
+    def __init__(
+        self,
+        taxis: NDArray,
+        haxis: NDArray,
+        pxaxis: NDArray,
+        nfft: int,
+        flims: Optional[Tuple[int, int]] = None,
+        kind: Optional[str] = "linear",
+        engine: Optional[str] = "numpy",
+        num_threads_per_blocks: Tuple[int, int] = (32, 32),
+        dtype: Optional[DTypeLike] = "float64",
+        name: Optional[str] = "C",
+    ) -> None:
+        # engine
+        if engine not in ["numpy", "numba", "cuda"]:
+            raise KeyError("engine must be numpy or numba or cuda")
+        if engine == "numba" and jit_message is not None:
+            engine = "numpy"
+
+        # dimensions and super
+        dims = len(pxaxis), len(taxis)
+        dimsd = len(haxis), len(taxis)
+        super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name)
+
+        # other input params
+        self.taxis, self.haxis = taxis, haxis
+        self.nh, self.nt = self.dimsd
+        self.px = pxaxis
+        self.npx, self.nfft = self.dims[0], nfft
+        self.dt = taxis[1] - taxis[0]
+        self.dh = haxis[1] - haxis[0]
+        self.f = np.fft.rfftfreq(self.nfft, d=self.dt)
+        self.nfft2 = self.f.size
+        self.cdtype = get_complex_dtype(dtype)
+
+        self.flims = flims
+        if flims is None:
+            self.flims = (0, self.nfft2)
+
+        if kind == "parabolic":
+            self.haxis = self.haxis**2
+
+        # create additional input parameters for engine=cuda
+        if engine == "cuda":
+            self.num_threads_per_blocks = num_threads_per_blocks
+            (
+                num_threads_per_blocks_hpx,
+                num_threads_per_blocks_f,
+            ) = num_threads_per_blocks
+            num_blocks_px = (
+                self.dims[0] + num_threads_per_blocks_hpx - 1
+            ) // num_threads_per_blocks_hpx
+            num_blocks_h = (
+                self.dimsd[0] + num_threads_per_blocks_hpx - 1
+            ) // num_threads_per_blocks_hpx
+            num_blocks_f = (
+                self.dims[1] + num_threads_per_blocks_f - 1
+            ) // num_threads_per_blocks_f
+            self.num_blocks_matvec = (num_blocks_h, num_blocks_f)
+            self.num_blocks_rmatvec = (num_blocks_px, num_blocks_f)
+        self._register_multiplications(engine)
+
+    def _register_multiplications(self, engine: str) -> None:
+        if engine == "numba" and jit_message is None:
+            self._matvec = self._matvec_numba
+            self._rmatvec = self._rmatvec_numba
+        elif engine == "cuda":
+            self._matvec = self._matvec_cuda
+            self._rmatvec = self._rmatvec_cuda
+        else:
+            self._matvec = self._matvec_numpy
+            self._rmatvec = self._rmatvec_numpy
+
+    @reshaped
+    def _matvec_numpy(self, x: NDArray) -> NDArray:
+        ncp = get_array_module(x)
+        x = ncp.fft.rfft(x, n=self.nfft, axis=-1)
+
+        H, PX, F = ncp.meshgrid(
+            self.haxis, self.px, self.f[self.flims[0] : self.flims[1]], indexing="ij"
+        )
+        y = ncp.zeros((self.nh, self.nfft2), dtype=self.cdtype)
+        y[:, self.flims[0] : self.flims[1]] = ncp.einsum(
+            "ijk,jk->ik",
+            ncp.exp(-1j * 2 * ncp.pi * F * PX * H),
+            x[:, self.flims[0] : self.flims[1]],
+        )
+        y = ncp.real(ncp.fft.irfft(y, n=self.nfft, axis=-1))[:, : self.nt]
+        return y
+
+    @reshaped
+    def _rmatvec_numpy(self, y: NDArray) -> NDArray:
+        ncp = get_array_module(y)
+        y = ncp.fft.rfft(y, n=self.nfft, axis=-1)
+
+        PX, H, F = ncp.meshgrid(
+            self.px, self.haxis, self.f[self.flims[0] : self.flims[1]], indexing="ij"
+        )
+        x = ncp.zeros((self.npx, self.nfft2), dtype=self.cdtype)
+        x[:, self.flims[0] : self.flims[1]] = ncp.einsum(
+            "ijk,jk->ik",
+            ncp.exp(1j * 2 * ncp.pi * F * PX * H),
+            y[:, self.flims[0] : self.flims[1]],
+        )
+        x = ncp.real(ncp.fft.irfft(x, n=self.nfft, axis=-1))[:, : self.nt]
+        return x
+
+    @reshaped
+    def _matvec_cuda(self, x: NDArray) -> NDArray:
+        ncp = get_array_module(x)
+        y = ncp.zeros((self.nh, self.nfft2), dtype=self.cdtype)
+
+        x = ncp.fft.rfft(x, n=self.nfft, axis=-1)
+        y = _radon_inner_2d_cuda(
+            x,
+            y,
+            ncp.asarray(self.f),
+            self.px,
+            self.haxis,
+            self.flims[0],
+            self.flims[1],
+            self.npx,
+            self.nh,
+            num_blocks=self.num_blocks_matvec,
+            num_threads_per_blocks=self.num_threads_per_blocks,
+        )
+        y = ncp.real(ncp.fft.irfft(y, n=self.nfft, axis=-1))[:, : self.nt]
+        return y
+
+    @reshaped
+    def _rmatvec_cuda(self, y: NDArray) -> NDArray:
+        ncp = get_array_module(y)
+        x = ncp.zeros((self.npx, self.nfft2), dtype=self.cdtype)
+
+        y = ncp.fft.rfft(y, n=self.nfft, axis=-1)
+        x = _aradon_inner_2d_cuda(
+            x,
+            y,
+            ncp.asarray(self.f),
+            self.px,
+            self.haxis,
+            self.flims[0],
+            self.flims[1],
+            self.npx,
+            self.nh,
+            num_blocks=self.num_blocks_rmatvec,
+            num_threads_per_blocks=self.num_threads_per_blocks,
+        )
+        x = ncp.real(ncp.fft.irfft(x, n=self.nfft, axis=-1))[:, : self.nt]
+        return x
+
+    @reshaped
+    def _matvec_numba(self, x: NDArray) -> NDArray:
+        y = np.zeros((self.nh, self.nfft2), dtype=self.cdtype)
+
+        x = sp.fft.rfft(x, n=self.nfft, axis=-1)
+        y = _radon_inner_2d(
+            x,
+            y,
+            self.f,
+            self.px,
+            self.haxis,
+            self.flims[0],
+            self.flims[1],
+            self.npx,
+            self.nh,
+        )
+        y = np.real(sp.fft.irfft(y, n=self.nfft, axis=-1))[:, : self.nt]
+        return y
+
+    @reshaped
+    def _rmatvec_numba(self, y: NDArray) -> NDArray:
+        x = np.zeros((self.npx, self.nfft2), dtype=self.cdtype)
+
+        y = sp.fft.rfft(y, n=self.nfft, axis=-1)
+        x = _aradon_inner_2d(
+            x,
+            y,
+            self.f,
+            self.px,
+            self.haxis,
+            self.flims[0],
+            self.flims[1],
+            self.npx,
+            self.nh,
+        )
+        x = np.real(sp.fft.irfft(x, n=self.nfft, axis=-1))[:, : self.nt]
+        return x
diff --git a/pylops/signalprocessing/radon2d.py b/pylops/signalprocessing/radon2d.py
index d045eb8a..a76696b4 100644
--- a/pylops/signalprocessing/radon2d.py
+++ b/pylops/signalprocessing/radon2d.py
@@ -218,7 +218,7 @@ def Radon2D(
     -----
     The Radon2D operator applies the following linear transform in adjoint mode
     to the data after reshaping it into a 2-dimensional array of
-    size :math:`[n_x \times n_t]` in adjoint mode:
+    size :math:`[n_x \times n_t]`:
 
     .. math::
         m(p_x, t_0) = \int{d(x, t = f(p_x, x, t))} \,\mathrm{d}x
diff --git a/pylops/utils/estimators.py b/pylops/utils/estimators.py
index 887fac3f..5e73f686 100644
--- a/pylops/utils/estimators.py
+++ b/pylops/utils/estimators.py
@@ -87,7 +87,7 @@ def trace_hutchinson(
         Operator trace.
 
     Raises
-    -------
+    ------
     ValueError
         If ``neval`` is smaller than 3.
 
@@ -194,7 +194,7 @@ def trace_hutchpp(
         Operator trace.
 
     Raises
-    -------
+    ------
     ValueError
         If ``neval`` is smaller than 3.
 
@@ -294,7 +294,7 @@ def trace_nahutchpp(
         Operator trace.
 
     Raises
-    -------
+    ------
     ValueError
         If ``neval`` not large enough to accomodate ``c1`` and ``c2``.
 
diff --git a/pytests/test_fourierradon.py b/pytests/test_fourierradon.py
new file mode 100755
index 00000000..36f3f7b6
--- /dev/null
+++ b/pytests/test_fourierradon.py
@@ -0,0 +1,152 @@
+import multiprocessing
+
+import numpy as np
+import pytest
+from numpy.testing import assert_array_almost_equal
+
+from pylops.optimization.sparsity import fista
+from pylops.signalprocessing import FourierRadon2D
+from pylops.utils import dottest
+
+par1 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 1e-2,
+    "pxmax": 2e-2,
+    "centeredh": True,
+    "kind": "linear",
+    "interp": True,
+    "engine": "numpy",
+}  # linear, centered, linear interp, numpy
+par2 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 1e-2,
+    "pxmax": 2e-2,
+    "centeredh": False,
+    "kind": "linear",
+    "interp": True,
+    "engine": "numpy",
+}  # linear, uncentered, linear interp, numpy
+par3 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 1e-2,
+    "pxmax": 2e-2,
+    "centeredh": True,
+    "kind": "linear",
+    "interp": True,
+    "engine": "numba",
+}  # linear, centered, linear interp, numba
+par4 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 1e-2,
+    "pxmax": 2e-2,
+    "centeredh": False,
+    "kind": "linear",
+    "interp": False,
+    "engine": "numba",
+}  # linear, uncentered, linear interp, numba
+par5 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 8e-3,
+    "pxmax": 7e-3,
+    "centeredh": True,
+    "kind": "parabolic",
+    "interp": False,
+    "engine": "numpy",
+}  # parabolic, centered, no interp, numpy
+par6 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 8e-3,
+    "pxmax": 7e-3,
+    "centeredh": False,
+    "kind": "parabolic",
+    "interp": True,
+    "engine": "numba",
+}  # parabolic, uncentered, interp, numba
+par7 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 9e-2,
+    "pxmax": 8e-2,
+    "centeredh": True,
+    "kind": "hyperbolic",
+    "interp": True,
+    "engine": "numpy",
+}  # hyperbolic, centered, interp, numpy
+par8 = {
+    "nt": 11,
+    "nhx": 21,
+    "nhy": 10,
+    "npx": 21,
+    "npy": 17,
+    "pymax": 7e-2,
+    "pxmax": 8e-2,
+    "centeredh": False,
+    "kind": "hyperbolic",
+    "interp": False,
+    "engine": "numba",
+}  # hyperbolic, uncentered, interp, numba
+
+
+def test_unknown_engine():
+    """Check error is raised if unknown engine is passed"""
+    with pytest.raises(KeyError):
+        _ = FourierRadon2D(None, None, None, None, engine="foo")
+
+
+@pytest.mark.parametrize(
+    "par", [(par1), (par2), (par3), (par4), (par5), (par6), (par7), (par8)]
+)
+def test_FourierRadon2D(par):
+    """Dot-test and sparse inverse for FourierRadon2D operator"""
+    dt, dh = 0.005, 1
+    t = np.arange(par["nt"]) * dt
+    h = np.arange(par["nhx"]) * dh
+    px = np.linspace(0, par["pxmax"], par["npx"])
+    nfft = int(2 ** np.ceil(np.log2(par["nt"])))
+
+    x = np.zeros((par["npx"], par["nt"]))
+    x[2, par["nt"] // 2] = 1
+
+    Rop = FourierRadon2D(
+        t,
+        h,
+        px,
+        nfft,
+        kind=par["kind"],
+        engine=par["engine"],
+        dtype="float64",
+    )
+    assert dottest(Rop, par["nhx"] * par["nt"], par["npx"] * par["nt"], rtol=1e-3)
+
+    y = Rop * x.ravel()
+
+    if par["engine"] == "numba":  # as numpy is too slow here...
+        xinv, _, _ = fista(Rop, y, niter=200, eps=3e0)
+        assert_array_almost_equal(x.ravel(), xinv, decimal=1)
diff --git a/pytests/test_radon.py b/pytests/test_radon.py
index cc4b8054..26d08ba1 100755
--- a/pytests/test_radon.py
+++ b/pytests/test_radon.py
@@ -1,5 +1,3 @@
-import multiprocessing
-
 import numpy as np
 import pytest
 from numpy.testing import assert_array_almost_equal
@@ -176,63 +174,59 @@ def test_Radon2D(par):
     "par", [(par1), (par2), (par3), (par4), (par5), (par6), (par7), (par8)]
 )
 def test_Radon3D(par):
-    """Dot-test,  forward and adjoint consistency check
+    """Dot-test, forward and adjoint consistency check
     (for onthefly parameter), and sparse inverse for Radon3D operator
     """
-    if (
-        par["engine"] == "numpy" or multiprocessing.cpu_count() >= 4
-    ):  # avoid timeout in travis for numba
-
-        dt, dhy, dhx = 0.005, 1, 1
-        t = np.arange(par["nt"]) * dt
-        hy = np.arange(par["nhy"]) * dhy
-        hx = np.arange(par["nhx"]) * dhx
-        py = np.linspace(0, par["pymax"], par["npy"])
-        px = np.linspace(0, par["pxmax"], par["npx"])
-        x = np.zeros((par["npy"], par["npx"], par["nt"]))
-        x[3, 2, par["nt"] // 2] = 1
-
-        Rop = Radon3D(
-            t,
-            hy,
-            hx,
-            py,
-            px,
-            centeredh=par["centeredh"],
-            interp=par["interp"],
-            kind=par["kind"],
-            onthefly=False,
-            engine=par["engine"],
-            dtype="float64",
-        )
-        R1op = Radon3D(
-            t,
-            hy,
-            hx,
-            py,
-            px,
-            centeredh=par["centeredh"],
-            interp=par["interp"],
-            kind=par["kind"],
-            onthefly=True,
-            engine=par["engine"],
-            dtype="float64",
-        )
-
-        assert dottest(
-            Rop,
-            par["nhy"] * par["nhx"] * par["nt"],
-            par["npy"] * par["npx"] * par["nt"],
-            rtol=1e-3,
-        )
-        y = Rop * x.ravel()
-        y1 = R1op * x.ravel()
-        assert_array_almost_equal(y, y1, decimal=4)
-
-        xadj = Rop.H * y
-        xadj1 = R1op.H * y
-        assert_array_almost_equal(xadj, xadj1, decimal=4)
-
-        if Rop.engine == "numba":  # as numpy is too slow here...
-            xinv, _, _ = fista(Rop, y, niter=200, eps=3e0)
-            assert_array_almost_equal(x.ravel(), xinv, decimal=1)
+    dt, dhy, dhx = 0.005, 1, 1
+    t = np.arange(par["nt"]) * dt
+    hy = np.arange(par["nhy"]) * dhy
+    hx = np.arange(par["nhx"]) * dhx
+    py = np.linspace(0, par["pymax"], par["npy"])
+    px = np.linspace(0, par["pxmax"], par["npx"])
+    x = np.zeros((par["npy"], par["npx"], par["nt"]))
+    x[3, 2, par["nt"] // 2] = 1
+
+    Rop = Radon3D(
+        t,
+        hy,
+        hx,
+        py,
+        px,
+        centeredh=par["centeredh"],
+        interp=par["interp"],
+        kind=par["kind"],
+        onthefly=False,
+        engine=par["engine"],
+        dtype="float64",
+    )
+    R1op = Radon3D(
+        t,
+        hy,
+        hx,
+        py,
+        px,
+        centeredh=par["centeredh"],
+        interp=par["interp"],
+        kind=par["kind"],
+        onthefly=True,
+        engine=par["engine"],
+        dtype="float64",
+    )
+
+    assert dottest(
+        Rop,
+        par["nhy"] * par["nhx"] * par["nt"],
+        par["npy"] * par["npx"] * par["nt"],
+        rtol=1e-3,
+    )
+    y = Rop * x.ravel()
+    y1 = R1op * x.ravel()
+    assert_array_almost_equal(y, y1, decimal=4)
+
+    xadj = Rop.H * y
+    xadj1 = R1op.H * y
+    assert_array_almost_equal(xadj, xadj1, decimal=4)
+
+    if Rop.engine == "numba":  # as numpy is too slow here...
+        xinv, _, _ = fista(Rop, y, niter=200, eps=3e0)
+        assert_array_almost_equal(x.ravel(), xinv, decimal=1)