Scipy的谱图与多个Numpy FFT的比较

Question

我试图优化一些我获得的代码，其中对时间序列（以列表形式给出）进行滑动窗口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）），每个频率bin中的值相差很大，有时相差几个数量级。

我错在哪里了？

- ESilk

---

一个结果是另一个的缩放版本吗？可能存在归一化方面的差异。但这并不重要，重要的是输出的形状。 - Cris Luengo

Answer 1

扩展性

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

另一方面，calc_new的实现使用scipy.signal.spectrogram，最终使用np.fft.rfft，但也根据接收到的scaling和return_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。

频率bin的数量

对于偶数点nperseg，您的初始实现calc_old仅保留nperseg//2个频率bin。

另一方面，完整的非冗余半谱应该给出nperseg//2 + 1个频率bin（有nperseg-2个带有对称性的bin加上0Hz和Nyquist速率的两个非对称bin，所以保留非冗余部分会让你得到(nperseg-2)//2 + 2 == nperseg//2 + 1）。这是scipy.signal.spectrogram返回的内容。

换句话说，您的初始实现calc_old缺少Nyquist频率bin。

时间步数

calc_old实现在最后一个时间步跳过了少于bf.WINDOW_LEN个样本的时间步计算。只有当len(raw_data)-bf.WINDOW_STRIDE_LEN是bf.WINDOW_LEN的精确倍数时，才不会跳过这些样本。我猜测这在您特定的输入序列中并不成立。

相比之下，scipy.signal.spectrogram会填充数据以使用所有输入样本进行频谱图计算，这可能会导致与您的calc_old实现相比多出一个时间步。

- SleuthEye

---

这证实了我和同事怀疑并一直在调查的所有内容，特别是时间步骤填充，我没有在任何地方看到有记录？现在我们得到了匹配的结果（而且速度更快），谢谢！ - ESilk

Answer 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')

这个小竖直差异是从哪里来的呢？谢谢！

- Valentin Melnikov