mne中计算Hjorth的函数

MNE中没有直接计算Hjorth函数的函数，但可以通过以下步骤计算：

1. 导入需要使用的库和数据

import numpy as np
from mne import io, Epochs
from mne.time_frequency import psd_array_multitaper

2. 读取数据并创建Epochs

raw = io.read_raw_edf('data.edf', preload=True)
events = mne.find_events(raw)
event_id = {'event': 1}
tmin, tmax = -1, 1
epochs = Epochs(raw, events, event_id, tmin, tmax, baseline=None)

3. 计算每个通道的功率谱密度

freqs, psds = psd_array_multitaper(epochs.get_data(), sfreq=raw.info['sfreq'], fmin=0, fmax=50)

4. 计算每个通道的Hjorth参数

def hjorth(x):
    """
    Compute Hjorth mobility and complexity of a time series from either two
    cases:
        1.  x, the time series of type list (length N)
        2.  x, the power spectral density (PSD) of a time series of type numpy
            array (length N/2+1)

    In case 1, the function returns the tuple (m, c), where m is the mobility
    and c is the complexity.

    In case 2, the function returns the tuple (m, c, freq), where freq is the
    frequency vector of length N/2+1.

    In both cases, the units of the output values depend on the input values.
    In case 1, if x has unit u, m and c will have units of u/Hz and u/Hz^2,
    respectively. In case 2, if x has unit u^2/Hz, m and c will have units of
    1/Hz and u^2/Hz^3, respectively.

    This implementation is based on the algorithm described in the following
    paper:

    Hjorth, B. (1970). EEG analysis based on time domain properties.
    Electroencephalography and Clinical Neurophysiology, 29(3), 306-310.

    Parameters
    ----------
    x : list | np.ndarray
        The time series or PSD.

    Returns
    -------
    m : float
        The Hjorth mobility.
    c : float
        The Hjorth complexity.
    freq : np.ndarray | None
        The frequency vector if input is a PSD, otherwise None.
    """
    if isinstance(x, list):
        dx = np.diff(x)
        ddx = np.diff(dx)
        m2 = np.mean(dx**2)
        m4 = np.mean(ddx**2)
        return np.sqrt(m2 / x.var()), np.sqrt(m4 / m2) / x.var()
    else:
        w = 2 * np.pi * np.fft.fftfreq(len(x), d=1. / (2 * (len(x) - 1)))
        w = w[:len(x) // 2 + 1]
        pw = np.abs(x[:len(x) // 2 + 1])**2
        m2w = np.trapz(pw * w**2, w)
        m0w = np.trapz(pw, w)
        m4w = np.trapz(pw * w**4, w)
        return np.sqrt(m2w / m0w), np.sqrt(m4w / m2w) / m0w, w

hjorth_params = []
for psd in psds:
    hjorth_params.append(hjorth(psd)[0]) # 计算Hjorth mobility

5. 将Hjorth参数添加到原始数据中

for i, ch_name in enumerate(raw.info['ch_names']):
    raw.info['chs'][i]['hjorth_mobility'] = hjorth_params[i]