Python 的意外 FFT 结果

Unexpected FFT Results with Python

我正在分析本质上是呼吸波形的东西,它由 3 种不同的形状构成(数据来自 MRI,因此使用了多个回波时间;如果您想快速了解一些背景知识,请参阅 here) .以下是某些上下文中绘制的两个波形的几个片段:

对于每个波形,我正在尝试执行 DFT 以发现主频率或呼吸频率。

我的问题是,当我绘制我执行的 DFT 时,我得到了两件事之一,这取决于我使用的 FFT 库。此外,它们都不能代表我的期望。我知道数据并不总是像我们想要的那样,但我的数据中显然有波形,所以我会期望离散傅里叶变换在合理的地方产生频率峰值。作为参考,我希望大约 80 到 130 Hz。

我的数据存储在 pandas 数据框中,每个回波时间的数据都在一个单独的系列中。我正在应用选择的 FFT 函数(参见下面的代码)并为每个函数接收不同的结果。

SciPy (fftpack)

import pandas as pd
import scipy.fftpack

# temporary copy to maintain data structure
lead_pts_fft_df = lead_pts_df.copy()

# apply a discrete fast Fourier transform to each data series in the data frame
lead_pts_fft_df.magnitude = lead_pts_df.magnitude.apply(scipy.fftpack.fft)
lead_pts_fft_df

NumPy:

import pandas as pd
import numpy as np 

# temporary copy to maintain data structure
lead_pts_fft_df = lead_pts_df.copy()

# apply a discrete fast Fourier transform to each data series in the data frame
lead_pts_fft_df.magnitude = lead_pts_df.magnitude.apply(np.fft.fft)
lead_pts_fft_df

其余相关代码:

ECHO_TIMES = [0.080, 0.200, 0.400] # milliseconds

f_s = 1 / (0.006) # 0.006 = time between samples
freq = np.linspace(0, 29556, 29556) * (f_s / 29556) # (29556 = length of data)

# generate subplots
fig, axes = plt.subplots(3, 1, figsize=(18, 16))

for idx in range(len(ECHO_TIMES)):
    axes[idx].plot(freq, lead_pts_fft_df.magnitude[idx].real, # real part
                   freq, lead_pts_fft_df.magnitude[idx].imag, # imaginary part

    axes[idx].legend() # apply legend to each set of axes

# show the plot
plt.show()

Post-DFT (SciPy fftpack):

Post-DFT (NumPy)

Here is a link to the dataset (on Dropbox) used to create these plots, and here 是 Jupyter Notebook 的 link。

编辑:

我使用了发布的建议并获取了数据的大小(绝对值),并且还绘制了对数 y 轴。新的结果发布在下面。看来我的信号有一些环绕。我没有使用正确的频率标度吗?更新后的代码和图表如下。

# generate subplots
fig, axes = plt.subplots(3, 1, figsize=(18, 16))

for idx in range(len(ECHO_TIMES)):
    axes[idx].plot(freq[1::], np.log(np.abs(lead_pts_fft_df.magnitude[idx][1::])),
                   label=lead_pts_df.index[idx], # apply labels
                   color=plot_colors[idx]) # colors
    axes[idx].legend() # apply legend to each set of axes

# show the plot
plt.show()

我已经解决了我的问题。

Cris Luengo was very helpful with this link,这帮助我确定了如何正确缩放我的频点并适当地绘制 DFT。

另外:我忘记考虑信号中存在的偏移量。它不仅会导致 DFT 中 0 Hz 处的巨大峰值出现问题,而且还会导致转换信号中的大部分噪声。我利用 scipy.signal.detrend() 消除了这一点,得到了一个非常合适的 DFT。

# import DFT and signal packages from SciPy
import scipy.fftpack
import scipy.signal

# temporary copy to maintain data structure; original data frame is NOT changed due to ".copy()"
lead_pts_fft_df = lead_pts_df.copy()

# apply a discrete fast Fourier transform to each data series in the data frame AND detrend the signal
lead_pts_fft_df.magnitude = lead_pts_fft_df.magnitude.apply(scipy.signal.detrend)
lead_pts_fft_df.magnitude = lead_pts_fft_df.magnitude.apply(np.fft.fft)
lead_pts_fft_df

相应地安排频点:

num_projections = 29556
N = num_projections
T = (6 * 3) / 1000 # 6*3 b/c of the nature of my signal: 1 pt for each waveform collected every third acquisition
xf = np.linspace(0.0, 1.0 / (2.0 * T), num_projections / 2)

然后剧情:

# generate subplots
fig, axes = plt.subplots(3, 1, figsize=(18, 16))

for idx in range(len(ECHO_TIMES)):
    axes[idx].plot(xf, 2.0/num_projections * np.abs(lead_pts_fft_df.magnitude[idx][:num_projections//2]),
                   label=lead_pts_df.index[idx], # apply labels
                   color=plot_colors[idx]) # colors
    axes[idx].legend() # apply legend to each set of axes

# show the plot
plt.show()