any reason psd.spl function does not have a kwargs to getresp?

[cluffspace]

Maybe i'm missing the big pictures here but it seems very straight forward to request the psdmap using getresp in psdmod and apply the conversion to dB to the psdmap and get P by finding the maximum therefore both SPL, SPL spectrogram and OASPL could be available.

[twmacro]

Hi @cluffspace,

I like your thinking here ... a nice efficiency boost if you want all that info. To make this work without breaking backward compatibility, I'm thinking I could add a parameter to psd.spl called "get_psd_info" (for example). If truthy, return the output of psdmod in a namespace or dict. Further, if get_psd_info is specifically the string "getmap", it would turn on that feature of psdmod also. This is after a quick look only, but does that make sense? Or did I misinterpret what you're saying?

Regards,
Tim

[cluffspace]

does my code make sense?
def spl(
x: np.ndarray,
sr: float,
nperseg: Optional[int] = None,
overlap: float = 0.5,
window: str | tuple | np.ndarray = "hann",
timeslice: float | str = 1.0,
tsoverlap: float | str = 0.5,
fs: int = 3,
pref: float = 2.9e-9,
frange: tuple[float, float] = (25.0, np.inf),
extendends: bool = True,
getmap: bool = False
) -> tuple[np.ndarray, np.ndarray, float] | tuple[np.ndarray, np.ndarray, float, np.ndarray, np.ndarray]:
"""
Sound pressure level estimation using PSD with optional spectrogram.

Modified version of pyyeti.psd.spl that supports getmap parameter
to return SPL spectrogram (time-frequency response).

Parameters
----------
x : np.ndarray
    Vector of pressure values (1D).
sr : float
    Sample rate (Hz).
nperseg : int, optional
    Length of each segment for the FFT. Defaults to
    ``int(sr / 5)`` for 5 Hz frequency step in PSD.
overlap : float, optional
    Amount of overlap in windows, e.g., 0.5 would be 50% overlap.
window : str or tuple or array_like, optional
    Passed to :func:`scipy.signal.welch`.
timeslice : float or str
    If scalar, length in seconds for each slice. If string,
    integer number of points for each slice.
tsoverlap : float or str
    If scalar, fraction of each time-slice to overlap [0, 1).
    If string, integer number of points to overlap.
fs : int, optional
    Specifies output frequency scale:
    - 0: linear scale
    - 1: full octave scale
    - 3: 1/3 octave scale
    - 6: 1/6 octave scale
pref : float, optional
    Reference pressure. Default 2.9e-9 psi (20e-6 Pa).
frange : tuple, optional
    (fmin, fmax) frequency bounds. Default (25.0, inf).
extendends : bool, optional
    Passed to :func:`pyyeti.psd.rescale` if octave scale requested.
getmap : bool, optional
    If True, return SPL spectrogram (Smap and t). Default False.

Returns
-------
f : np.ndarray
    Frequency vector (Hz).
spl : np.ndarray
    Maximax SPL vector in dB (1D).
oaspl : float
    Overall sound pressure level (dB).
Smap : np.ndarray (only if getmap=True)
    SPL spectrogram (freq × time). Rows are frequencies,
    columns are time slices.
t : np.ndarray (only if getmap=True)
    Time vector for columns of Smap (seconds).

Notes
-----
This function is a modified version of :func:`pyyeti.psd.spl` that
adds spectrogram output when ``getmap=True``.

Math follows pyyeti.psd.spl:
    V = (PSD * Δf) / P_ref^2
    SPL = 10 * log10(V)
    OASPL = 10 * log10(sum(V))

When getmap=False:
    - Calls pyyeti.psd.psdmod() without getmap
    - Returns (f, spl, oaspl) - matches pyyeti.psd.spl

When getmap=True:
    - Calls pyyeti.psd.psdmod(getmap=True) to get PSD map
    - Converts PSD map to SPL map
    - Returns (f, spl, oaspl, Smap, t)

See Also
--------
pyyeti.psd.spl : Original function without getmap support
pyyeti.psd.psdmod : Underlying PSD calculation with getmap
"""
from pyyeti.psd import psdmod, rescale, _set_frange

