Repository: OverLordGoldDragon/ssqueezepy

Issue:
Batched `icwt`
### Description

Hello, 
I am applying `cwt` and its inverse function `icwt()` to a 2D signal with shape of (channel, time).
`cwt` works without problem, however, `icwt()` gives a `ValueError` as following. I suppose, the inverse function does not support 2D signal yet? It would really nice to have it though :)
![bug_report](https://github.com/user-attachments/assets/b2c99df3-0d4c-48a0-b24d-f9903717aa7c)




### Code for reproduction

```python
import numpy as np
from ssqueezepy import cwt, icwt, Wavelet

fs=2048
duration = 5 #s
nv = 32 
cwt_scales = 'log:maximal'
wavelet = Wavelet('gmw') 

n_time = np.arange(0, duration, 1/fs)
N = len(n_time)
source_1 = np.sin(3*2*np.pi*n_time+2) + 8 #3hz
source_2 = np.sin(2*2*np.pi*n_time+10) + 9 #2hz

ob_1 = (source_1 + source_2) * 1
ob_2 = source_1 * 1/2 + source_2

x = np.array([ob_1,ob_2])
x_mean = x.mean()

Wx, scales = cwt(x, wavelet=wavelet,fs=fs,scales=cwt_scales, nv=nv) 
x_inv = icwt(Wx, wavelet=wavelet,scales=cwt_scales, nv=nv, x_mean=x_mean)
```

Hints from issue comments:
@OverLordGoldDragon: Fair request.

In the meantime, one option is a simple workaround: `Wx.real.sum(axis=1) + x.mean(axis=1)` will be accurate within a constant scaling and shifting factor. Assumes L1 norm, and log `scales` (`log`, `log:minimal`, `log:maximal` are fine).

Notes on the default `scales='log-piecewise'`: 1) its `icwt` is incorrect, I just noticed; 2) my above workaround won't work on it.

Alternatively, a complete match with `icwt`, with mentioned assumptions:

```python
import numpy as np
from ssqueezepy import cwt, icwt
from ssqueezepy.utils import adm_ssq

x = np.random.randn(2, 2048)

nv = 32
wavelet = Wavelet('gmw')
Wx, scales = cwt(x, wavelet, nv=nv, scales='log')

const = (2 / adm_ssq(wavelet)) * np.log(2 ** (1 / nv))
xi = Wx.real.sum(axis=1) * const + x.mean(axis=1, keepdims=True)

# test workaround vs `icwt`
xi0 = icwt(Wx[0], wavelet, scales=scales, nv=nv, x_mean=x[0].mean())
assert np.allclose(xi0, xi[0])

# print MSE
print("MSE = %.3g" % np.mean(np.abs(x - xi)**2))
```
```
MSE = 0.00626
```@OverLordGoldDragon: Was simpler than I thought, implemented. But only `one_int=True` this time.

Solution patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 95773af..04bb83a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,14 @@
+### 0.6.6
+
+#### FEATURES
+ - `icwt` now supports batched `Wx` (3D `Wx`, i.e. `cwt(x)` upon 2D `x`, `(n_inputs, n_times)`)
+
+#### FIXES
+ - `icwt` with `scaletype='linear'`: fix constant scaling factor
+ - scipy deprecation: `scipy.integrate.trapz` -> `scipy.integrate.trapezoid`
+ - numpy deprecation: fix `int()` upon 1D array (of size 1)
+ - numpy deprecation: `np.cfloat` -> `np.complex128`
+
 ### 0.6.5
 
 #### FIXES
