如何使用 Python 对原始信号应用 FFT
How to apply FFT on raw signal using Python
我是处理信号的新手,需要你的帮助。
我从我的 TI AFE4490 获得了一个 10 秒的原始 PPG(光体积描记图)信号。我的硬件已校准,我每秒使用 250 个样本来记录这些信号。我最后获得了2500点
你可以看到下面的图像,点和代码。
顶部:我的原始 PPG 信号 - 底部:尝试应用 FFT:
代码:
RED, IR, nSamples, sRate = getAFESignal()
period = 1/sRate
plt.figure(1)
plt.subplot(2,1,1)
x = np.linspace(0.0, nSamples*period, nSamples)
y = IR
plt.xlabel("Time (s)")
plt.ylabel("Voltage (V)")
plt.plot(x,y)
plt.subplot(2,1,2)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*period), nSamples//2)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Gain")
plt.plot(xf, 2.0/nSamples * np.abs(yf[0:nSamples//2]))
plt.grid()
plt.show()
函数getAFEsignal()
只是一个读取.txt文件并将所有内容放入两个numpy数组的函数。
您可以在此处找到 .txt 文件:Raw signal file
如您所见,我没有正确应用 FFT,我需要这个来发现我需要过滤的频率。你知道我做错了什么吗?是否可以对这个信号应用 FFT?
好消息是您的 FFT 计算没有问题。您在时域中显示的数据具有相当强的低频分量。相应地,您获得的频域图在 0Hz 附近显示出明显的尖峰。
主要问题在于如何绘制结果。为了根据对时域波形的直观感知更好地查看您可能希望在频域中看到的内容,您需要重新调整每个轴的比例。特别是,在显示的时间尺度上,您可能会注意到持续时间约为 0.25 秒的模式
最多可能几秒钟。这对应于大约 0-5Hz 的频率范围。关注该范围而不是显示整个 0-125Hz 频谱是有意义的。这可以通过这样设置 x 轴限制来实现:
plt.xlim(0,5) # set x-axis limits from 0 to 5Hz
与 y 轴类似,您需要考虑具有小振幅的频率分量(在线性标度上变得更难注意到)仍然可以对时域信号产生非常明显的贡献。因此,通常需要以对数分贝标度显示频域振幅。这可以按如下方式完成:
plt.plot(xf, 20*np.log10(2.0/nSamples * np.abs(yf[0:nSamples//2])))
最后,如果您想更好地了解某些特定频率分量的贡献而不受其他频率分量 spectral leakage 的干扰,您可能需要考虑在计算快速傅立叶变换。例如,如果你想消除恒定信号偏差的影响,缓慢的 ~0.1Hz 变化和频率大于 10Hz 的噪声,你可以使用如下内容:
import scipy.signal
b,a = signal.butter(4, [0.25/sRate, 10/sRate], btype='bandpass')
y = signal.filtfilt(b,a,signal.detrend(IR, type='constant'))
我是处理信号的新手,需要你的帮助。
我从我的 TI AFE4490 获得了一个 10 秒的原始 PPG(光体积描记图)信号。我的硬件已校准,我每秒使用 250 个样本来记录这些信号。我最后获得了2500点
你可以看到下面的图像,点和代码。
顶部:我的原始 PPG 信号 - 底部:尝试应用 FFT:
代码:
RED, IR, nSamples, sRate = getAFESignal()
period = 1/sRate
plt.figure(1)
plt.subplot(2,1,1)
x = np.linspace(0.0, nSamples*period, nSamples)
y = IR
plt.xlabel("Time (s)")
plt.ylabel("Voltage (V)")
plt.plot(x,y)
plt.subplot(2,1,2)
yf = fft(y)
xf = np.linspace(0.0, 1.0/(2.0*period), nSamples//2)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Gain")
plt.plot(xf, 2.0/nSamples * np.abs(yf[0:nSamples//2]))
plt.grid()
plt.show()
函数getAFEsignal()
只是一个读取.txt文件并将所有内容放入两个numpy数组的函数。
您可以在此处找到 .txt 文件:Raw signal file
如您所见,我没有正确应用 FFT,我需要这个来发现我需要过滤的频率。你知道我做错了什么吗?是否可以对这个信号应用 FFT?
好消息是您的 FFT 计算没有问题。您在时域中显示的数据具有相当强的低频分量。相应地,您获得的频域图在 0Hz 附近显示出明显的尖峰。
主要问题在于如何绘制结果。为了根据对时域波形的直观感知更好地查看您可能希望在频域中看到的内容,您需要重新调整每个轴的比例。特别是,在显示的时间尺度上,您可能会注意到持续时间约为 0.25 秒的模式 最多可能几秒钟。这对应于大约 0-5Hz 的频率范围。关注该范围而不是显示整个 0-125Hz 频谱是有意义的。这可以通过这样设置 x 轴限制来实现:
plt.xlim(0,5) # set x-axis limits from 0 to 5Hz
与 y 轴类似,您需要考虑具有小振幅的频率分量(在线性标度上变得更难注意到)仍然可以对时域信号产生非常明显的贡献。因此,通常需要以对数分贝标度显示频域振幅。这可以按如下方式完成:
plt.plot(xf, 20*np.log10(2.0/nSamples * np.abs(yf[0:nSamples//2])))
最后,如果您想更好地了解某些特定频率分量的贡献而不受其他频率分量 spectral leakage 的干扰,您可能需要考虑在计算快速傅立叶变换。例如,如果你想消除恒定信号偏差的影响,缓慢的 ~0.1Hz 变化和频率大于 10Hz 的噪声,你可以使用如下内容:
import scipy.signal
b,a = signal.butter(4, [0.25/sRate, 10/sRate], btype='bandpass')
y = signal.filtfilt(b,a,signal.detrend(IR, type='constant'))