如何仅对信号中的较低频率进行高分辨率 FFT?

信息处理 fft 傅里叶变换 Python 频谱图
2022-02-08 23:09:43

数据及问题描述:

我有一个以 1000 Hz 采样的信号。我在 120 Hz 处对其进行低通滤波,并希望制作低于此阈值的频率的频谱图。我正在使用 scipy 函数fftpack.fftfreq 来获取傅立叶系数,然后fftpack.fft 进行实际变换。

信号很长,我希望使用 50 毫秒窗口的频谱图长度约为 5 秒。我还过滤掉负频率进行绘图。给定此窗口中的样本数(50),我得到 50 个傅立叶系数。但是,其中一半是负数,另一部分介于 120 Hz 和 500 Hz 之间(自然,因为采样率为 1000 Hz)。

这给我留下了相当低分辨率的频谱图。我在感兴趣的频率区域(0、20、40、60、100 和 120 Hz)中只有大约 6 个块,然后有 18 个块在高达 500 Hz 的频率中显示低活动,因为它们被过滤掉了。

问题:

例如,我怎么能在 FFT 中指定 0 到 120 Hz 之间的 50 个频率区间?

我尝试了什么:

np.linspace我尝试通过使用来指定频率(而不是)来做这样的事情fft.fftfreq,但这引入了一些其他错误,即频谱图总是看起来像一面镜子,在最高和最低频率处具有高功率,在中间具有低功率图表,与范围无关。老实说,我不确定为什么 linspace 而不是 fftfreq 会发生这种情况。它们都返回浮点数组。

代码示例包括在下面。提前感谢您的帮助!干杯

X = samples_filt  # long 1-D vector of low-pass filtered samples
fs = 1000
time = .05  # window length in sec
N = int(fs * time)  # num samples

tot_len = 5 * fs  # 5 sec of the whole signal
X = X[:tot_len]

f = fftpack.fftfreq(N, 1.0/fs)
# f = np.ceil(np.linspace(0, max_freq, 52)[1:51])  # introduced error described above
mask = (f > 0)  # mask for positive freqs

n_max = int(np.ceil(X.shape[0] / N))  # the number of segments of length N in the sample array data
f_values = np.sum(1 * mask)  # how many values meet mask reqs
spectogram_data = np.zeros((f_values, n_max))

window = sp.signal.blackman(N)  # taper used to improve contrast of spectrogram

for n in range(0, n_max):
    subdata = X[(N * n):(N * (n + 1))]
    F = fftpack.fft(subdata * window)
    spectogram_data[:, n] = np.log(abs(F[mask]))

fig, ax = plt.subplots(1, 1, figsize=(8, 6))
p = ax.imshow(spectogram_data, origin='lower', 
    extent=(0, X.shape[0] / fs, 0, max(f)), 
    aspect='auto',
    cmap=mpl.cm.RdBu_r)
cb = fig.colorbar(p, ax=ax)

样本输出:

在此处输入图像描述

例如,如果有相同的图表,其中包含 24 个频率区间,最高可达 120 Hz,那就太好了。

3个回答

对于给定的窗口大小,您将获得多少线性频率分辨率有一个基本限制,至少使用传统的线性技术并且不对源信号做任何假设。

有些技术与 FFT 几乎相同,但只返回频率区间的子集。它们通常用于成本原因。

为方便起见,您可能会考虑将信号重新采样为 2x120=240Hz,然后根据您需要的时间分辨率进行超过 50ms 窗口的频谱图。

频率分辨率为:d_f = fs/L

对于 d_s = 120/50=2.4 Hz,您可能需要以下内容:

L = fs/d_f=240/2.4=100 个样本

即窗口长度:100/240 = 0.42 秒

例如,我怎么能在 FFT 中指定 0 到 120 Hz 之间的 50 个频率区间?

如果您想要的分辨率,那么您需要的样本长度至少使用窗口,您将被限制为Δf=120Hz50Δt=1Δf0.416¯sΔt=50msΔf=150ms=20Hz

以 1kHz 录制并将 3/4 的结果用于过滤实际上并没有那么糟糕。你可以搞乱预过滤,但你最终可能会发现这样做涉及 FFT、IFFT 和重叠相加——所以你不妨去掉中间人,对“整体”进行 FFT “ 数据。

如果您担心从 FFT 中获取所有额外的东西,那么 FFT 算法的一些风格旨在获取真实(非复杂)样本的向量并仅返回 FFT 的正频率侧。它们运行得更快一些,因为对于真实的输入数据,可以利用一些对称性。

任何体面的 FFT 包(即 FFTW)都有它们——你只需要在文档中挖掘一下就可以找到它们。

这可以通过 STFT 的行实现来实现,而不是标准的列实现。它目前是 ssqueezepy 中的 TODO,有一个已经工作的变体断言与标准相等。

但是请注意,只是更多的频率箱并不一定会产生更高的频率分辨率。与此答案相同的推理,但沿频率轴-加上一个从根本上由普通 DFT 的分辨率强加的。