Scipy 频谱图与多个 Numpy FFT

Question

我正在尝试优化给出的一些代码，其中 FFT 取自时间序列（以列表形式给出）的滑动 window，并且每个结果都累积到一个列表中。原代码如下：

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

总而言之，旧代码 windows 时间序列，删除均值，应用 Hann windowing 函数，采用 FFT（同时零填充，如 ZERO_PAD_LEN>WINDOW_LEN ), 然后取实数一半的绝对值并平方得到功率谱 (V**2 的单位)。然后它将 window 移动 WINDOW_STRIDE_LEN，并重复该过程，直到 window 超出数据末尾。这有 WINDOW_OVERLAP.

的重叠

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

我在这里错过了什么？

Answer 1

缩放

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

另一方面，实现 calc_new 使用 scipy.signal.spectrogram which ultimately uses np.fft.rfft，但也根据收到的 scaling 和 return_onesided 参数缩放结果。更具体地说：

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

因此，您可能期望新实施 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 个频率仓（有 nperseg-2 个具有相应对称性的仓加上 0Hz 处的 2 个不对称仓和奈奎斯特率，所以保留非冗余部分留给你(nperseg-2)//2 + 2 == nperseg//2 + 1）。那就是scipy.signal.spectrogramreturns.

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

时间步数

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

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

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

这种小的纵向分歧可能是从哪里产生的？谢谢！

Scipy 频谱图与多个 Numpy FFT

Scipy Spectrogram vs. multiple Numpy FFT's

python

numpy

scipy

fft

spectrogram