频谱图#
- ShortTimeFFT.spectrogram(x, y=None, detr=None, *, p0=None, p1=None, k_offset=0, padding='zeros', axis=-1)[source]#
计算频谱图或互频谱图。
频谱图是短时傅里叶变换（STFT）的绝对平方，即对于给定的
S[q,p]，它是abs(S[q,p])**2，因此总是非负的。对于两个STFTSx[q,p], Sy[q,p]，互频谱图定义为Sx[q,p] * np.conj(Sy[q,p])，并且是复数值的。这是一个调用stft/stft_detrend的便捷函数，因此所有参数都在那里讨论。- 参数:
- xnp.ndarray
输入信号，可以是实值或复值数组。对于复数值，属性
fft_mode必须设置为 'twosided' 或 'centered'。- ynp.ndarray
第二个输入信号，与 x 具有相同的形状。如果为
None，则假定为 x。对于复数值，属性fft_mode必须设置为 'twosided' 或 'centered'。- detr‘linear’ | ‘constant’ | Callable[[np.ndarray], np.ndarray] | None
如果设置为 'constant'，则减去均值；如果设置为 'linear'，则从每个段中移除线性趋势。这通过调用
detrend实现。如果 detr 是一个带有一个参数的函数，则 detr 将应用于每个段。对于None（默认），不移除趋势。- p0int | None
要计算的切片范围的第一个元素。如果为
None，则将其设置为p_min，即最小可能的切片。- p1int | None
数组的末尾。如果为
None，则使用 p_max(n)。- k_offsetint
x 中第一个样本（t = 0）的索引。
- padding‘zeros’ | ‘edge’ | ‘even’ | ‘odd’
当滑动窗口超出输入 x 的下端或上端时，添加的值的类型。如果设置为默认的 'zeros'，则添加零。对于 'edge'，使用 x 的第一个或最后一个值。'even' 通过在第一个或最后一个样本上反射信号来填充，'odd' 额外乘以 -1。
- axisint
计算 STFT 时 x 的轴。如果未给出，则使用最后一个轴。
- 返回:
- S_xynp.ndarray
如果
x is y或 y 为None，则返回一个非负实值数组。其维度总是比 x 大一。最后一个轴总是代表频谱图的时间切片。axis 定义频率轴（默认为倒数第二个轴）。例如，对于一维 x，返回一个复杂的二维数组，其中轴 0 代表频率，轴 1 代表时间切片。
另请参阅
stft执行短时傅里叶变换。
stft_detrend从每个段中减去趋势的 STFT。
scipy.signal.ShortTimeFFT此方法所属的类。
说明
互频谱图可以解释为互谱密度的时间-频率模拟量（请参阅
csd）。互频谱图 Sxy 的绝对平方 |Sxy|² 除以频谱图 Sxx 和 Syy 可以解释为相干频谱图Cxy := abs(Sxy)**2 / (Sxx*Syy)，这是coherence的时间-频率模拟量。如果 STFT 被参数化为酉变换，即利用
from_win_equals_dual，则标量积的值以及能量都得以保留。示例
以下示例展示了以 20 Hz 采样、频率 \(f_i(t)\)（在图中用绿色虚线标记）变化的方波的频谱图。所用的高斯窗长为 50 个样本或 2.5 秒。对于
ShortTimeFFT，选择了参数mfft=800（过采样因子 16）和hop间隔 2，以产生足够多的点。由于功率谱密度以 dB 为单位，绘图的颜色映射按对数比例缩放。信号 x 的时间范围由垂直虚线标记，阴影区域表示边界效应的存在。
>>> import matplotlib.pyplot as plt >>> import numpy as np >>> from scipy.signal import square, 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 = 5e-3*(t_x - t_x[N // 3])**2 + 1 # varying frequency >>> x = square(2*np.pi*np.cumsum(f_i)*T_x) # the signal ... >>> g_std = 12 # standard deviation for Gaussian window in samples >>> win = gaussian(50, std=g_std, sym=True) # symmetric Gaussian wind. >>> SFT = ShortTimeFFT(win, hop=2, fs=1/T_x, mfft=800, scale_to='psd') >>> Sx2 = SFT.spectrogram(x) # calculate absolute square of STFT ... >>> 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"Spectrogram ({SFT.m_num*SFT.T:g}$\,s$ Gaussian " + ... rf"window, $\sigma_t={g_std*SFT.T:g}\,$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)) >>> Sx_dB = 10 * np.log10(np.fmax(Sx2, 1e-4)) # limit range to -40 dB >>> im1 = ax1.imshow(Sx_dB, origin='lower', aspect='auto', ... extent=SFT.extent(N), cmap='magma') >>> ax1.plot(t_x, f_i, 'g--', alpha=.5, label='$f_i(t)$') >>> fig1.colorbar(im1, label='Power Spectral Density ' + ... r"$20\,\log_{10}|S_x(t, f)|$ in dB") ... >>> # 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=.3) >>> for t_ in [0, N * SFT.T]: # mark signal borders with vertical line ... ax1.axvline(t_, color='c', linestyle='--', alpha=0.5) >>> ax1.legend() >>> fig1.tight_layout() >>> plt.show()
对数缩放揭示了方波的奇次谐波，这些谐波在 10 Hz 的奈奎斯特频率处反射。这种混叠也是图中噪声伪影的主要来源。