Python SciPy windows.dpss用法及代碼示例

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

用法:

scipy.signal.windows.dpss(M, NW, Kmax=None, sym=True, norm=None, return_ratios=False)#

計算離散長球體序列 (DPSS)。

DPSS(或 Slepian 序列)通常用於多錐功率譜密度估計(參見 [1])。序列中的第一個窗口可用於最大化主瓣中的能量集中，也稱為Slepian窗口。

參數 ：：

M： int

窗口長度。

NW： 浮點數

對應於2*NW = BW/f0 = BW*M*dt的標準化半帶寬，其中dt取為1。

Kmax： int |無，可選

要返回的 DPSS 窗口數(訂單 0 到 Kmax-1 )。如果無(默認)，則僅返回一個形狀為 (M,) 的窗口，而不是形狀為 (Kmax, M) 的窗口數組。

sym： 布爾型，可選

當為 True(默認)時，生成一個對稱窗口，用於濾波器設計。當為 False 時，生成一個周期窗口，用於頻譜分析。

norm： {2，‘approximate’, ‘subsample’} |無，可選

如果‘approximate’或‘subsample’，則窗口按最大值歸一化，並且使用M**2/(M**2+NW)(“approximate”)或基於FFT的子樣本移位(“subsample”)，詳情請參閱注釋。如果無，則在 Kmax=None 時使用 “approximate”，否則使用 2(使用 l2 範數)。

return_ratios： 布爾型，可選

如果為 True，除了窗口之外，還返回濃度比率。

返回 ：：

v： ndarray，形狀(Kmax，M)或(M，)

DPSS 窗口。如果 Kmax 為無，則為一維。

r： ndarray，形狀(Kmax，)或浮點數，可選

窗口的濃度比。僅在 return_ratios 計算結果為 True 時返回。如果 Kmax 為無，則為 0D。

注意：

此計算使用 [2] 中給出的三對角特征向量公式。

Kmax=None 的默認歸一化，即窗口生成模式，隻需使用 l-infinity 範數將創建一個具有兩個單位值的窗口，這會在偶數和奇數階之間產生輕微的歸一化差異。對偶數樣本數的 M**2/float(M**2+NW) 進行近似校正用於抵消這種影響(請參見下麵的示例)。

對於非常長的信號(例如，1e6 個元素)，計算窗口數量級更短並使用插值(例如，scipy.interpolate.interp1d) 以獲得長度的錐度M，但這通常不會保持錐度之間的正交性。

參考：

[1]

Percival DB，Walden WT。物理應用的頻譜分析：多錐度和傳統的單變量技術。劍橋大學出版社； 1993 年。

[2]

Slepian, D. Prolate 球麵波函數、傅裏葉分析和不確定性 V：離散情況。貝爾係統技術雜誌，第 57 卷(1978 年)，1371430。

[3]

凱撒，JF，謝弗 RW。關於使用 I0-Sinh 窗口進行頻譜分析。 IEEE 聲學、語音和信號處理匯刊。 ASSP-28(1)：105-107； 1980 年。

例子：

我們可以將該窗口與 kaiser 進行比較，它是作為更容易計算的替代方案而發明的 [3](示例改編自 這裏 )：

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.signal import windows, freqz
>>> M = 51
>>> fig, axes = plt.subplots(3, 2, figsize=(5, 7))
>>> for ai, alpha in enumerate((1, 3, 5)):
...     win_dpss = windows.dpss(M, alpha)
...     beta = alpha*np.pi
...     win_kaiser = windows.kaiser(M, beta)
...     for win, c in ((win_dpss, 'k'), (win_kaiser, 'r')):
...         win /= win.sum()
...         axes[ai, 0].plot(win, color=c, lw=1.)
...         axes[ai, 0].set(xlim=[0, M-1], title=r'$\alpha$ = %s' % alpha,
...                         ylabel='Amplitude')
...         w, h = freqz(win)
...         axes[ai, 1].plot(w, 20 * np.log10(np.abs(h)), color=c, lw=1.)
...         axes[ai, 1].set(xlim=[0, np.pi],
...                         title=r'$\beta$ = %0.2f' % beta,
...                         ylabel='Magnitude (dB)')
>>> for ax in axes.ravel():
...     ax.grid(True)
>>> axes[2, 1].legend(['DPSS', 'Kaiser'])
>>> fig.tight_layout()
>>> plt.show()

以下是前四個窗口的示例，以及它們的濃度比：

>>> M = 512
>>> NW = 2.5
>>> win, eigvals = windows.dpss(M, NW, 4, return_ratios=True)
>>> fig, ax = plt.subplots(1)
>>> ax.plot(win.T, linewidth=1.)
>>> ax.set(xlim=[0, M-1], ylim=[-0.1, 0.1], xlabel='Samples',
...        title='DPSS, M=%d, NW=%0.1f' % (M, NW))
>>> ax.legend(['win[%d] (%0.4f)' % (ii, ratio)
...            for ii, ratio in enumerate(eigvals)])
>>> fig.tight_layout()
>>> plt.show()

使用標準\(l_{\infty}\)norm 會為偶數產生兩個統一值M, 但奇數隻有一個單位值M.這會產生不均勻的窗口功率，可以通過近似校正來抵消M**2/float(M**2+NW), 可以使用norm='approximate'(這與norm=None當Kmax=None，就像這裏的情況一樣)。或者，速度越慢norm='subsample'可以使用，它使用頻域中的子樣本移位 (FFT) 來計算校正：

>>> Ms = np.arange(1, 41)
>>> factors = (50, 20, 10, 5, 2.0001)
>>> energy = np.empty((3, len(Ms), len(factors)))
>>> for mi, M in enumerate(Ms):
...     for fi, factor in enumerate(factors):
...         NW = M / float(factor)
...         # Corrected using empirical approximation (default)
...         win = windows.dpss(M, NW)
...         energy[0, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
...         # Corrected using subsample shifting
...         win = windows.dpss(M, NW, norm='subsample')
...         energy[1, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
...         # Uncorrected (using l-infinity norm)
...         win /= win.max()
...         energy[2, mi, fi] = np.sum(win ** 2) / np.sqrt(M)
>>> fig, ax = plt.subplots(1)
>>> hs = ax.plot(Ms, energy[2], '-o', markersize=4,
...              markeredgecolor='none')
>>> leg = [hs[-1]]
>>> for hi, hh in enumerate(hs):
...     h1 = ax.plot(Ms, energy[0, :, hi], '-o', markersize=4,
...                  color=hh.get_color(), markeredgecolor='none',
...                  alpha=0.66)
...     h2 = ax.plot(Ms, energy[1, :, hi], '-o', markersize=4,
...                  color=hh.get_color(), markeredgecolor='none',
...                  alpha=0.33)
...     if hi == len(hs) - 1:
...         leg.insert(0, h1[0])
...         leg.insert(0, h2[0])
>>> ax.set(xlabel='M (samples)', ylabel=r'Power / $\sqrt{M}$')
>>> ax.legend(leg, ['Uncorrected', r'Corrected: $\frac{M^2}{M^2+NW}$',
...                 'Corrected (subsample)'])
>>> fig.tight_layout()