我正在尝试优化给出的一些代码,其中 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))并且每个频率仓中的值有很大不同——几个数量级在某些情况下,数量级。
我在这里错过了什么?
最佳答案
缩放
calc_old
中的实现使用 np.fft.fft
的输出直接没有任何缩放。
另一方面,实现 calc_new
使用 scipy.signal.spectrogram
最终使用 np.fft.rfft
而且还根据接收到的 scaling
和 return_onesided
参数缩放结果。更具体地说:
- 对于默认的
return_onesided=True
(因为您没有在calc_new
中提供明确的值),每个 bin 的值加倍以计算总能量,包括对称容器。 - 对于所提供的
scaling='spectrum'
,这些值进一步按1.0/win.sum()**2
因子缩放。对于选定的 Hann 窗口,对应于4/N**2
,其中N=bf.WINDOW_LEN
是窗口长度。
因此,您可能期望新实现 cald_new
为您提供一个与 calc_old
。或者,如果您希望您的第二个实现提供与 calc_old
相同的缩放比例,您应该将 scipy.signal.spectrogram
的结果乘以 0.125 * bf.WINDOW_LEN** 2
。
频点数
给定偶数个点 nperseg
,您的初始实现 calc_old
仅保留 nperseg//2
个频点。
另一方面,完整的非冗余半谱应该为您提供 nperseg//2 + 1
频率仓(有 nperseg-2
具有相应对称性的仓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_LEN
是
bf.WINDOW_LEN
。我猜您的特定输入序列不是这种情况。
相比之下,如果需要,scipy.signal.spectrogram
会用额外的样本填充数据,以便在频谱图计算期间使用所有输入样本,并且与您的相比,这可能会导致一个额外的时间步长calc_old
实现。
关于python - Scipy 频谱图与多个 Numpy FFT,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/51898882/