diff --git a/pywt/_cwt.py b/pywt/_cwt.py
index 0d69095b..ed69504a 100644
--- a/pywt/_cwt.py
+++ b/pywt/_cwt.py
@@ -24,9 +24,8 @@ def next_fast_len(n):
     return 2**ceil(np.log2(n))
 
 
-def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
+def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1, *, hop_size=1):
     """
-    cwt(data, scales, wavelet)
 
     One dimensional Continuous Wavelet Transform.
 
@@ -60,6 +59,14 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
     axis: int, optional
         Axis over which to compute the CWT. If not given, the last axis is
         used.
+    hop_size : int
+        Specifies the down-sampling factor applied on temporal axis during the transform.
+        The output is sampled every hop size samples, rather than at every consecutive sample.
+        For example:
+        A signal of length 1024 yields 1024 output samples when ``hop_size=1``;
+        512 output samples when ``hop_size=2``;
+        256 output samples when ``hop_size=4``.
+        ``hop_size`` must be a positive integer (≥1).
 
     Returns
     -------
@@ -73,8 +80,14 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
 
     Notes
     -----
-    Size of coefficients arrays depends on the length of the input array and
-    the length of given scales.
+    Size of coefficients arrays depends on the length of the input array,
+    the length of given scales and the given hop size.
+    
+    References
+    ----------
+    .. [1] Phan, D. T., Huynh, T. A., Pham, V. T., Tran, C. M., Mai, V. T., & Tran, N. Q. (2025). 
+           Optimal Scalogram for Computational Complexity Reduction in Acoustic Recognition Using Deep Learning. 
+           *arXiv preprint arXiv:2505.13017*. https://arxiv.org/abs/2505.13017
 
     Examples
     --------
@@ -83,7 +96,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
     >>> import matplotlib.pyplot as plt
     >>> x = np.arange(512)
     >>> y = np.sin(2*np.pi*x/32)
-    >>> coef, freqs=pywt.cwt(y,np.arange(1,129),'gaus1')
+    >>> coef, freqs=pywt.cwt(y,np.arange(1,129),'gaus1',hop_size=1)
     >>> plt.matshow(coef)
     >>> plt.show()
 
@@ -93,7 +106,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
     >>> t = np.linspace(-1, 1, 200, endpoint=False)
     >>> sig  = np.cos(2 * np.pi * 7 * t) + np.real(np.exp(-7*(t-0.4)**2)*np.exp(1j*2*np.pi*2*(t-0.4)))
     >>> widths = np.arange(1, 31)
-    >>> cwtmatr, freqs = pywt.cwt(sig, widths, 'mexh')
+    >>> cwtmatr, freqs = pywt.cwt(sig, widths, 'mexh', hop_size=2)
     >>> plt.imshow(cwtmatr, extent=[-1, 1, 1, 31], cmap='PRGn', aspect='auto',
     ...            vmax=abs(cwtmatr).max(), vmin=-abs(cwtmatr).max())
     >>> plt.show()
@@ -114,7 +127,8 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
         raise AxisError("axis must be a scalar.")
 
     dt_out = dt_cplx if wavelet.complex_cwt else dt
-    out = np.empty((np.size(scales),) + data.shape, dtype=dt_out)
+    data_sampled = data[..., ::hop_size]
+    out = np.empty((np.size(scales),) + data_sampled.shape, dtype=dt_out)
     precision = 10
     int_psi, x = integrate_wavelet(wavelet, precision=precision)
     int_psi = np.conj(int_psi) if wavelet.complex_cwt else int_psi
@@ -187,7 +201,7 @@ def cwt(data, scales, wavelet, sampling_period=1., method='conv', axis=-1):
             # restore original data shape and axis position
             coef = coef.reshape(data_shape_pre)
             coef = coef.swapaxes(axis, -1)
-        out[i, ...] = coef
+        out[i, ...] = coef[..., ::hop_size]
 
     frequencies = scale2frequency(wavelet, scales, precision)
     if np.isscalar(frequencies):
diff --git a/pywt/tests/test_cwt_wavelets.py b/pywt/tests/test_cwt_wavelets.py
index f142404c..0e0d7f7b 100644
--- a/pywt/tests/test_cwt_wavelets.py
+++ b/pywt/tests/test_cwt_wavelets.py
@@ -381,7 +381,7 @@ def test_cwt_complex(dtype, tol, method):
     scales = np.arange(1, 32)
 
     # real-valued tranfsorm as a reference
-    [cfs, f] = pywt.cwt(sst, scales, wavelet, dt, method=method)
+    [cfs, f] = pywt.cwt(sst, scales, wavelet, dt, method=method, hop_size=1)
 
     # verify same precision
     assert_equal(cfs.real.dtype, sst.dtype)
@@ -390,7 +390,7 @@ def test_cwt_complex(dtype, tol, method):
     # and imaginary components
     sst_complex = sst + 1j*sst
     [cfs_complex, f] = pywt.cwt(sst_complex, scales, wavelet, dt,
-                                method=method)
+                                method=method, hop_size=1)
     assert_allclose(cfs + 1j*cfs, cfs_complex, atol=tol, rtol=tol)
     # verify dtype is preserved
     assert_equal(cfs_complex.dtype, sst_complex.dtype)
@@ -409,10 +409,10 @@ def test_cwt_batch(axis, method):
     scales = np.arange(1, 32)
 
     # non-batch transform as reference
-    [cfs1, f] = pywt.cwt(sst1, scales, wavelet, dt, method=method, axis=axis)
+    [cfs1, f] = pywt.cwt(sst1, scales, wavelet, dt, method=method, axis=axis, hop_size=1)
 
     shape_in = sst.shape
-    [cfs, f] = pywt.cwt(sst, scales, wavelet, dt, method=method, axis=axis)
+    [cfs, f] = pywt.cwt(sst, scales, wavelet, dt, method=method, axis=axis, hop_size=1)
 
     # shape of input is not modified
     assert_equal(shape_in, sst.shape)