python - 在 Python 中以 Hz 为单位绘制 FFT 频率
问题描述
我正在实施本文中的方法: https ://dspace.mit.edu/bitstream/handle/1721.1/66243/Picard_Noncontact%20Automated.pdf?sequence=1&isAllowed=y
主要思想是使用来自 10 秒视频的一组帧 (N=300) 进行心脏脉搏测量,因此帧速率等于 30 fps。
red = [item[:,:,0] for item in imgs]
green = [item[:,:,1] for item in imgs]
blue = [item[:,:,2] for item in imgs]
red_avg = [item.mean() for item in red]
green_avg = [item.mean() for item in green]
blue_avg = [item.mean() for item in blue]
red_mean, red_std = np.array(red_avg).mean(), np.array(red_avg).std()
green_mean, green_std = np.array(green_avg).mean(), np.array(green_avg).std()
blue_mean, blue_std = np.array(blue_avg).mean(), np.array(blue_avg).std()
red_avg = [(item - red_mean)/red_std for item in red_avg]
green_avg = [(item - green_mean)/green_std for item in green_avg]
blue_avg = [(item - blue_mean)/blue_std for item in blue_avg]
data = np.vstack([signal.detrend(red_avg), signal.detrend(green_avg), signal.detrend(blue_avg)]).reshape(300,3)
from sklearn.decomposition import FastICA
transformer = FastICA(n_components=3)
X_transformed = transformer.fit_transform(data)
from scipy.fftpack import fft
first = X_transformed.T[0]
second = X_transformed.T[1]
third = X_transformed.T[2]
ff = np.fft.fft(first)
fs = np.fft.fft(second)
ft = np.fft.fft(third)
imgs
- 是具有 300 个图像像素值的数组的初始列表。如您所见,我将所有帧拆分为 RGB 通道,因此有痕迹x_i(t)
,其中i = 1,2,3
标准化后,我们去除所有轨迹并将它们堆叠以进一步应用 ICA,然后对所有三个分量进行 FFT。
然后该方法声称我们需要绘制功率与频率 (Hz) 的关系图,并选择最有可能是心脏脉冲的分量。
最后,我们对选定的源信号应用快速傅里叶变换 (FFT) 以获得功率谱。脉冲频率被指定为对应于操作频带内频谱的最高功率的频率。对于我们的实验,我们将操作范围设置为 [0.75, 4] Hz(对应于 [45, 240] bpm)以提供广泛的心率测量。
这是我尝试可视化频率的方法:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
data = ft
print(fs.size)
ps = np.abs(np.fft.fft(data))**2
sampling_rate = 30
freqs = np.fft.fftfreq(data.size, 1/sampling_rate)
idx = np.argsort(freqs)
#print(idx)
plt.plot(freqs[idx], ps[idx])
我得到的完全不同,因为频率范围从 -15 美元到 15 美元,我不知道这是否以赫兹为单位。
上面的三个图像是我在执行代码以可视化频率和信号功率时得到的。
我将不胜感激任何帮助或建议。
解决方案
您应该真正学习如何将图像/视频用作 nD-张量。这样做,您可以用更简洁的代码替换所有数据争吵:
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft
from sklearn.decomposition import FastICA
images = [np.random.rand(640, 480, 3) for _ in range(30)]
# Create tensor with all images
images = np.array(images)
images.shape
# Take average of all pixels, for each image and each channel individually
avgs = np.mean(images, axis=(1, 2))
mean, std = np.mean(avgs), np.std(avgs)
# Normalize all average channels
avgs = (avgs - mean) / std
# Detrend across images
avgs = scipy.signal.detrend(avgs, axis=0)
transformer = FastICA(n_components=3)
X_transformed = transformer.fit_transform(avgs)
X_ff = np.fft.fft(X_transformed, axis=0)
plt.plot(np.abs(X_ff) ** 2)
稍微回答一下您的问题:我认为您错误地采用了第三个 PCA 组件的傅立叶光谱的傅立叶光谱。
FFT(FFT(PCA[:, 2]))
虽然您打算只进行一次 FFT:
FFT(PCA[:, 2]))
关于 -15...15 轴:您已将采样频率设置为 30Hz(或视频术语中的 30fps)。这意味着您可以在视频中检测到高达 15Hz 的任何内容。
在傅里叶理论中,存在一种叫做“负频率”的东西。现在,由于我们主要分析真实信号(而不是复杂信号),因此负频率始终与正频率相同。这意味着您的频谱始终是对称的,您可以忽略左半部分。
但是,由于您已经进行了两次 FFT,因此您正在查看复杂信号的 FFT,该信号确实具有负频率。这就是为什么您的光谱不对称且令人困惑的原因。
此外,我相信您会混淆重塑和转置。在 PCA 之前,您正在组装数据,例如
np.vstack([red, green, blue]) # shape = (3, 300)
你想转置得到(300, 3)
. 如果您改为重塑,则不会交换行和列,而是以不同的形状解释相同的数据。
推荐阅读
- reactjs - 是否可以从 HTML 默认播放器中删除控件(
- c++ - 如果参数与数据成员的名称相同怎么办?
- java - Java 流减少功能。累加器的关联性
- python - 如何将python中的多行返回到html?
- antd - 如何在 antdTable 中使用“defaultSortOrder”?
- sql-server - WAS 6.1 的 SQL Server 2016 问题
- path - 如何使用相对路径为 MSTest 的 OpenCover 编写批处理文件,使其变得独立于机器?
- php - 如何在 foreach 循环中创建数组 - php?
- excel - 使用用户窗体调整列表框大小
- sql - 根据年份过滤日期