diff --git a/ssqueezepy/__init__.py b/ssqueezepy/__init__.py
index ce69f5d..2184efa 100644
--- a/ssqueezepy/__init__.py
+++ b/ssqueezepy/__init__.py
@@ -27,7 +27,7 @@ SOFTWARE.
 """
 
 
-__version__ = '0.6.5'
+__version__ = '0.6.6-dev'
 __title__ = 'ssqueezepy'
 __author__ = 'John Muradeli'
 __license__ = __doc__
diff --git a/ssqueezepy/_cwt.py b/ssqueezepy/_cwt.py
index eed8b05..4339b9d 100644
--- a/ssqueezepy/_cwt.py
+++ b/ssqueezepy/_cwt.py
@@ -327,6 +327,10 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True,
         Wx: np.ndarray
             CWT computed via `ssqueezepy.cwt`.
 
+              - 2D: (n_scales, n_times)
+              - 3D: (n_inputs, n_scales, n_times).
+                Doesn't support `one_int=False`.
+
         wavelet: str / tuple[str, dict] / `wavelets.Wavelet`
             Wavelet sampled in Fourier frequency domain.
                 - str: name of builtin wavelet. `ssqueezepy.wavs()`
@@ -352,6 +356,8 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True,
 
         x_mean: float. mean of original `x` (not picked up in CWT since it's an
             infinite scale component). Default 0.
+            Note: if `Wx` is 3D, `x_mean` should be 1D (`x.mean()` along samples
+            axis).
 
         padtype: str
             Pad scheme to apply on input, in case of `one_int=False`.
@@ -365,7 +371,9 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True,
 
     # Returns:
         x: np.ndarray
-            The signal, as reconstructed from Wx.
+            The signal(s), as reconstructed from Wx.
+
+            If `Wx` is 3D, `x` has shape `(n_inputs, n_times)`.
 
     # References:
         1. One integral inverse CWT. John Muradeli.
@@ -394,7 +402,7 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True,
         synsq_cwt_iw.m
     """
     #### Prepare for inversion ###############################################
-    na, n = Wx.shape
+    *_, na, n = Wx.shape
     x_len = x_len or n
     if not isinstance(scales, np.ndarray) and nv is None:
         nv = 32  # must match forward's; default to `cwt`'s
@@ -414,8 +422,8 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True,
                   padtype=padtype, rpadded=rpadded, l1_norm=l1_norm)
 
         idx = logscale_transition_idx(scales)
-        x  = icwt(Wx[:idx], scales=scales[:idx], **kw)
-        x += icwt(Wx[idx:], scales=scales[idx:], **kw)
+        x  = icwt(Wx[..., :idx, :], scales=scales[:idx], **kw)
+        x += icwt(Wx[..., idx:, :], scales=scales[idx:], **kw)
         return x
     ##########################################################################
 
@@ -423,6 +431,8 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True,
     if one_int:
         x = _icwt_1int(Wx, scales, scaletype, l1_norm)
     else:
+        if Wx.ndim == 3:
+            raise NotImplementedError("batched `Wx` requires `one_int=True`.")
         x = _icwt_2int(Wx, scales, scaletype, l1_norm,
                        wavelet, x_len, padtype, rpadded)
 
