How to extract the Freezing Index from Accelerometer data for FoG detecting?

Question

I'm working on a project to detect Freezing of Gait episodes for Parkinson's Disease patients, and according to this paper and other ones, extracting the freezing index which is "The power in the “freeze” band (3- 8Hz) divided by the power in the locomotor band (0.5-3Hz)" will be a good feature.

Here is the explanation: One standard feature which is extracted from the raw signals is the Freezing Index (FI), defined as the ratio between the power contained in the so-called freezing and locomotion frequency bands (3-8 Hz and 0.5-3 Hz respectively). This feature is convenient since it requires only FFT-computation.

But I cannot how to implement it in Python.

I have a dataframe like this: Then I did something like this to extract features from the sensors time-series data:

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]

Now, I would like to extract the Freezing Index, How can I do this? and I need to someone to explain it to me because I don't fully understand it!

Answer 1

I found this article is very useful. There is a function in this article that computes the average power of the signal in a specific frequency band. And since the Freezing Index is the ratio between the power contained in the so-called freezing and locomotion frequency bands (3-8 Hz and 0.5-3 Hz respectively), we can use this function to get the power in each frequency band and divide them.

Here is the function:

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

Then, I created this simple function to return the 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

And here is how I called it in the rolling window function:

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]

I hope it's the right solution, and I wish it helps.