Scipy 频谱图与多个 Numpy FFT

Posted

技术标签:

【中文标题】Scipy 频谱图与多个 Numpy FFT【英文标题】:Scipy Spectrogram vs. multiple Numpy FFT's 【发布时间】:2019-01-24 16:51:46 【问题描述】:

我正在尝试优化我给出的一些代码,其中 FFT 取自时间序列上的滑动窗口(作为列表给出),并且每个结果都累积到一个列表中。原代码如下:

def calc_old(raw_data):
    FFT_old = list()

    for i in range(0, len(raw_data), bf.WINDOW_STRIDE_LEN):
        if (i + bf.WINDOW_LEN) >= len(raw_data):
            # Skip the windows that would extend beyond the end of the data
            continue

        data_tmp = raw_data[i:i+bf.WINDOW_LEN]
        data_tmp -= np.mean(data_tmp)
        data_tmp = np.multiply(data_tmp, np.hanning(len(data_tmp)))
        fft_data_tmp = np.fft.fft(data_tmp, n=ZERO_PAD_LEN)
        fft_data_tmp = abs(fft_data_tmp[:int(len(fft_data_tmp)/2)])**2
        FFT_old.append(fft_data_tmp)

还有新代码:

def calc_new(raw_data):
    data = np.array(raw_data)  # Required as the data is being handed in as a list
    f, t, FFT_new = spectrogram(data,
                                fs=60.0,
                                window="hann",
                                nperseg=bf.WINDOW_LEN,
                                noverlap=bf.WINDOW_OVERLAP,
                                nfft=bf.ZERO_PAD_LEN,
                                scaling='spectrum')

总而言之,旧代码对时间序列进行窗口化,去除均值,应用 Hann 窗口函数,采用 FFT(在零填充时,如 ZERO_PAD_LEN>WINDOW_LEN),然后取实数一半的绝对值和将其平方以产生功率谱(V**2 的单位)。然后它将窗口移动WINDOW_STRIDE_LEN,并重复该过程,直到窗口超出数据的末尾。这与WINDOW_OVERLAP 有重叠。

据我所知,频谱图应该对我给出的论点做同样的事情。然而,每个轴的 FFT 的结果维度相差 1(例如,旧代码是 MxN,新代码是 (M+1)x(N+1)),并且每个频率仓中的值大不相同——几个数量级数量级,在某些情况下。

我在这里错过了什么?

【问题讨论】:

一个结果是另一个结果的缩放版本吗?标准化可能存在差异。但这并不重要,重要的是输出的形状。 【参考方案1】:

缩放

calc_old 中的实现直接使用来自np.fft.fft 的输出,没有任何缩放。

另一方面,实现calc_new 使用scipy.signal.spectrogram,它最终使用np.fft.rfft,但也根据收到的scalingreturn_onesided 参数缩放结果。更具体地说:

对于默认的 return_onesided=True(因为您没有在 calc_new 中提供明确的值),每个 bin 的值加倍以计算包括对称 bin 在内的总能量。 对于提供的scaling='spectrum',这些值将进一步按系数1.0/win.sum()**2 缩放。对于选定的 Hann 窗口,对应于 4/N**2,其中 N=bf.WINDOW_LEN 是窗口长度。

因此,您可能期望新的实现 cald_new 为您提供一个与 calc_old 相比按整体因子 8/bf.WINDOW_LEN**2 缩放的结果。或者,如果您希望您的第二个实现提供与calc_old 相同的缩放比例,您应该将scipy.signal.spectrogram 的结果乘以0.125 * bf.WINDOW_LEN**2

频率箱数

给定偶数点 nperseg,您的初始实现 calc_old 仅保留 nperseg//2 频率区间。

另一方面,完整的非冗余半频谱应该为您提供nperseg//2 + 1 频率 bin(有 nperseg-2 bin 具有相应的对称性加上 2 个 0Hz 和 Nyquist 速率的不对称 bin,因此保留非冗余部分留给你(nperseg-2)//2 + 2 == nperseg//2 + 1)。这就是scipy.signal.spectrogram 返回的内容。

换句话说,您的初始实现 calc_old 缺少奈奎斯特频率仓。

时间步数

如果留给最后一个时间步计算的样本少于bf.WINDOW_LEN,则calc_old 中的实现会跳过最后一个时间步。它不会跳过这些样本,仅当len(raw_data)-bf.WINDOW_STRIDE_LENbf.WINDOW_LEN。我猜你的特定输入序列不是这种情况。

相比之下,scipy.signal.spectrogram 会在需要时使用额外样本填充数据,以便在频谱图计算期间使用所有输入样本,与您的 calc_old 实现相比,这可能会导致额外的时间步长。

【讨论】:

这证实了我和我的同事所怀疑和正在调查的一切,尤其是时间步长填充,我在任何地方都没有看到记录?我们现在得到了匹配的结果(而且速度更快),谢谢!【参考方案2】:

可能有人知道原因,为什么在将频谱图结果与手动 FFT 进行比较时会出现一些差异?

# Parametrs of signal and its preprocessing
sample_rate = 4
window_size = 512 * sample_rate
detrend = 'linear'
tukey_alpha = 0.25
[![enter image description here][1]][1]
# Spectrogram
f, t, S = scipy.signal.spectrogram(signal, sample_rate, nperseg=window_size, noverlap=sample_rate, scaling='spectrum', mode='magnitude', detrend=detrend)

# FFT on the leftmost window of signal
windowed_signal = signal[:window_size]
windowed_signal = scipy.signal.detrend(windowed_signal, type=detrend)
windowed_signal *= scipy.signal.windows.tukey(window_size, tukey_alpha)
A = np.fft.rfft(windowed_signal)
freqs = np.fft.rfftfreq(window_size) * sample_rate
positive_freqs_n = int(np.ceil(window_size / 2.))
freqs_slice = slice(0, positive_freqs_n)
magnitudes = np.abs(A)[freqs_slice] / window_size

# Plotting
plt.plot(f, S.T[0], label='First window from scipy.signal.spectrogram')
plt.plot(freqs[freqs_slice], magnitudes, alpha=0.5, label='np.fft on first window')
plt.yscale('log')
plt.xscale('log')
plt.legend()
plt.xlabel('Frequency, Hz')
plt.ylabel('Magnitude')

这种小的纵向分歧可能来自哪里? 谢谢!

【讨论】:

以上是关于Scipy 频谱图与多个 Numpy FFT的主要内容,如果未能解决你的问题,请参考以下文章

为啥 scipy 和 numpy fft 图看起来不同?

如何使用 python 获得像 audacity 这样的音频频谱?

来自近乎完美的正弦曲线的奇怪 FFT 频谱

FFT频谱分析(补零频谱泄露栅栏效应加窗细化频谱混叠),MatlabC语言代码

关于用FFT分析信号频谱的问题

FFT 频谱未正确显示