关于肌电脉冲编码的一些代码

一、基于阈值的表示法

def TBR_encod(data_sample, sampling_rate, f=1, refractory=5):
    '''
    Input:
        - EMG data sample, sampling rate of EMG electrodes, 'f' factor for threshold tuning, refractory period in ms
    Output:
        - arrays of UP and DOWN spike times (size=dependt. on no. of spikes)
        - arrays of spike indeces w.r.t. to the data_sample array (i.e., at what sample point does spike occur)
        - spike data (size=len(data_sample, values ={-1, 0, 1})
        - threshold (for signal reconstruction)
    Comment:
        - converts EMG data into spike trains using the temporal contrast method (also called threshold-based representation)
    '''
    refractory_samples = round((refractory/1000)*sampling_rate) #number of samples corresponding to a refractory period
    # plot_EMG(data_sample, no_electrodes=2, pose='1')
    #save delta values between consecutive samples
    diff = np.zeros(len(data_sample))
    #determining threshold (V_th) based on mean and std of data
    for t in range(len(data_sample)-1):
        diff[t] = data_sample[t+1]-data_sample[t]
    V_th = np.mean(diff)+f*np.std(diff) #negative threshold [uV]

    UP_spikes = []
    DOWN_spikes = []
    spikes = np.zeros(len(data_sample)) #both UP and DOWN
    it = iter(range(len(data_sample)-1)) #set up an iterator
    #extracting spikes based on threshold
    for t in it:
        delta = data_sample[t+1]-data_sample[t]
        if delta > V_th:
            UP_spikes.append(t)
            spikes[t] = 1
            #enforce refractory period by skipping over samples
            for _ in range(refractory_samples):
                try:
                    t = next(it)
                except: #to prevent code from stopping when t reaches end of iterator
                    break
        if delta < -V_th:
            DOWN_spikes.append(t)
            spikes[t] = -1
            #enforce refractory period by skipping over samples
            for _ in range(refractory_samples):
                try:
                    t = next(it)
                except: #to prevent code from stopping when t reaches end of iterator
                    break
    
    #converting from sample indeces to spike times
    UP_spike_times = np.array(UP_spikes)/sampling_rate
    DOWN_spike_times = np.array(DOWN_spikes)/sampling_rate
    
    
    return UP_spike_times, DOWN_spike_times, UP_spikes, DOWN_spikes, spikes, V_th
'''
这个函数 TBR_encod 实现了一种称为“基于阈值的表示法”（Threshold-Based Representation, TBR）的编码方法，用于将肌电图（EMG）数据转换成尖峰时间序列。这种编码方法通过检测信号的变化（delta）来确定尖峰，并使用一定的阈值和不应期来生成上尖峰（UP spike）和下尖峰（DOWN spike）。
具体的步骤如下：

计算不应期样本数：将不应期（以毫秒为单位）转换为对应的采样数。
计算差值：计算相邻样本之间的差值。
确定阈值：基于差值的均值和标准差确定阈值（V_th）。
提取尖峰：
遍历数据样本，如果差值超过阈值，则记录为上尖峰，并跳过不应期内的样本。
如果差值低于负阈值，则记录为下尖峰，并跳过不应期内的样本。
转换尖峰时间：将尖峰的样本索引转换为时间。
函数参数和返回值说明
输入参数：

data_sample: EMG数据样本。
sampling_rate: EMG电极的采样率。
f: 调整阈值的因子。
refractory: 不应期（毫秒）。
输出：
UP_spike_times: 上尖峰的时间数组。
DOWN_spike_times: 下尖峰的时间数组。
UP_spikes: 上尖峰的样本索引数组。
DOWN_spikes: 下尖峰的样本索引数组。
spikes: 尖峰数据数组（大小与数据样本相同，值为{-1, 0, 1}）。
V_th: 用于信号重建的阈值。
主要编码思路
使用差值和阈值检测信号的显著变化。
记录尖峰的位置并跳过不应期内的样本以避免重复检测。
将尖峰的样本索引转换为时间，以便后续处理和分析。
这个函数主要用于分析肌电信号，通过将连续的EMG数据转换为离散的尖峰序列，可以更有效地捕捉肌电活动的瞬时变化。
'''

二、步进编码

