[FEATURE] add the enhanced spectrogram

pytta/__init__.py

@@ -113,7 +113,7 @@
     list_devices, print_devices, get_device_from_user,\
     fft_convolve, find_delay, corr_coef, resample,\
     peak_time, save, load, plot_time, plot_time_dB, plot_freq,\
-    plot_bars, plot_spectrogram, plot_waterfall, SPL
+    plot_bars, plot_spectrogram, plot_waterfall, SPL, enhanced_spectrogram
 
 # Must go on v0.1.0 (DEPRECATED)
 from .functions import fft_degree
@@ -146,6 +146,7 @@
            'resample',
            'corr_coef',
            'peak_time',
+           'enhanced_spectrogram',
            'save',
            'load',
            'weighting',

pytta/_plot.py

@@ -26,6 +26,7 @@
 import scipy.signal as ss
 import copy as cp
 from plotly import io
+from pytta.utils.spectrogram import enhanced_map
 io.renderers.default = 'browser'
 
 
@@ -804,7 +805,7 @@ def _curve_data_extractor_bars(analyses):
 
 def spectrogram(sigObjs, winType, winSize,
                 overlap, xLabel, yLabel, xLim, yLim,
-                title, decimalSep):
+                title, decimalSep, enhanced=None):
     """
     Plots a signal spectrogram in frequency domain.
 
@@ -846,6 +847,10 @@ def spectrogram(sigObjs, winType, winSize,
             may be dot or comma.
 
             >>> decimalSep = ',' # in Brazil
+        * enhanced (None), (str):
+            'readable' or 'sharp' plots the enhanced spectrogram
+            (pytta.enhanced_spectrogram) on a logarithmic frequency axis;
+            winType, winSize and overlap are then not used.
 
     Return:
     --------
@@ -867,13 +872,23 @@ def spectrogram(sigObjs, winType, winSize,
                         projection='rectilinear')
         ax.set_snap(False)
 
-        _spectrogram, _specTime, _specFreq = \
-            _calc_spectrogram(data['timeSignal'], data['timeVector'],
-                              data['samplingRate'], overlap, winType,
-                              winSize, data['dBRef'])
-
-        pcmesh = ax.pcolormesh(_specTime, _specFreq, _spectrogram,
-                            cmap=plt.jet(), vmin=-120)
+        if enhanced:
+            _specTime, _specFreq, meanSquare = \
+                enhanced_map(data['timeSignal'], data['samplingRate'],
+                             enhanced)
+            _spectrogram = 10 * np.log10(meanSquare / data['dBRef'] ** 2
+                                         + 1e-12)
+            pcmesh = ax.pcolormesh(_specTime, _specFreq, _spectrogram,
+                                   cmap=plt.jet(), shading='nearest',
+                                   vmin=_spectrogram.max() - 80)
+            ax.set_yscale('log')
+        else:
+            _spectrogram, _specTime, _specFreq = \
+                _calc_spectrogram(data['timeSignal'], data['timeVector'],
+                                  data['samplingRate'], overlap, winType,
+                                  winSize, data['dBRef'])
+            pcmesh = ax.pcolormesh(_specTime, _specFreq, _spectrogram,
+                                cmap=plt.jet(), vmin=-120)
 
         if xLim is None:
             xLim = (data['timeVector'][0], data['timeVector'][-1])
@@ -883,8 +898,15 @@ def spectrogram(sigObjs, winType, winSize,
             yLim = (data['minFreq'], data['maxFreq'])
         ax.set_ylim(yLim)
 
-        ax.yaxis.set_major_locator(ticker.MaxNLocator(min_n_ticks=8))
-        ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useOffset=True))
+        if enhanced:
+            octaves = [f for f in (31.5, 63, 125, 250, 500, 1000, 2000, 4000,
+                                   8000, 16000) if yLim[0] <= f <= yLim[1]]
+            ax.set_yticks(octaves)
+            ax.set_yticklabels([str(f) for f in octaves])
+            ax.yaxis.set_minor_formatter(ticker.NullFormatter())
+        else:
+            ax.yaxis.set_major_locator(ticker.MaxNLocator(min_n_ticks=8))
+            ax.yaxis.set_major_formatter(ticker.ScalarFormatter(useOffset=True))
         ax.xaxis.set_major_locator(ticker.MaxNLocator(min_n_ticks=10))
         ax.xaxis.set_major_formatter(ticker.ScalarFormatter(useOffset=True))
         for item in (ax.get_xticklabels() + ax.get_yticklabels()):

pytta/classes/signal.py

@@ -738,7 +738,8 @@ def plot_freq(self, smooth:bool=False, xLabel:str=None, yLabel:str=None,
     def plot_spectrogram(self, winType:str='hann', winSize:int=1024,
                          overlap:float=0.5, xLabel:str=None, yLabel:str=None,
                          yLim:list=None, xLim:list=None, title:str=None,
-                         decimalSep:str=','):
+                         decimalSep:str=',',
+                         enhanced:str=None):
         """Plots a signal spectrogram.
 
         xLabel, yLabel, and title are saved for the next plots when provided.
