如何从用于 FoG 检测的加速度计数据中提取冻结指数?
How to extract the Freezing Index from Accelerometer data for FoG detecting?
我正在开展一个项目来检测帕金森病患者的步态冻结发作,并根据这个 paper 和其他人,提取冻结指数,即“功率在“冻结”带(3-
8Hz) 除以运动频带中的功率 (0.5-3Hz)" 将是一个很好的特征。
解释如下:从原始信号中提取的一个标准特征是冻结指数 (FI),定义为所谓的冻结中包含的功率与运动频率之间的比率频段(分别为 3-8 Hz 和 0.5-3 Hz)。此功能很方便,因为它只需要 FFT 计算。
但是我不知道如何在Python中实现它。
我有一个这样的数据框:
然后我做了这样的事情来从传感器时间序列数据中提取特征:
win_size=200
step_size=40
for col in df.columns:
if col != 'Label':
df[col + '_avg'] = df[col].rolling(win_size).mean()[step_size - 1::step_size]
现在,我想提取冷冻指数,我该怎么做?我需要有人向我解释,因为我不完全理解!
我发现这个 article 非常有用。这篇文章中有一个函数可以计算特定频段内信号的平均功率。并且由于冷冻指数是所谓的冷冻和运动频带(分别为 3-8 Hz 和 0.5-3 Hz)中包含的功率之间的比率,我们可以使用此函数来获取每个频带中的功率并除以他们。
函数如下:
def bandpower(data, sf, band, window_sec=None, relative=False):
"""Compute the average power of the signal x in a specific frequency band.
Parameters
----------
data : 1d-array
Input signal in the time-domain.
sf : float
Sampling frequency of the data.
band : list
Lower and upper frequencies of the band of interest.
window_sec : float
Length of each window in seconds.
If None, window_sec = (1 / min(band)) * 2
relative : boolean
If True, return the relative power (= divided by the total power of the signal).
If False (default), return the absolute power.
Return
------
bp : float
Absolute or relative band power.
"""
from scipy.signal import welch
from scipy.integrate import simps
band = np.asarray(band)
low, high = band
# Define window length
if window_sec is not None:
nperseg = window_sec * sf
else:
nperseg = (2 / low) * sf
# Compute the modified periodogram (Welch)
freqs, psd = welch(data, sf, nperseg=nperseg)
# Frequency resolution
freq_res = freqs[1] - freqs[0]
# Find closest indices of band in frequency vector
idx_band = np.logical_and(freqs >= low, freqs <= high)
# Integral approximation of the spectrum using Simpson's rule.
bp = simps(psd[idx_band], dx=freq_res)
if relative:
bp /= simps(psd, dx=freq_res)
return bp
然后,我为 return FI 创建了这个简单的函数:
def freeze_index(data, sf, band1, band2, win_sec=None, relative=False):
fi = bandpower(data, sf, band1, win_sec, relative) / bandpower(data, sf, band2, win_sec, relative)
return fi
下面是我在滚动 window 函数中调用它的方式:
for col in df.columns:
if col != 'Label':
df[col + '_fi'] = df[col].rolling(win_size).apply(freeze_index, args=(fs, [3, 8], [0.5, 3], win_size,))[step_size - 1::step_size]
我希望这是正确的解决方案,希望对您有所帮助。
我正在开展一个项目来检测帕金森病患者的步态冻结发作,并根据这个 paper 和其他人,提取冻结指数,即“功率在“冻结”带(3- 8Hz) 除以运动频带中的功率 (0.5-3Hz)" 将是一个很好的特征。
解释如下:从原始信号中提取的一个标准特征是冻结指数 (FI),定义为所谓的冻结中包含的功率与运动频率之间的比率频段(分别为 3-8 Hz 和 0.5-3 Hz)。此功能很方便,因为它只需要 FFT 计算。
但是我不知道如何在Python中实现它。
我有一个这样的数据框:
win_size=200
step_size=40
for col in df.columns:
if col != 'Label':
df[col + '_avg'] = df[col].rolling(win_size).mean()[step_size - 1::step_size]
现在,我想提取冷冻指数,我该怎么做?我需要有人向我解释,因为我不完全理解!
我发现这个 article 非常有用。这篇文章中有一个函数可以计算特定频段内信号的平均功率。并且由于冷冻指数是所谓的冷冻和运动频带(分别为 3-8 Hz 和 0.5-3 Hz)中包含的功率之间的比率,我们可以使用此函数来获取每个频带中的功率并除以他们。
函数如下:
def bandpower(data, sf, band, window_sec=None, relative=False):
"""Compute the average power of the signal x in a specific frequency band.
Parameters
----------
data : 1d-array
Input signal in the time-domain.
sf : float
Sampling frequency of the data.
band : list
Lower and upper frequencies of the band of interest.
window_sec : float
Length of each window in seconds.
If None, window_sec = (1 / min(band)) * 2
relative : boolean
If True, return the relative power (= divided by the total power of the signal).
If False (default), return the absolute power.
Return
------
bp : float
Absolute or relative band power.
"""
from scipy.signal import welch
from scipy.integrate import simps
band = np.asarray(band)
low, high = band
# Define window length
if window_sec is not None:
nperseg = window_sec * sf
else:
nperseg = (2 / low) * sf
# Compute the modified periodogram (Welch)
freqs, psd = welch(data, sf, nperseg=nperseg)
# Frequency resolution
freq_res = freqs[1] - freqs[0]
# Find closest indices of band in frequency vector
idx_band = np.logical_and(freqs >= low, freqs <= high)
# Integral approximation of the spectrum using Simpson's rule.
bp = simps(psd[idx_band], dx=freq_res)
if relative:
bp /= simps(psd, dx=freq_res)
return bp
然后,我为 return FI 创建了这个简单的函数:
def freeze_index(data, sf, band1, band2, win_sec=None, relative=False):
fi = bandpower(data, sf, band1, win_sec, relative) / bandpower(data, sf, band2, win_sec, relative)
return fi
下面是我在滚动 window 函数中调用它的方式:
for col in df.columns:
if col != 'Label':
df[col + '_fi'] = df[col].rolling(win_size).apply(freeze_index, args=(fs, [3, 8], [0.5, 3], win_size,))[step_size - 1::step_size]
我希望这是正确的解决方案,希望对您有所帮助。