当前位置: 首页>>代码示例 >>用法及示例精选 >>正文


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 窗口数(订单 0Kmax-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()
scipy-signal-windows-dpss-1_00_00.png

以下是前四个窗口的示例,以及它们的浓度比:

>>> 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()
scipy-signal-windows-dpss-1_01_00.png

使用标准\(l_{\infty}\) norm 会为偶数产生两个统一值M, 但奇数只有一个单位值M.这会产生不均匀的窗口功率,可以通过近似校正来抵消M**2/float(M**2+NW), 可以使用norm='approximate'(这与norm=NoneKmax=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()
scipy-signal-windows-dpss-1_02_00.png

相关用法


注:本文由纯净天空筛选整理自scipy.org大神的英文原创作品 scipy.signal.windows.dpss。非经特殊声明,原始代码版权归原作者所有,本译文未经允许或授权,请勿转载或复制。