Python SciPy signal.ShortTimeFFT用法及代码示例

本文简要介绍 python 语言中 scipy.signal.ShortTimeFFT 的用法。

用法:

class  scipy.signal.ShortTimeFFT(win, hop, fs, *, fft_mode='onesided', mfft=None, dual_win=None, scale_to=None, phase_shift=0)#

提供参数化离散Short-time 傅立叶变换 (stft) 及其逆变换 (istft)。

斯夫特通过滑动窗口计算顺序 FFT (win)通过输入信号hop增量。它可用于量化频谱随时间的变化。

斯夫特由 complex-valued 矩阵 S[q,p] 表示，其中 p-th 列表示 FFT，窗口以时间 t[p] = p * 为中心delta_t= p *hop*T其中T是输入信号的采样间隔。 q-th 行表示频率 f[q] = q * 处的值delta_f和delta_f= 1 /(mfft*T) 是 FFT 的 bin 宽度。

逆 STFT伊斯特夫通过反转 STFT 的步骤来计算：对 S[q,p] 的 p-th 切片进行 IFFT，并将结果与所谓的双窗口相乘(参见dual_win)。将结果平移 p *delta_t并将结果添加到先前的移位结果以重建信号。如果只知道双窗口并且 STFT 是可逆的，from_dual可以用来实例化这个类。

由于时间 t = 0 的惯例是在输入信号的第一个样本处，因此 STFT 值通常具有负时隙。因此，像 p_min 或 k_min 这样的负索引并不像标准 Python 索引那样指示从数组末尾向后计数，而是指示在 t = 0 的左侧。

更详细的信息可以在 SciPy 用户指南的 Short-Time 傅里叶变换部分找到。

请注意，除了 scale_to (使用 scaling )之外，初始化程序的所有参数都具有相同的命名属性。

参数 ：：

win： np.ndarray

该窗口必须是实数或 complex-valued 一维数组。

hop： int

样本增量，窗口在每一步中移动。

fs： 浮点数

输入信号和窗口的采样频率。它与采样间隔 T 的关系是T = 1 / fs。

fft_mode： ‘twosided’, ‘centered’, ‘onesided’，'oneside2X'

要使用的 FFT 模式(默认‘onesided’)。有关详细信息，请参阅属性 fft_mode 。

mfft: int | None：

如果需要零填充 FFT，则使用 FFT 的长度。如果None(默认)，则使用窗口 win 的长度。

dual_win： np.ndarray |没有任何

win 的双窗口。如果设置为 None ，则根据需要进行计算。

scale_to： ‘magnitude’, ‘psd’ |没有任何

如果不是 None(默认)，则窗口函数会缩放，因此每个 STFT 列代表 ‘magnitude’ 或功率谱密度 (‘psd’) 频谱。此参数将属性 scaling 设置为相同的值。有关详细信息，请参阅方法 scale_to 。

phase_shift： int |没有任何

如果设置，则为每个频率 f 添加线性相位 phase_shift / mfft * f 。默认值 0 确保第 0 个切片(其中以 t=0 为中心)没有相移。有关更多详细信息，请参阅属性 phase_shift 。

例子：

