完整信号上的相似 FFT 代码或块上的平均值之间的差异

Question

我有 2 段代码应该会产生相同的结果。在第一个中，我考虑在 64 个通道上以 500 Hz 采样的 2500 个样本信号；在第二个中，我将这个信号分成 5 个块，每个块有 500 个样本。

# Snippet 1: 
data = get_data() # data shape is 64x2500

# Snippet 2:
data = get_data().reshape(64, 5, 500) # data shape is 64x5x500

其中 get_data() 只是此示例的占位符，您可以随意生成一个随机样本数组以用于这两段代码。实际上，为了测试我的代码，我用 copy.deepcopy().

复制了一个样本信号

我的目标是通过对频段和所有 64 个通道进行平均来提取 2 个频段的频率功率。

α: (8, 13) 赫兹
增量：(1, 4) 赫兹

为此，我通过 numpy 应用快速傅立叶变换：（下面的代码应该适用于具有两个数组维度的两个片段。）

import numpy as np

alpha = (8, 13)
delta = (1, 4)
fs = 500. # sampling frequency

frequencies = np.fft.rfftfreq(data.shape[-1], 1.0/fs)
alpha_band = np.where(np.logical_and(frequencies>=alpha[0], 
                                     frequencies<=alpha[1]))[0]
delta_band = np.where(np.logical_and(frequencies>=delta[0], 
                                     frequencies<=delta[1]))[0]

fftval = np.abs(np.fft.rfft(data, axis=-1) / fs)

然后对绝对值进行平方。此外，64 个通道中的每一个通道还通过乘法对其施加了权重。为了测试代码，我将所有权重设置为 1。

weights = np.ones(shape=(64,))
alpha = np.average(np.multiply(np.square(fftval[..., alpha_band]).T, weights))
delta = np.average(np.multiply(np.square(fftval[..., delta_band]).T, weights))

我的理解是，对块进行平均不应改变结果。然而，这是我在相同数据上得到的 2 个输出：

# Snippet 1, 64x2500
alpha
Out: 5.294149596657551e-13

delta
Out: 9.372696436349552e-13

alpha/delta
Out: 0.564848081084284

# Snippet 2: 64x5x500
alpha
Out: 6.326672955916193e-12

delta
Out: 7.584706602278469e-13

alpha/delta
Out: 8.34135489699185

有谁知道我在这里做错了什么以及为什么两个结果完全不同？我错过了什么？

Answer 1

只有为您采样频率和分块的“完美”信号，分块频谱和全信号频谱的平均值才会相同。

您应该对信号应用 windowing（例如 np.hamming），特别是如果您正在寻找如此低频的组件。

fft（dft）是循环的，起点和终点相遇，如果有低频分量不重复多次（甚至一次），起点和终点连接时会不连续。这种影响称为频谱泄漏。它“创造”了干扰整体结果的虚假高频成分。 window 功能的选择取决于您的需要。所以不是完全直接但测试一些，看看它是否改进！

windowing https://download.ni.com/evaluation/pxi/Understanding%20FFTs%20and%20Windowing.pdf

上有很好的资源

重复测试，但在片段 1 中对整个信号执行 windowing。在片段 2 中，对每个块执行 windowing。可能即使现在结果也不完全相同，但希望更接近！

Answer 2

我认为问题出在这里：

fftval[..., alpha_band]

对于代码段 2 没问题，fft 值与频率相关 [8., 9., 10., 11., 12., 13.]
但是对于片段 1，您的 alpha 波段的频率是 [1.6, 1.8, 2. , 2.2, 2.4, 2.6] Hz
还有

fftval = np.abs(np.fft.rfft(data, axis=-1) / data.shape[-1])

完整信号上的相似 FFT 代码或块上的平均值之间的差异

Difference between similar FFT code on complete signal or averaged on chunks

python

numpy

fft