@@ -430,11 +440,12 @@ def icwt(Wx, wavelet='gmw', scales='log-piecewise', nv=None, one_int=True,
     Cpsi = (adm_ssq(wavelet) if one_int else
             adm_cwt(wavelet))
     if scaletype == 'log':
-        # Eq 4.67 in [1]; Theorem 4.5 in [1]; below Eq 14 in [2]
+        # Eq 4.67 in [3]; Theorem 4.5 in [3]; below Eq 14 in [5]
         # ln(2**(1/nv)) == ln(2)/nv == diff(ln(scales))[0]
         x *= (2 / Cpsi) * np.log(2 ** (1 / nv))
     else:
-        x *= (2 / Cpsi)
+        # unclear why the `pi/4` here but it improves inversion
+        x *= (2 / Cpsi) * np.pi / 4
 
     x += x_mean  # CWT doesn't capture mean (infinite scale)
     return x
@@ -466,7 +477,7 @@ def _icwt_2int(Wx, scales, scaletype, l1_norm, wavelet, x_len,
 def _icwt_1int(Wx, scales, scaletype, l1_norm):
     """One-integral iCWT; assumes analytic wavelet."""
     norm = _icwt_norm(scaletype, l1_norm)
-    return (Wx.real / norm(scales)).sum(axis=0)
+    return (Wx.real / norm(scales)).sum(axis=-2)
 
 
 def _icwt_norm(scaletype, l1_norm):
diff --git a/ssqueezepy/_ssq_cwt.py b/ssqueezepy/_ssq_cwt.py
index a1bfd14..43cb021 100644
--- a/ssqueezepy/_ssq_cwt.py
+++ b/ssqueezepy/_ssq_cwt.py
@@ -581,7 +581,7 @@ def phase_cwt_num(Wx, dt, difforder=4, gamma=None):
 
     # epsilon from Daubechies, H-T Wu, et al.
     # gamma from Brevdo, H-T Wu, et al.
-    gamma = gamma or 10 * (EPS64 if Wx.dtype == np.cfloat else EPS32)
+    gamma = gamma or 10 * (EPS64 if Wx.dtype == np.complex128 else EPS32)
     w[np.abs(Wx) < gamma] = np.inf
 
     # see `phase_cwt`, though negatives may no longer be in minority
diff --git a/ssqueezepy/ridge_extraction.py b/ssqueezepy/ridge_extraction.py
index 2cc2cb7..2fbbf93 100644
--- a/ssqueezepy/ridge_extraction.py
+++ b/ssqueezepy/ridge_extraction.py
@@ -110,8 +110,8 @@ def extract_ridges(Tf, scales, penalty=2., n_ridges=1, bw=15, transform='cwt',
 
         return ridge_idxs_fw_bw
 
-    eps   = EPS64      if Tf.dtype == np.cfloat else EPS32
-    dtype = np.float64 if Tf.dtype == np.cfloat else np.float32
+    eps   = EPS64      if Tf.dtype == np.complex128 else EPS32
+    dtype = np.float64 if Tf.dtype == np.complex128 else np.float32
     scales, eps, penalty = [np.asarray(x, dtype=dtype)
                             for x in (scales, eps, penalty)]
 
diff --git a/ssqueezepy/utils/cwt_utils.py b/ssqueezepy/utils/cwt_utils.py
index cc0d2dd..7e2dfdf 100644
--- a/ssqueezepy/utils/cwt_utils.py
+++ b/ssqueezepy/utils/cwt_utils.py
@@ -280,7 +280,7 @@ def infer_scaletype(scales):
     if np.mean(np.abs(np.diff(np.log(scales), 2, axis=0))) < th_log:
         scaletype = 'log'
         # ceil to avoid faulty float-int roundoffs
-        nv = int(np.round(1 / np.diff(np.log2(scales), axis=0)[0]))
+        nv = int(np.round(1 / np.diff(np.log2(scales), axis=0)[0].squeeze()))
 
     elif np.mean(np.abs(np.diff(scales, 2, axis=0))) < th_lin:
         scaletype = 'linear'
@@ -620,11 +620,11 @@ def integrate_analytic(int_fn, nowarn=False):
         # (.001 to .1 may not be negligible, however; better captured by logspace)
         t = np.logspace(-15, -1, 1000)
         arr = int_fn(t)
-        return integrate.trapz(arr, t)
+        return integrate.trapezoid(arr, t)
 
     int_nz = _integrate_near_zero()
     arr, t = _find_convergent_array()
-    return integrate.trapz(arr, t) + int_nz
+    return integrate.trapezoid(arr, t) + int_nz
 
 
 def find_max_scale_alt(wavelet, N, min_cutoff=.1, max_cutoff=.8):
diff --git a/ssqueezepy/utils/fft_utils.py b/ssqueezepy/utils/fft_utils.py
index 80cb199..7ae6656 100644
--- a/ssqueezepy/utils/fft_utils.py
+++ b/ssqueezepy/utils/fft_utils.py
@@ -275,7 +275,7 @@ class FFT():
         shape_out = self._get_output_shape(x, axis, real, inverse, n)
 
         # dtypes
-        double = x.dtype in (np.float64, np.cfloat)
+        double = x.dtype in (np.float64, np.complex128)
         cdtype = 'complex128' if double else 'complex64'
         rdtype = 'float64'    if double else 'float32'
         dtype_in  = rdtype if (real and not inverse) else cdtype
diff --git a/ssqueezepy/utils/stft_utils.py b/ssqueezepy/utils/stft_utils.py
index 2d92c93..208d424 100644
--- a/ssqueezepy/utils/stft_utils.py
+++ b/ssqueezepy/utils/stft_utils.py
@@ -208,8 +208,8 @@ def window_resolution(window):
     apsi2   = np.abs(window)**2
     apsih2s = np.abs(psihs)**2
 
-    var_w = integrate.trapz(ws**2 * apsih2s, ws) / integrate.trapz(apsih2s, ws)
-    var_t = integrate.trapz(t**2  * apsi2, t)    / integrate.trapz(apsi2, t)
+    var_w = integrate.trapezoid(ws**2 * apsih2s, ws) / integrate.trapezoid(apsih2s, ws)
+    var_t = integrate.trapezoid(t**2  * apsi2, t)    / integrate.trapezoid(apsi2, t)
 
     std_w, std_t = np.sqrt(var_w), np.sqrt(var_t)
     harea = std_w * std_t
@@ -226,11 +226,11 @@ def window_area(window, time=True, frequency=False):
 
     if time:
         t = np.arange(-len(window)/2, len(window)/2, step=1)
-        at = integrate.trapz(np.abs(window)**2, t)
+        at = integrate.trapezoid(np.abs(window)**2, t)
     if frequency:
         ws = fftshift(_xifn(1, len(window)))
         apsih2s = np.abs(fftshift(fft(window)))**2
-        aw = integrate.trapz(apsih2s, ws)
+        aw = integrate.trapezoid(apsih2s, ws)
 
     if time and frequency:
         return at, aw
diff --git a/ssqueezepy/wavelets.py b/ssqueezepy/wavelets.py
index 8eacf12..cf1e5e1 100644
--- a/ssqueezepy/wavelets.py
+++ b/ssqueezepy/wavelets.py
@@ -703,8 +703,8 @@ def center_frequency(wavelet, scale=None, N=1024, kind='energy', force_int=None,
             scale = (4/pi) * wc_ct
 
         w, psih, apsih2 = _params(wavelet, scale, N)
-        wc = (integrate.trapz(apsih2 * w) /
-              integrate.trapz(apsih2))
+        wc = (integrate.trapezoid(apsih2 * w) /
+              integrate.trapezoid(apsih2))
 
         if use_formula:
             wc *= (scale / scale_orig)
@@ -794,8 +794,8 @@ def freq_resolution(wavelet, scale=10, N=1024, nondim=True, force_int=True,
     wce = center_frequency(wavelet, scale, force_int=force_int, kind='energy')
 
     apsih2 = np.abs(psih)**2
-    var_w = (integrate.trapz((w - wce)**2 * apsih2, w) /
-             integrate.trapz(apsih2, w))
+    var_w = (integrate.trapezoid((w - wce)**2 * apsih2, w) /
+             integrate.trapezoid(apsih2, w))
 
     std_w = np.sqrt(var_w)
     if use_formula:
@@ -910,8 +910,8 @@ def time_resolution(wavelet, scale=10, N=1024, min_decay=1e3, max_mult=2,
     psi = asnumpy(ifft(psih * (-1)**np.arange(Nt)))
 
     apsi2 = np.abs(psi)**2
-    var_t = (integrate.trapz(t**2 * apsi2, t) /
-             integrate.trapz(apsi2, t))
+    var_t = (integrate.trapezoid(t**2 * apsi2, t) /
+             integrate.trapezoid(apsi2, t))
 
     std_t = np.sqrt(var_t)
     if use_formula:


Test patch:
diff --git a/tests/z_all_test.py b/tests/z_all_test.py
index 5d86db4..53bfe79 100644
--- a/tests/z_all_test.py
+++ b/tests/z_all_test.py
@@ -361,7 +361,7 @@ def test_trigdiff():
     dWx = dWx[:, n1:n1+N]
 
     mae = np.mean(np.abs(dWx - dWx2))
-    th = 1e-15 if dWx.dtype == np.cfloat else 1e-7
+    th = 1e-15 if dWx.dtype == np.complex128 else 1e-7
     assert mae < th, mae