以下示例显示了具有不同频率 \(f_i(t)\) 的正弦的 STFT 幅度(在图中用红色虚线标记)：

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.signal import ShortTimeFFT
>>> from scipy.signal.windows import gaussian
...
>>> T_x, N = 1 / 20, 1000  # 20 Hz sampling rate for 50 s signal
>>> t_x = np.arange(N) * T_x  # time indexes for signal
>>> f_i = 1 * np.arctan((t_x - t_x[N // 2]) / 2) + 5  # varying frequency
>>> x = np.sin(2*np.pi*np.cumsum(f_i)*T_x) # the signal

所使用的高斯窗口为 50 个样本或 2.5 秒长。参数mfft=200在ShortTimeFFT导致频谱过采样 4 倍：

>>> g_std = 8  # standard deviation for Gaussian window in samples
>>> w = gaussian(50, std=g_std, sym=True)  # symmetric Gaussian window
>>> SFT = ShortTimeFFT(w, hop=10, fs=1/T_x, mfft=200, scale_to='magnitude')
>>> Sx = SFT.stft(x)  # perform the STFT

在图中，信号 x 的时间范围由垂直虚线标记。请注意，SFT 生成的值超出了 x 的时间范围。左侧和右侧的阴影区域表示该区域中的窗口切片未完全位于 x 时间范围内而导致的边框效应：

>>> fig1, ax1 = plt.subplots(figsize=(6., 4.))  # enlarge plot a bit
>>> t_lo, t_hi = SFT.extent(N)[:2]  # time range of plot
>>> ax1.set_title(rf"STFT ({SFT.m_num*SFT.T:g}$\,s$ Gaussian window, " +
...               rf"$\sigma_t={g_std*SFT.T}\,$s)")
>>> ax1.set(xlabel=f"Time $t$ in seconds ({SFT.p_num(N)} slices, " +
...                rf"$\Delta t = {SFT.delta_t:g}\,$s)",
...         ylabel=f"Freq. $f$ in Hz ({SFT.f_pts} bins, " +
...                rf"$\Delta f = {SFT.delta_f:g}\,$Hz)",
...         xlim=(t_lo, t_hi))
...
>>> im1 = ax1.imshow(abs(Sx), origin='lower', aspect='auto',
...                  extent=SFT.extent(N), cmap='viridis')
>>> ax1.plot(t_x, f_i, 'r--', alpha=.5, label='$f_i(t)$')
>>> fig1.colorbar(im1, label="Magnitude $|S_x(t, f)|$")
...
>>> # Shade areas where window slices stick out to the side:
>>> for t0_, t1_ in [(t_lo, SFT.lower_border_end[0] * SFT.T),
...                  (SFT.upper_border_begin(N)[0] * SFT.T, t_hi)]:
...     ax1.axvspan(t0_, t1_, color='w', linewidth=0, alpha=.2)
>>> for t_ in [0, N * SFT.T]:  # mark signal borders with vertical line:
...     ax1.axvline(t_, color='y', linestyle='--', alpha=0.5)
>>> ax1.legend()
>>> fig1.tight_layout()
>>> plt.show()

重建信号伊斯特夫很简单，但请注意长度x1应指定，因为 SFT 长度增加hop脚步：

>>> SFT.invertible  # check if invertible
True
>>> x1 = SFT.istft(Sx, k1=N)
>>> np.allclose(x, x1)
True

可以计算信号部分的 SFT：

>>> p_q = SFT.nearest_k_p(N // 2)
>>> Sx0 = SFT.stft(x[:p_q])
>>> Sx1 = SFT.stft(x[p_q:])

将连续的 STFT 部件组装在一起时，需要考虑重叠：

>>> p0_ub = SFT.upper_border_begin(p_q)[1] - SFT.p_min
>>> p1_le = SFT.lower_border_end[1] - SFT.p_min
>>> Sx01 = np.hstack((Sx0[:, :p0_ub],
...                   Sx0[:, p0_ub:] + Sx1[:, :p1_le],
...                   Sx1[:, p1_le:]))
>>> np.allclose(Sx01, Sx)  # Compare with SFT of complete signal
True

还可以计算信号部分的itsft：

>>> y_p = SFT.istft(Sx, N//3, N//2)
>>> np.allclose(y_p, x[N//3:N//2])
True

属性 ：：

T

输入信号和窗口的采样间隔。

delta_f

STFT 频率仓的宽度。

delta_t

STFT 的时间增量。

dual_win

规范双窗口。

f

STFT 的频率值。

f_pts

沿频率轴的点数。

fac_magnitude

乘以 STFT 值的因子，将每个频率切片缩放为幅度谱。

fac_psd

将 STFT 值乘以将每个频率切片缩放到功率谱密度 (PSD) 的因子。

fft_mode

使用的 FFT 模式(‘twosided’, ‘centered’, ‘onesided’ 或“oneside2X”)。

fs

输入信号和窗口的采样频率。

hop

滑动窗口的信号样本的时间增量。

invertible

检查 STFT 是否可逆。

k_min

STFT 的最小可能信号索引。

lower_border_end

第一信号索引和第一切片索引不受预填充的影响。

m_num

窗口 win 中的样本数。

m_num_mid

窗口 win 的中心索引。

mfft

使用的 FFT 输入长度 - 可能大于窗口长度 m_num 。

onesided_fft

如果使用单侧 FFT，则返回 True。

p_min

尽可能小的切片索引。

phase_shift

如果设置，则将线性相位 phase_shift / mfft * f 添加到频率 f 的每个 FFT 切片。

scaling

标准化应用于窗口函数(‘magnitude’, ‘psd’ 或 None )。

win

窗口函数为实数或 complex-valued 一维数组。