if nperseg is None:
    nperseg = int(sr / 5)

# Always compute PSD with map - psdmod calculates it internally regardless
F, _, Pmap, t = psdmod(
    x, sr,
    nperseg=nperseg,
    window=window,
    timeslice=timeslice,
    tsoverlap=tsoverlap,
    noverlap=int(overlap * nperseg),
    getmap=True
)

# Apply frequency range and octave scaling
s, e = _set_frange(frange, 0.0, sr / 2)

if fs != 0:
    # Octave band scaling - rescale operates columnwise on matrices
    Pmap_scaled, F_scaled, _, _ = rescale(Pmap, F, n_oct=fs, frange=(s, e), extendends=extendends)
else:
    # Linear frequency scale
    df = F[1]
    Pmap_ms = Pmap * df  # Convert to mean-square per band

    # Apply frequency range filter
    pv = (F >= s) & (F <= e)
    F_scaled = F[pv]
    Pmap_scaled = Pmap_ms[pv, :]

# Convert PSD map to SPL map and compute maximax/OASPL
v_map = Pmap_scaled / pref**2
Smap = 10 * np.log10(v_map)

# Maximax SPL: take max over time (axis=1) for each frequency
spl = np.max(Smap, axis=1)

# OASPL: sum mean-square values over frequency, then convert to dB
oaspl = 10 * np.log10(np.sum(np.max(v_map, axis=1)))

if not getmap:
    return F_scaled, spl, oaspl
else:
    return F_scaled, spl, oaspl, Smap, t

[twmacro]

Hi @cluffspace,

Oh, I see. I did misunderstand what you were looking for. This is a nice update! I had to make some small mods to your routine for when octave scaling is used however. Here is the final code that passed all my tests (for brevity, I chopped off the top part of the routine ... mostly the documentation):

# Always compute PSD with map - psdmod calculates it internally regardless
F, P, Pmap, t = psdmod(
    x,
    sr,
    nperseg=nperseg,
    window=window,
    timeslice=timeslice,
    tsoverlap=tsoverlap,
    noverlap=int(overlap * nperseg),
    getmap=True,
)

# Apply frequency range and octave scaling
s, e = _set_frange(frange, 0.0, sr / 2)

if fs != 0:
    # Octave band scaling - rescale operates columnwise on matrices
    # - augment map with overall PSD so we only call rescale one time:
    Pmap_augmented = np.column_stack((P, Pmap))

    _, F_scaled, _, Pmap_scaled_augmented = rescale(
        Pmap_augmented, F, n_oct=fs, frange=(s, e), extendends=extendends
    )

    Pmap_scaled = Pmap_scaled_augmented[:, 1:]
    P_scaled = Pmap_scaled_augmented[:, 0]
else:
    # Linear frequency scale
    df = F[1]
    Pmap_ms = Pmap * df  # Convert to mean-square per band

    # Apply frequency range filter
    pv = (F >= s) & (F <= e)
    F_scaled = F[pv]
    Pmap_scaled = Pmap_ms[pv, :]
    P_scaled = P[pv] * df

# Convert PSD map to SPL map and compute maximax/OASPL
v_map = Pmap_scaled / pref**2
Smap = 10 * np.log10(v_map)

# SPL is computed from overall mean-square values per band (P_scaled):
v = P_scaled / pref**2
spl = 10 * np.log10(v)

# OASPL: sum mean-square values over frequency, then convert to dB
oaspl = 10 * np.log10(np.sum(v))

if not getmap:
    return F_scaled, spl, oaspl
else:
    return F_scaled, spl, oaspl, Smap, t

I had to augment the map with the PSD in order for the new routine to exactly match the old routine (when octave scaling is used; for linear scales, it works fine either way). Please let me know if this makes sense to you and I'll include this in the next update.

Regards,
Tim