scipy.signal.

stft#

scipy.signal.stft(x, fs=1.0, window='hann', nperseg=256, noverlap=None, nfft=None, detrend=False, return_onesided=True, boundary='zeros', padded=True, axis=-1, scaling='spectrum')[source]#

計算短時傅立葉轉換 (legacy function)。

STFT 可用來量化非平穩訊號的頻率和相位內容隨時間的變化。

舊版

此函數被視為舊版，將不再接收更新。雖然我們目前沒有移除它的計畫，但我們建議新程式碼改用更現代的替代方案。ShortTimeFFT 是一個更新的 STFT / ISTFT 實作，具有更多功能。關於實作之間的比較，可以在短時傅立葉轉換章節的 SciPy 使用者指南中找到。

參數:

xarray_like: 時間序列的量測值
fsfloat, optional: x 時間序列的取樣頻率。預設值為 1.0。
windowstr 或 tuple 或 array_like, optional: 想要使用的視窗。如果 window 是字串或 tuple，則會將其傳遞給 get_window 以產生視窗值，預設情況下這些值是 DFT 偶數。有關視窗和所需參數的列表，請參閱 get_window。如果 window 是 array_like，它將直接用作視窗，且其長度必須為 nperseg。預設值為 Hann 視窗。
npersegint, optional: 每個區段的長度。預設值為 256。
noverlapint, optional: 區段之間重疊的點數。如果 None，noverlap = nperseg // 2。預設值為 None。當指定時，必須滿足 COLA 約束（請參閱下面的「註解」）。
nfftint, optional: FFT 的長度，如果需要零填充 FFT。如果 None，則 FFT 長度為 nperseg。預設值為 None。
detrendstr 或 function 或 False, optional: 指定如何對每個區段進行去趨勢。如果 detrend 是字串，則會將其作為 type 參數傳遞給 detrend 函數。如果它是一個函數，它會接收一個區段並返回一個去趨勢的區段。如果 detrend 是 False，則不進行去趨勢。預設值為 False。
return_onesidedbool, optional: 如果 True，則為實數資料傳回單邊頻譜。如果 False，則傳回雙邊頻譜。預設值為 True，但對於複數資料，始終傳回雙邊頻譜。
boundarystr 或 None, optional: 指定是否在輸入訊號的兩端擴展，以及如何產生新值，以便將第一個加窗區段置於第一個輸入點的中心。這樣做的好處是，當使用的視窗函數從零開始時，可以重建第一個輸入點。有效選項為 ['even', 'odd', 'constant', 'zeros', None]。預設值為 'zeros'，用於零填充擴展。例如 [1, 2, 3, 4] 會擴展為 [0, 1, 2, 3, 4, 0]，當 nperseg=3 時。
paddedbool, optional: 指定是否在輸入訊號的末尾進行零填充，使訊號正好適合整數個視窗區段，以便輸出中包含訊號的所有部分。預設值為 True。如果在邊界擴展之後進行填充，如果 boundary 不是 None，且 padded 為 True，如同預設值。
axisint, optional: 計算 STFT 的軸；預設值為最後一個軸（即 axis=-1）。
scaling: {‘spectrum’, ‘psd’}: 預設的 'spectrum' 縮放允許將 Zxx 的每個頻率線解釋為幅度頻譜。'psd' 選項將每條線縮放到功率頻譜密度 - 它允許通過數值積分 abs(Zxx)**2 來計算訊號的能量。

在版本 1.9.0 中新增。

返回:

fndarray: 取樣頻率陣列。
tndarray: 區段時間陣列。
Zxxndarray: x 的 STFT。預設情況下，Zxx 的最後一個軸對應於區段時間。

另請參閱

istft: 反短時傅立葉轉換
ShortTimeFFT: 提供更多功能的更新 STFT/ISTFT 實作。
check_COLA: 檢查是否滿足恆定重疊相加 (COLA) 約束
check_NOLA: 檢查是否滿足非零重疊相加 (NOLA) 約束
welch: Welch 方法的功率頻譜密度。
spectrogram: Welch 方法的頻譜圖。
csd: Welch 方法的交叉頻譜密度。
lombscargle: 用於不均勻取樣資料的 Lomb-Scargle 週期圖

註解

為了能夠通過 istft 中的反 STFT 反轉 STFT，訊號加窗必須遵守「非零重疊相加」(NOLA) 約束，並且輸入訊號必須具有完整的加窗覆蓋率（即 (x.shape[axis] - nperseg) % (nperseg-noverlap) == 0）。padded 參數可用於實現此目的。

給定時域訊號 \(x[n]\)、視窗 \(w[n]\) 和跳躍大小 \(H\) = nperseg - noverlap，時間索引 \(t\) 的加窗幀由下式給出

\[x_{t}[n]=x[n]w[n-tH]\]

重疊相加 (OLA) 重建方程式由下式給出

\[x[n]=\frac{\sum_{t}x_{t}[n]w[n-tH]}{\sum_{t}w^{2}[n-tH]}\]

NOLA 約束確保出現在 OLA 重建方程式分母中的每個正規化項均為非零。可以使用 check_NOLA 測試 window、nperseg 和 noverlap 的選擇是否滿足此約束。

在版本 0.19.0 中新增。

參考文獻

[1]

Oppenheim, Alan V., Ronald W. Schafer, John R. Buck “Discrete-Time Signal Processing”, Prentice Hall, 1999.

[2]

Daniel W. Griffin, Jae S. Lim “Signal Estimation from Modified Short-Time Fourier Transform”, IEEE 1984, 10.1109/TASSP.1984.1164317

範例

>>> import numpy as np
>>> from scipy import signal
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng()

產生一個測試訊號，一個 2 Vrms 正弦波，其頻率在 3kHz 附近緩慢調變，並受到指數衰減幅度的白雜訊干擾，以 10 kHz 取樣。

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2 * np.sqrt(2)
>>> noise_power = 0.01 * fs / 2
>>> time = np.arange(N) / float(fs)
>>> mod = 500*np.cos(2*np.pi*0.25*time)
>>> carrier = amp * np.sin(2*np.pi*3e3*time + mod)
>>> noise = rng.normal(scale=np.sqrt(noise_power),
...                    size=time.shape)
>>> noise *= np.exp(-time/5)
>>> x = carrier + noise

計算並繪製 STFT 的幅度。

>>> f, t, Zxx = signal.stft(x, fs, nperseg=1000)
>>> plt.pcolormesh(t, f, np.abs(Zxx), vmin=0, vmax=amp, shading='gouraud')
>>> plt.title('STFT Magnitude')
>>> plt.ylabel('Frequency [Hz]')
>>> plt.xlabel('Time [sec]')
>>> plt.show()

../../_images/scipy-signal-stft-1_00_00.png

比較訊號 x 的能量與其 STFT 的能量

>>> E_x = sum(x**2) / fs  # Energy of x
>>> # Calculate a two-sided STFT with PSD scaling:
>>> f, t, Zxx = signal.stft(x, fs, nperseg=1000, return_onesided=False,
...                         scaling='psd')
>>> # Integrate numerically over abs(Zxx)**2:
>>> df, dt = f[1] - f[0], t[1] - t[0]
>>> E_Zxx = sum(np.sum(Zxx.real**2 + Zxx.imag**2, axis=0) * df) * dt
>>> # The energy is the same, but the numerical errors are quite large:
>>> np.isclose(E_x, E_Zxx, rtol=1e-2)
True