Python SciPy signal.istft用法及代碼示例

本文簡要介紹 python 語言中 scipy.signal.istft 的用法。

用法:

scipy.signal.istft(Zxx, fs=1.0, window='hann', nperseg=None, noverlap=None, nfft=None, input_onesided=True, boundary=True, time_axis=-1, freq_axis=-2, scaling='spectrum')#

執行逆短時傅裏葉變換 (iSTFT)。

遺產

此函數被視為遺留函數，將不再接收更新。這也可能意味著它將在未來的 SciPy 版本中被刪除。 ShortTimeFFT 是一個較新的 STFT /ISTFT 實現，具有更多函數。實現之間的比較可以在 SciPy 用戶指南的 Short-Time 傅裏葉變換部分找到。

參數 ：：

Zxx： array_like

要重建的信號的 STFT。如果傳遞的是純實數數組，它將被強製轉換為複雜數據類型。

fs： 浮點數，可選

時間序列的采樣頻率。默認為 1.0。

window： str 或 tuple 或 數組，可選

想要使用的窗口。如果窗戶是一個字符串或元組，它被傳遞給scipy.signal.get_window生成窗口值，默認為DFT-even。看scipy.signal.get_window獲取窗口列表和所需參數。如果窗戶是數組，它將直接用作窗口，其長度必須為nperseg。默認為 Hann 窗口。必須與用於生成 STFT 的窗口匹配以實現忠實的反演。

nperseg： 整數，可選

每個 STFT 段對應的數據點數。如果每段的數據點數是奇數，或者如果 STFT 是通過填充的，則必須指定此參數nfft > nperseg.如果None，值取決於形狀Zxx和input_onesided.如果input_onesided是True,nperseg=2*(Zxx.shape[freq_axis] - 1).否則，nperseg=Zxx.shape[freq_axis].默認為None.

noverlap： 整數，可選

段之間重疊的點數。如果沒有，則為段長度的一半。默認為無。指定時，必須滿足 COLA 約束(請參閱下麵的注釋)，並且應與用於生成 STFT 的參數匹配。默認為無。

nfft： 整數，可選

每個 STFT 段對應的 FFT 點數。如果 STFT 通過以下方式填充，則必須指定此參數nfft > nperseg.如果None, 默認值與nperseg，上麵詳述，但有一個異常：如果input_onesided是真的並且nperseg==2*Zxx.shape[freq_axis] - 1,nfft也具有該值。這種情況允許使用 odd-length 未填充的 STFT 正確反轉nfft=None.默認為None.

input_onesided： 布爾型，可選

如果True，將輸入數組解釋為單邊 FFT，例如由scipy.signal.stft和return_onesided=True和numpy.fft.rfft.如果False，將輸入解釋為雙邊 FFT。默認為True.

boundary： 布爾型，可選

指定輸入信號是否在其邊界處通過提供非擴展None boundary參數scipy.signal.stft.默認為True.

time_axis： 整數，可選

STFT的時間段所在的位置；默認值為最後一個軸(即 axis=-1 )。

freq_axis： 整數，可選

STFT的頻率軸所在的位置；默認值為倒數第二個軸(即 axis=-2 )。

scaling: {‘spectrum’, ‘psd’}：

默認的 ‘spectrum’ 縮放允許每個頻率線Zxx被解釋為幅度譜。 ‘psd’ 選項將每條線縮放到功率譜密度 - 它允許通過數字積分來計算信號能量abs(Zxx)**2.

返回 ：：

t： ndarray

輸出數據時間數組。

x： ndarray

Zxx 的 iSTFT。

注意：

為了通過 istft 實現 STFT 的反轉，信號加窗必須遵守 “nonzero overlap add” (NOLA) 的約束：

\[\sum_{t}w^{2}[n-tH] \ne 0\]

這確保了出現在 overlap-add 重建方程的分母中的歸一化因子

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

不為零。可以使用 check_NOLA 函數檢查 NOLA 約束。

已修改(通過掩蔽或其他方式)的 STFT 不能保證對應於完全可實現的信號。該函數通過[2]中詳述的最小二乘估計算法實現 iSTFT，該算法產生的信號可最小化返回信號的 STFT 與修改後的 STFT 之間的均方誤差。

參考：

[1]

Oppenheim、Alan V.、Ronald W. Schafer、John R. Buck “Discrete-Time 信號處理”，Prentice Hall，1999 年。

[2]

Daniel W. Griffin, Jae S. Lim “信號估計來自修改的 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()

生成一個測試信號，一個 50Hz 的 2 Vrms 正弦波，被 1024 Hz 采樣的 0.001 V**2/Hz 白噪聲破壞。

>>> fs = 1024
>>> N = 10*fs
>>> nperseg = 512
>>> amp = 2 * np.sqrt(2)
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / float(fs)
>>> carrier = amp * np.sin(2*np.pi*50*time)
>>> noise = rng.normal(scale=np.sqrt(noise_power),
...                    size=time.shape)
>>> x = carrier + noise

計算 STFT，並繪製其大小

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

將小於或等於載波幅度 10% 的分量歸零，然後通過逆 STFT 轉換回時間序列

>>> Zxx = np.where(np.abs(Zxx) >= amp/10, Zxx, 0)
>>> _, xrec = signal.istft(Zxx, fs)

將清潔後的信號與原始和真實的載波信號進行比較。

>>> plt.figure()
>>> plt.plot(time, x, time, xrec, time, carrier)
>>> plt.xlim([2, 2.1])
>>> plt.xlabel('Time [sec]')
>>> plt.ylabel('Signal')
>>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
>>> plt.show()

請注意，清理後的信號不會像原始信號那樣突然開始，因為還刪除了一些瞬態係數：

>>> plt.figure()
>>> plt.plot(time, x, time, xrec, time, carrier)
>>> plt.xlim([0, 0.1])
>>> plt.xlabel('Time [sec]')
>>> plt.ylabel('Signal')
>>> plt.legend(['Carrier + Noise', 'Filtered via STFT', 'True Carrier'])
>>> plt.show()