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 一維數組。