@@ -778,6 +779,10 @@ def plot_spectrogram(self, winType:str='hann', winSize:int=1024,
                 may be dot or comma.
 
                 >>> decimalSep = ',' # in Brazil
+            * enhanced (None), (str):
+                'readable' or 'sharp' plots the enhanced spectrogram
+                (pytta.enhanced_spectrogram) on a logarithmic frequency
+                axis; winType, winSize and overlap are then not used.
 
         Return:
         --------
@@ -807,7 +812,7 @@ def plot_spectrogram(self, winType:str='hann', winSize:int=1024,
 
         figs = plot.spectrogram((self,), winType, winSize,
                                 overlap, xLabel, yLabel, xLim, yLim,
-                                title, decimalSep)
+                                title, decimalSep, enhanced)
         return figs
 
     def calib_voltage(self, chIndex, refSignalObj, refVrms=1, refFreq=1000):

pytta/functions.py

@@ -17,6 +17,7 @@
     >>> pytta.corr_coef(signalObj1, signalObj2)
     >>> pytta.resample(signalObj, newSamplingRate)
     >>> pytta.peak_time(signalObj1, signalObj2, ..., signalObjN)
+    >>> pytta.enhanced_spectrogram(signalObj)
     >>> pytta.plot_time(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.plot_time_dB(signalObj1, signalObj2, ..., signalObjN)
     >>> pytta.plot_freq(signalObj1, signalObj2, ..., signalObjN)
@@ -49,6 +50,7 @@
 import copy as cp
 from warnings import warn
 from pytta import _plot as plot
+from pytta.utils.spectrogram import enhanced_map
 # For backwards compatibility purposes. Planned to get out of here
 from pytta.utils.maths import fft_degree as new_fft_degree
 
@@ -274,6 +276,62 @@ def peak_time(signal):
     else:
         return peaks_time[0]
 
+def enhanced_spectrogram(signalObj, smoothing='readable', n_frames=2000,
+                         f_min=20, bins_per_octave=192):
+    """
+    Enhanced spectrogram: a spectrogram with no window to choose.
+
+    Nine Blackman-Harris windows from 8 to 512 ms (geometric spacing) are
+    reassigned (Auger and Flandrin, 1995), the energy of each frame is held
+    over the frame step of its window, each map is smoothed, and the maps are
+    combined by geometric mean, so that only the energy that all the windows
+    place at the same point of the time-frequency plane remains (the idea of
+    Cheung and Lim, 1991, applied to reassigned maps).
+
+    Parameters (default), (type):
+    -----------------------------
+        * signalObj (), (SignalObj):
+            signal to analyse. It is not modified;
+        * smoothing ('readable'), (str):
+            'readable': Gaussian smoothing of 4 ms in time and 1.45 % of the
+            frequency, for continuous lines; 'sharp': 1 ms and 1 Hz, for the
+            thinnest lines;
+        * n_frames (2000), (int):
+            largest number of time columns. The time step is 1 ms or more;
+        * f_min (20), (float):
+            lowest frequency of the output, in Hz;
+        * bins_per_octave (192), (int):
+            frequency resolution of the logarithmic output grid.
+
+    Return:
+    -------
+        * t: time of each column, in s, shape (M,);
+        * f: frequency of each row, in Hz, shape (F,);
+        * L: level of each cell in dB re the channel's dBRef (dB SPL for a
+          channel in Pa), shape (F, M, numChannels). A cell holds the mean
+          square value of its share of the signal, so the cells of a column
+          add up to the level of that instant; a steady tone spread by the
+          smoothing over several cells reads lower in each of them (about
+          10 dB with 'readable'). Silence sits at -120 dB.
+
+    Author: Sergio Aguirre, 2026. Released in PyTTa under the MIT license.
+    The same method is in the SQAT toolbox (GPL-3.0) as
+    gui/SQAT_GUI_enhanced_stft.m (https://github.com/ggrecow/SQAT).
+
+    >>> t, f, L = pytta.enhanced_spectrogram(signalObj, smoothing='sharp')
+    """
+    if not isinstance(signalObj, SignalObj):
+        raise TypeError('signalObj must be a SignalObj.')
+    L = []
+    for chIndex, chNum in enumerate(signalObj.channels.mapping):
+        t, f, meanSquare = enhanced_map(signalObj.timeSignal[:, chIndex],
+                                        signalObj.samplingRate, smoothing,
+                                        n_frames, f_min, bins_per_octave)
+        dBRef = signalObj.channels[chNum].dBRef
+        L.append(10 * np.log10(meanSquare / dBRef ** 2 + 1e-12))
+    return t, f, np.stack(L, axis=2)
+
+
 def fft_degree(*args,**kwargs):
     """
     DEPRECATED
@@ -533,7 +591,8 @@ def plot_bars(*analyses, xLabel:str=None, yLabel:str=None,
 def plot_spectrogram(*sigObjs, winType:str='hann', winSize:int=1024,
                      overlap:float=0.5, xLabel:str=None, yLabel:str=None,
                      yLim:list=None, xLim:list=None, title:str=None,
-                     decimalSep:str=','):
+                     decimalSep:str=',',
+                     enhanced:str=None):
     """
     Plots provided SignalObjs spectrogram.
 
@@ -575,6 +634,10 @@ def plot_spectrogram(*sigObjs, winType:str='hann', winSize:int=1024,
             may be dot or comma.
 
             >>> decimalSep = ',' # in Brazil
+        * enhanced (None), (str):
+            'readable' or 'sharp' plots the enhanced spectrogram
+            (pytta.enhanced_spectrogram) on a logarithmic frequency axis;
+            winType, winSize and overlap are then not used.
 
     Return:
     --------
@@ -586,7 +649,7 @@ def plot_spectrogram(*sigObjs, winType:str='hann', winSize:int=1024,
     if len(realSigObjs) > 0:
         figs = plot.spectrogram(realSigObjs, winType, winSize,
                                 overlap, xLabel, yLabel, xLim, yLim,
-                                title, decimalSep)
+                                title, decimalSep, enhanced)
         return figs
     else:
         return

pytta/utils/__init__.py

@@ -15,6 +15,7 @@
     * colore
     * freq
     * maths
+    * spectrogram
 
 
 Created on Tue May 5 00:31:42 2020

pytta/utils/spectrogram.py

@@ -0,0 +1,144 @@
+"""
+Enhanced spectrogram on NumPy arrays.
+
+Author: Sergio Aguirre, 2026. Released in PyTTa under the MIT license.
+The same method is in the SQAT toolbox (GPL-3.0) as
+gui/SQAT_GUI_enhanced_stft.m (https://github.com/ggrecow/SQAT).
+The SignalObj interface is pytta.enhanced_spectrogram.
+"""
+import numpy as np
+from scipy.ndimage import convolve1d
+
+__all__ = ['enhanced_map']
+
+
+def enhanced_map(x, fs, smoothing='readable', n_frames=2000, f_min=20,
+                 bins_per_octave=192):
+    """
+    Enhanced spectrogram of one channel, as a map of mean square values.
+
+    Nine Blackman-Harris windows from 8 to 512 ms (geometric spacing) are
+    reassigned (Auger and Flandrin, 1995), the energy of each frame is held
+    over the frame step of its window, each map is smoothed, and the maps are
+    combined by geometric mean, so that only the energy that all the windows
+    place at the same point of the time-frequency plane remains (the idea of
+    Cheung and Lim, 1991, applied to reassigned maps).
+
+    Returns t (M,) in s, f (F,) in Hz on a logarithmic grid, and the map
+    (F, M) of mean square values in the unit of x squared: the cells of a
+    column add up to the mean square value of that instant. See
+    pytta.enhanced_spectrogram for the parameters.
+    """
+    if f_min >= fs / 2:
+        raise ValueError(f'f_min ({f_min} Hz) must be below fs/2 '
+                         f'({fs / 2} Hz).')
+    if smoothing not in ('readable', 'sharp'):
+        raise ValueError("smoothing must be 'readable' or 'sharp'.")
+    x = np.asarray(x, dtype=float).ravel()
+    numSamples = x.size
+    durations = 8e-3 * 64 ** (np.arange(9) / 8)
+    windows = 2 * _round_half_away(durations * fs / 2).astype(int)
+    hop = max(1, int(_round_half_away(0.001 * fs)),
+              int(np.ceil(numSamples / n_frames)))
+    steps = np.maximum(1, _round_half_away(windows / 8).astype(int))
+    numColumns = len(range(0, numSamples, hop))
+    numRows = int(np.ceil(np.log2((fs / 2) / f_min) * bins_per_octave)) + 1
+    f = f_min * 2 ** (np.arange(numRows) / bins_per_octave)
+    cellWidth = f * (2 ** (1 / bins_per_octave) - 1)
+    if smoothing == 'readable':
+        sigmaTime = 0.004 * fs / hop
+        sigmaFreq = 2 * (2 ** (1 / 96) - 1) * f / cellWidth
+    else:
+        sigmaTime = 0.001 * fs / hop
+        sigmaFreq = 1 / cellWidth
+    timeKernel = _gaussian_kernel(max(sigmaTime, 0.3))
+    sigmaFreq = np.maximum(sigmaFreq, 0.3)
+    logSum = np.zeros((numColumns, numRows))
+    for N, step in zip(windows, steps):
+        energy = np.concatenate(([0.], np.cumsum(np.concatenate(
+            (np.zeros(N), x ** 2, np.zeros(N))))))
+        frameMax = np.max(energy[N:] - energy[:-N])
+        R = _reassigned_map(x, fs, N, step, hop, numColumns, numRows,
+                            f_min, bins_per_octave, frameMax) / hop
+        numHold = 2 * int(np.floor(step / hop / 2)) + 1
+        if numHold > 1:
+            R = convolve1d(R, np.ones(numHold) / numHold, axis=0,
+                           mode='constant')
+        R = convolve1d(R, timeKernel, axis=0, mode='constant')
+        R = _smooth_frequency_rows(R, sigmaFreq)
+        with np.errstate(divide='ignore'):
+            logSum += np.log(R + 1e-6 * R.max())
+    t = np.arange(numColumns) * hop / fs
+    return t, f, np.exp(logSum / windows.size).T
+
+
+def _round_half_away(x):
+    """Round halves away from zero, as MATLAB does."""
+    x = np.asarray(x, dtype=float)
+    return np.sign(x) * np.floor(np.abs(x) + 0.5)
+
+
+def _gaussian_kernel(sigma):
+    half = int(np.ceil(3 * sigma))
+    kernel = np.exp(-np.arange(-half, half + 1) ** 2 / (2 * sigma ** 2))
+    return kernel / kernel.sum()
+
+
+def _smooth_frequency_rows(R, sigma):
+    """Gaussian smoothing along the frequency axis (columns of R) with a
+    width that changes with frequency: one kernel per geometric level."""
+    level = _round_half_away(np.log(sigma) / np.log(1.25)).astype(int)
+    out = R.copy()
+    numFreqs = R.shape[1]
+    for lev in np.unique(level):
+        cols = np.flatnonzero(level == lev)
+        kernel = _gaussian_kernel(1.25 ** lev)
+        half = kernel.size // 2
+        c1 = max(0, cols[0] - half)
+        c2 = min(numFreqs, cols[-1] + half + 1)
+        smoothed = convolve1d(R[:, c1:c2], kernel, axis=1, mode='constant')
+        out[:, cols] = smoothed[:, cols - c1]
+    return out
+
+
+def _reassigned_map(x, fs, N, step, hop, numColumns, numRows, f_min,
+                    bins_per_octave, frameMax):
+    """Reassigned spectrogram of one unit-energy Blackman-Harris window of N
+    samples and frame step `step`, accumulated on the output grid."""
+    half = N // 2
+    padded = np.concatenate((np.zeros(half), x, np.zeros(half + N)))
+    n = np.arange(-half, half)
+    phase = 2 * np.pi * (n + half) / N
+    a = (0.35875, 0.48829, 0.14128, 0.01168)
+    w = a[0] - a[1] * np.cos(phase) + a[2] * np.cos(2 * phase) \
+        - a[3] * np.cos(3 * phase)
+    wd = a[1] * (2 * np.pi / N) * np.sin(phase) \
+        - a[2] * (4 * np.pi / N) * np.sin(2 * phase) \
+        + a[3] * (6 * np.pi / N) * np.sin(3 * phase)
+    norm = 1 / np.linalg.norm(w)
+    w, wd = w * norm, wd * norm
+    wt = n * w
+    freqBins = np.arange(half + 1) * fs / N
+    starts = np.arange(0, x.size, step)
+    R = np.zeros(numColumns * numRows)
+    chunk = max(1, int(2e6 // N))
+    for b in range(0, starts.size, chunk):
+        s = starts[b:b + chunk]
+        frames = padded[s[:, None] + np.arange(N)[None, :]]
+        Xh = np.fft.rfft(frames * w, axis=1)
+        Xt = np.fft.rfft(frames * wt, axis=1)
+        Xd = np.fft.rfft(frames * wd, axis=1)
+        power = np.abs(Xh) ** 2
+        r = np.conj(Xh) / np.maximum(power, np.finfo(float).eps)
+        tHat = s[:, None] + np.real(Xt * r)
+        fHat = freqBins[None, :] - fs / (2 * np.pi) * np.imag(Xd * r)
+        col = _round_half_away(tHat / hop)
+        with np.errstate(divide='ignore'):
+            row = _round_half_away(bins_per_octave * np.log2(
+                np.maximum(fHat, np.finfo(float).eps) / f_min))
+        ok = (power > 1e-10 * frameMax) & (col >= 0) & (col < numColumns) \
+            & (row >= 0) & (row < numRows) & (fHat > 0)
+        cells = (col[ok] + row[ok] * numColumns).astype(int)
+        R += np.bincount(cells, weights=power[ok] * (2 * step / N),
+                         minlength=R.size)
+    return R.reshape(numRows, numColumns).T