scipy.signal.ShortTimeFFT.

from_win_equals_dual#

classmethod ShortTimeFFT.from_win_equals_dual(desired_win, hop, fs, *, fft_mode='onesided', mfft=None, scale_to=None, phase_shift=0)[source]#

Create instance where the window and its dual are equal up to a scaling factor.

An instance is created were window and dual window are equal as well as being closest to the parameter desired_win in the least-squares sense, i.e., minimizing abs(win-desired_win)**2. Hence, win has the same length as desired_win. Then a scaling factor is applied accoring to the scale_to parameter.

All other parameters have the identical meaning as in the initializer.

To be able to calculate a valid window, desired_win needs to have a valid dual STFT window for the given hop interval. If this is not the case, a ValueError is raised.

Parameters:
desired_winnp.ndarray

A real-valued or complex-valued 1d array containing the sample of the desired window.

hopint

The increment in samples, by which the window is shifted in each step.

fsfloat

Sampling frequency of input signal and window. Its relation to the sampling interval T is T = 1 / fs.

fft_mode‘twosided’, ‘centered’, ‘onesided’, ‘onesided2X’

Mode of FFT to be used (default ‘onesided’). See property fft_mode for details.

mfft: int | None

Length of the FFT used, if a zero padded FFT is desired. If None (default), the length of the window win is used.

scale_to‘magnitude’ | ‘psd’ | ‘unitary’ | None

If not None (default) the window function is scaled, so each STFT column represents either a ‘magnitude’ or a power spectral density (‘psd’) spectrum, Alternatively, the STFT can be scaled to a`unitary` mapping, i.e., dividing the window by np.sqrt(mfft) and multiplying the dual window by the same amount.

phase_shiftint | None

If set, add a linear phase phase_shift / mfft * f to each frequency f. The default value of 0 ensures that there is no phase shift on the zeroth slice (in which t=0 is centered). See property phase_shift for more details.

See also

closest_STFT_dual_window

Calculate the STFT dual window of a given window closest to a desired dual window.

ShortTimeFFT.spectrogram

Calculate squared STFTs

ShortTimeFFT

Class this property belongs to.

Notes

The set of all possible windows with identical dual is defined by the set of linear constraints of Eq. (24) in the Short-Time Fourier Transform section of the SciPy User Guide. There it is also derived that ShortTimeFFT.dual_win == ShortTimeFFT.m_pts * ShortTimeFFT.win needs to hold for an STFT to be a unitary mapping.

A unitary mapping preserves the value of the scalar product, i.e.,

\[\langle x, y\rangle = \sum_k x[k]\, \overline{y[k]} \stackrel{\stackrel{\text{unitary}}{\downarrow}}{=} \sum_{q,p} S_x[q,p]\, \overline{S_y[q,p]} = \langle S_x[q,p], S_y[q,p]\rangle\ ,\]

with \(S_{x,y}\) being the STFT of \(x,y\). Hence, the energy \(E_x=T\sum_k |x[k]|^2\) of a signal is also preserved. This is also illustrated in the example below.

Thie reason of distinguishing between no scaling (i.e., parameter scale_to is None) and unitary scaling (i.e., scale_to = 'unitary') is due to the utilized FFT function not being unitary (i.e., using the default value 'backward' for the fft parameter norm).

Examples

The following example shows that an STFT can be indeed unitary:

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from scipy.signal import ShortTimeFFT, windows
...
>>> m, hop, std = 36, 8, 5
>>> desired_win = windows.gaussian(m, std, sym=True)
>>> SFT = ShortTimeFFT.from_win_equals_dual(desired_win, hop, fs=1/m,
...                                         fft_mode='twosided',
...                                         scale_to='unitary')
>>> np.allclose(SFT.dual_win, SFT.win * SFT.m_num)  # check if STFT is unitary
True
>>> x1, x2 = np.tile([-1, -1, 1, 1], 5), np.tile([1, -1, -1, 1], 5)
>>> np.sum(x1*x2) # scalar product is zero -> orthogonal signals
0
>>> np.sum(x1**2)  # scalar product of x1 with itself
20
>>> Sx11, Sx12 = SFT.spectrogram(x1), SFT.spectrogram(x1, x2)
>>> np.sum(Sx12)  # STFT scalar product is also zero
-4.163336342344337e-16+0j  # may vary
>>> np.sum(Sx11)  # == np.sum(x1**2)
19.999999999999996  # may vary
...
... # Do the plotting:
>>> fg1, (ax11, ax12) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 4))
>>> s_fac = np.sqrt(SFT.mfft)
>>> _ = fg1.suptitle(f"Scaled Unitary Window of {m} Sample Gaussian with " +
...                  rf"{hop=}, $\sigma={std}$, Scale factor: {s_fac}")
>>> ax11.set(ylabel="Amplitude", xlabel="Samples", xlim=(0, m))
>>> ax12.set(xlabel="Frequency Bins", ylabel="Magnitude Spectrum",
...          xlim=(0, 15), ylim=(1e-5, 1.5))
>>> u_win_str = rf"Unitary $\times{s_fac:g}$"
>>> for x_, n_ in zip((desired_win, SFT.win*s_fac), ('Desired', u_win_str)):
...     ax11.plot(x_, '.-', alpha=0.5, label=n_)
...     X_ = np.fft.rfft(x_) / np.sum(abs(x_))
...     ax12.semilogy(abs(X_), '.-', alpha=0.5, label=n_)
>>> for ax_ in (ax11, ax12):
...     ax_.grid(True)
...     ax_.legend()
>>> plt.show()
../../_images/scipy-signal-ShortTimeFFT-from_win_equals_dual-1_00_00.png

Note that fftmode='twosided' is used, since we need sum over the complete time frequency plane. Due to passing scale_to='unitary' the window SFT.win is scaled by 1/np.sqrt(SFT.mfft). Hence, SFT.win needs to be scaled by s_fac in the plot above.