如何在Python中计算Dimitrov频谱疲劳指数及解决inf值问题

问题描述 #

我有一段计算肌电信号平均功率频率的Python脚本，现在需要计算Dimitrov频谱疲劳指数。该指数公式和平均功率频率不同，采用-1阶矩与5阶矩的比值，而非1阶矩与0阶矩的比值。我将原函数最后一行代码修改为mean_freq[i] = np.dot(P,f.T**-1)/np.dot(P,f.T**5)，并重命名为FInms5函数后，得到的结果全是inf值，求解决办法。

原工作函数代码如下：

from scipy.signal import periodogram

def get_mean_freq(emg_sig, sfreq, epoch_duration = 0.5):
'''
Parameters
----------
emg_sig : array
pre-filtered emg data.
sfreq : int
emg sampling frequency, in Hz.
epoch_duration : float
epoch (time window) duration, in seconds.

Returns
-------
mean_freq: array
    mean frequency at each epoch
time_points: array
    time point at the center of each evaluated epoch
samples: array
    sample numbers at the center of each evaluated epoch

Method according to:
    https://stackoverflow.com/questions/37922928/difference-in-mean-frequency-in-python-and-matlab
'''
ons = range(
    0,
    len(emg_sig),
    int(epoch_duration*sfreq)
    )

mean_freq = np.empty((len(ons),))
samples =np.empty((len(ons),))
time_points = np.empty((len(ons),))
for i,on in enumerate(ons):# i,on = 0,ons[0]
    off = ons[i+1]-1 if i+1<len(ons) else len(emg_sig)
    processing_window = emg_sig[on:off]
    mid_point = (on + off) / 2
    samples[i] = mid_point
    time_points[i] = mid_point / sfreq
    f, Pxx_den = periodogram(np.array(processing_window), fs=float(sfreq))
    Pxx_den = np.reshape(Pxx_den, (1,-1))
    width = np.tile(f[2]-f[0], (1, Pxx_den.shape[1]))
    f = np.reshape(f, (1,-1))
    P = Pxx_den*width
    pwr = np.sum(P)
    mean_freq[i] = np.dot(P, f.T)/pwr

return mean_freq, time_points, samples

原因分析 #

出现inf值的核心原因是periodogram返回的频率数组f第一个元素是0Hz，计算f**-1时，0的-1次方直接得到无穷大，导致分子或整体计算结果变为inf。另外，肌电信号的直流分量（0Hz）通常没有生理意义，应该被排除。

解决方案 #

• 从频率数组f和功率数组Pxx_den中移除0Hz对应的分量，即从索引1开始截取数组；
• 调整矩的计算逻辑，使用处理后的数组进行点积运算。

修改后的函数代码：

from scipy.signal import periodogram
import numpy as np

def get_dimitrov_fi(emg_sig, sfreq, epoch_duration = 0.5):
    '''
    Parameters
    ----------
    emg_sig : array
        预处理后的肌电数据.
    sfreq : int
        肌电信号采样频率，单位Hz.
    epoch_duration : float
        分析时间窗长度，单位秒.

    Returns
    -------
    dimitrov_fi: array
        每个时间窗的Dimitrov频谱疲劳指数
    time_points: array
        每个分析时间窗的中心时间点
    samples: array
        每个分析时间窗的中心采样点编号
    '''
    ons = range(
        0,
        len(emg_sig),
        int(epoch_duration*sfreq)
        )

    dimitrov_fi = np.empty((len(ons),))
    samples =np.empty((len(ons),))
    time_points = np.empty((len(ons),))
    for i,on in enumerate(ons):
        off = ons[i+1]-1 if i+1<len(ons) else len(emg_sig)
        processing_window = emg_sig[on:off]
        mid_point = (on + off) / 2
        samples[i] = mid_point
        time_points[i] = mid_point / sfreq
        f, Pxx_den = periodogram(np.array(processing_window), fs=float(sfreq))
        
        # 移除0Hz分量
        f = f[1:]
        Pxx_den = Pxx_den[1:]
        
        Pxx_den = np.reshape(Pxx_den, (1,-1))
        width = np.tile(f[1]-f[0], (1, Pxx_den.shape[1]))
        f = np.reshape(f, (1,-1))
        P = Pxx_den*width
        
        # 计算-1阶矩和5阶矩的比值
        numerator = np.dot(P, f.T ** -1)
        denominator = np.dot(P, f.T ** 5)
        # 添加微小偏移量避免极端情况的除以0
        dimitrov_fi[i] = numerator / (denominator + 1e-10)

    return dimitrov_fi, time_points, samples

额外注意事项 #

• 确保输入的emg_sig已经完成预处理（如去除直流分量、带通滤波），避免残留的直流分量影响结果；
• 代码中添加了1e-10的偏移量，可避免极端情况下分母为0的问题。

内容的提问来源于stack exchange，提问作者Marcel Vogelsang