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.
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.
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
nperseg = (2 / low) * sf
# Compute the modified periodogram (Welch)
freqs, psd = welch(data, sf, nperseg=nperseg)
# Frequency resolution
freq_res = freqs - freqs
# 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)
bp /= simps(psd, dx=freq_res)
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)
And here is how I called it in the rolling window function:
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.