def SF_encod(data_sample, sampling_rate, f=1, refractory=5):
    '''
    Input:
        - EMG data sample, sampling rate of EMG electrodes, 'f' factor for threshold tuning, refractory period in ms
    Output:
        - arrays of UP and DOWN spike times (size=dependt. on no. of spikes)
        - arrays of spike indeces w.r.t. to the data_sample array (i.e., at what sample point does spike occur)
        - spike data (size=len(data_sample, values ={-1, 0, 1})
        - threshold (for signal reconstruction)
    Comment:
        - converts EMG data into spike trains using the Step-Forward (SF) Encoding
    '''
    refractory_samples = round((refractory/1000)*sampling_rate) #number of samples corresponding to a refractory period 计算不应期样本数
    # plot_EMG(data_sample, no_electrodes=2, pose='1')
    #save delta values between consecutive samples
    diff = np.zeros(len(data_sample))
    #determining threshold (V_th) based on mean and std of data
    for t in range(len(data_sample)-1):
        diff[t] = data_sample[t+1]-data_sample[t]
    V_th = np.mean(diff)+f*np.std(diff) #negative threshold [uV]
    
    initial = data_sample[0]
    spikes = np.zeros(len(data_sample))
    baseline = initial
    UP_spikes = []
    DOWN_spikes = []
    it = iter(range(1,len(data_sample)-1)) #set up an iterator
    #extracting spikes based on threshold
    for t in it:
        if data_sample[t] > baseline + V_th:
            spikes[t] = 1
            UP_spikes.append(t)
            baseline += V_th
            #enforce refractory period by skipping over samples
            for _ in range(refractory_samples):
                try:
                    t = next(it)
                except: #to prevent code from stopping when t reaches end of iterator
                    break
        elif data_sample[t] < baseline - V_th:
            spikes[t] = -1
            DOWN_spikes.append(t)
            baseline -= V_th
            for _ in range(refractory_samples):
                try:
                    t = next(it)
                except: #to prevent code from stopping when t reaches end of iterator
                    break
                
    #converting from sample indeces to spike times
    UP_spike_times = np.array(UP_spikes)/sampling_rate
    DOWN_spike_times = np.array(DOWN_spikes)/sampling_rate
    
    return UP_spike_times, DOWN_spike_times, UP_spikes, DOWN_spikes, spikes, V_th
'''

这个函数 SF_encod 实现了一种称为“步进编码”（Step-Forward Encoding, SF）的编码方法，用于将肌电图（EMG）数据转换成尖峰时间序列。这种编码方法通过逐步调整基线值来确定尖峰，并使用一定的阈值和不应期来生成上尖峰（UP spike）和下尖峰（DOWN spike）。

主要功能
计算不应期样本数：将不应期（以毫秒为单位）转换为对应的采样数。
计算差值：计算相邻样本之间的差值。
确定阈值：基于差值的均值和标准差确定阈值（V_th）。
提取尖峰：
设置初始基线值为数据的第一个样本。
遍历数据样本，如果当前样本值超过基线值加上阈值，则记录为上尖峰，并调整基线值，同时跳过不应期内的样本。
如果当前样本值低于基线值减去阈值，则记录为下尖峰，并调整基线值，同时跳过不应期内的样本。
转换尖峰时间：将尖峰的样本索引转换为时间。
函数参数和返回值说明
输入参数：

data_sample: EMG数据样本。
sampling_rate: EMG电极的采样率。
f: 调整阈值的因子。
refractory: 不应期（毫秒）。
输出：

UP_spike_times: 上尖峰的时间数组。
DOWN_spike_times: 下尖峰的时间数组。
UP_spikes: 上尖峰的样本索引数组。
DOWN_spikes: 下尖峰的样本索引数组。
spikes: 尖峰数据数组（大小与数据样本相同，值为{-1, 0, 1}）。
V_th: 用于信号重建的阈值。

'''

关于不应期样本数的解释：
不应期样本数（refractory samples）是指在神经信号处理中，用于防止连续尖峰（spikes）过于密集的一段时间。具体来说，在一个尖峰检测到后，会有一段时间内不会再检测到新的尖峰，这段时间称为不应期（refractory period）。不应期样本数则是不应期对应的采样数。
不应期在神经信号处理和神经科学中非常重要，因为它模拟了生物神经元在产生一个动作电位（尖峰）后需要一段时间恢复，这段时间内神经元不会再产生新的动作电位。
在编码过程中，不应期样本数的计算方法如下：
不应期（毫秒）：这个值通常是由用户指定的，以控制两个尖峰之间的最小时间间隔。
采样率（Hz）：信号的采样率，表示每秒钟采集的样本数。
不应期样本数：将不应期从毫秒转换为对应的采样数。
计算公式：
不应期样本数 = （不应期 / 1000) * 采样率

refractory_samples = round((refractory / 1000) * sampling_rate)
'''
在代码中，不应期样本数的计算如上：
假设不应期为5毫秒，采样率为1000 Hz（每秒1000个样本），则不应期样本数计算如下：
不应期样本数 = （5/1000）* 1000 = 5
这意味着在检测到一个尖峰后，会跳过接下来的5个样本，在这5个样本内不会再检测到新的尖峰。
不应期样本数在尖峰检测过程中用于跳过不应期内的样本，确保在这段时间内不会再记录新的尖峰，从而避免连续检测到过多尖峰，模拟生物神经元的行为。
'''