ENH: Clarify scipy.signal.spectrogram documentation (scaling description does not consider mode)

[KDMortimer]

Is your feature request related to a problem? Please describe.

I experienced some confusion when I was trying to compare results from scipy's spectrogram function and a custom spectrogram function I had written. It took some experimentation for me to figure out that while the documentation states the following:

scaling : { ‘density’, ‘spectrum’ }, optional

Selects between computing the power spectral density (‘density’) where Sxx has units of V\*\*2/Hz and computing the power spectrum (‘spectrum’) where Sxx has units of V\*\*2, if x is measured in V and fs is measured in Hz. Defaults to ‘density’.

This is actually only true if mode is set to 'psd'. I had selected 'magnitude' because I naively thought since I was using 'spectrum' and I wanted a power magnitude, 'psd' was inappropriate.

scipy/scipy/signal/spectral.py

Lines 752 to 762 in 47bb6fe

	if mode == 'psd':
	freqs, time, Sxx = _spectral_helper(x, x, fs, window, nperseg,
	noverlap, nfft, detrend,
	return_onesided, scaling, axis,
	mode='psd')
	else:
	freqs, time, Sxx = _spectral_helper(x, x, fs, window, nperseg,
	noverlap, nfft, detrend,
	return_onesided, scaling, axis,
	mode='stft')

scipy/scipy/signal/spectral.py

Lines 1810 to 1811 in 47bb6fe

	if mode == 'stft':
	scale = np.sqrt(scale)

scipy/scipy/signal/spectral.py

Lines 1834 to 1842 in 47bb6fe

	result = _fft_helper(x, win, detrend_func, nperseg, noverlap, nfft, sides)
	if not same_data:
	# All the same operations on the y data
	result_y = _fft_helper(y, win, detrend_func, nperseg, noverlap, nfft,
	sides)
	result = np.conjugate(result) * result_y
	elif mode == 'psd':
	result = np.conjugate(result) * result

For all other modes, _spectral_helper receives the mode parameter stft. This in turn means that when our spectrogram mode is not psd, the output is in V/sqrt(Hz) or V. It also silently ignores the one_sided parameter for scaling, but still omits negative frequency components. There is no mention of units in the mode help:

mode : str, optional

Defines what kind of return values are expected. Options are [‘psd’, ‘complex’, ‘magnitude’, ‘angle’, ‘phase’]. ‘complex’ is equivalent to the output of stft with no padding or boundary extension. ‘magnitude’ returns the absolute magnitude of the STFT. ‘angle’ and ‘phase’ return the complex angle of the STFT, with and without unwrapping, respectively.

From the documentation, I had incorrectly assumed that the Sxx output of scipy.signal.spectrogram always had the units described in the section on scaling. The use of psd as the mode which returns the correct units for the spectrum scaling (one that is explicitly not a density) was also part of the problem.

Describe the solution you'd like.

My preferred solution would be to extend the help for scaling to include reference to the mode.

E.g.,

When the mode is psd, this selects between computing the power spectral density (‘density’) where Sxx has units of V**2/Hz and computing the power spectrum (‘spectrum’) where Sxx has units of V**2, if x is measured in V and fs is measured in Hz. For other modes, this selects between linear spectral density ('density') where Sxx has units of V/Hz**(1/2) and computing the linear spectrum ('spectrum') where Sxx has units of V, if x is measured in V and fs is measured in Hz.

This also has the benefit of not requiring any changes to the underlying math or existing functions.

Describe alternatives you've considered.

Updating mode to use terminology less likely to cause confusion would be another solution. As far as I can tell, psd is actually extracting the absolute square of complex mode, or the square of the magnitude mode, and is relevant for both power spectral density (PSD) and power spectrum (PS) scalings. Thus you could describe it as 'AbsSq'. But without changing the scaling help, some inexperienced users might expect this output to have units of V**4. It would also cause compatibility problems.

Additional context (e.g. screenshots, GIFs)

It might also be useful to mention that return_onesided does not actually scale the results in any way when mode is not psd. The onesided scaling by 2 is correcting for omitted power, and other modes do not actually return power.

If you reverse the math, you could argue that the positive frequencies of a double sided linear spectrum or spectral density could be multiplied by the square root of two to obtain the equivalent magnitude of a single sided linear spectrum or spectral density, but as far as I know this factor of square root of 2 does not have any obvious physical meaning in linear space, compared to the simple 'you lost half the power by omitting negative frequencies' explanation in power space.

[DietBru]

Hi, for reference, I looked into the variation of the parameters: The test signal is a 10 second single cosine with amplitude 10 V. Hence its energy is 500 V²s and its power is 50 V².

The table shows the STFT's maximum and its spectral power depending on the input parameters:

Mode       Scaling    One-sided    max(S)       Power
---------  ---------  -----------  ---------  -------
psd        spectrum   True         50.0          10.8
psd        spectrum   False        25.0          10.8
psd        density    True         232.3         50
psd        density    False        116.1         50
complex    spectrum   True         5.0+0.0j       2
complex    spectrum   False        5.0+0.0j       4
complex    density    True         10.8+0.0j      4.3
complex    density    False        10.8+0.0j      8.6
magnitude  spectrum   True         5.0            2
magnitude  spectrum   False        5.0            4
magnitude  density    True         10.8           4.3
magnitude  density    False        10.8           8.6

Probably only those combinations, where either the power (50 V²) or the amplitude (5 V) for a given frequency can be derived, are useful in practice.

The table can be reproduced by the following code:

from itertools import product

import numpy as np
from scipy import signal
from tabulate import tabulate

n = 500  # number of signal samples
T = 10 / n  # sampling interval in seconds
a = 10  # amplitude

# Create a cosine signal:
f_X = np.fft.rfftfreq(n, T)
X = np.zeros(len(f_X), dtype=complex)
X[len(f_X) // 2] = n * a / 2
x = np.fft.irfft(X)
t_x = np.arange(n) * T
print(f"x is a {n*T} s cosine with amplitude {x[0]}.")
print(f"Energy is {sum(x**2)*T:g} and power is {sum(x**2) / n:g}")


def gen_table():
    scalings = {'density', 'spectrum'}
    modes = ['psd', 'complex', 'magnitude']
    one_sided = (True, False)
    for m_, s_, o_ in product(modes, scalings, one_sided):
        f_S, t_S, S = signal.spectrogram(x, 1 / T, scaling=s_, mode=m_,
                                         return_onesided=o_)
        delta_f = f_S[1] - f_S[0]
        P_V = sum(abs(S[:, len(t_S) // 2])) * delta_f
        yield m_, s_, o_, f"{np.max(S):.1f}", f"{P_V:.1f}"


ss = tabulate(gen_table(), ('Mode', 'Scaling', 'One-sided', 'max(S)', 'Power'))
print(ss)

A PR for extending the documentation will follow.