肌电信号时域特征和频域特征提取，特征带有标签_时域特征ssc公式

import numpy as np
from feature_utils import *
import math
from keras.utils import to_categorical
 
arr = normalized_filtered_emg
new_shape = (-1, 12)
emg = np.reshape(arr, new_shape)
featureData = []
featureLabel = []
classes = 49
timeWindow = 200
strideWindow = 100
 
for i in range(classes):
    index = []
    for j in range(label.shape[0]):
        if np.all(label[j,:]==i):
            index.append(j)
    iemg = emg[index,:]
    length = math.floor((iemg.shape[0] - timeWindow) / strideWindow)
    print("class ", i, ", number of sample: ", iemg.shape[0], length)
 
    for j in range(length):
        start_index = strideWindow * j
        end_index = start_index + timeWindow
        if end_index <= iemg.shape[0]:
            rms = featureRMS(iemg[start_index:end_index,:])
            mav = featureMAV(iemg[start_index:end_index,:])
            wl = featureWL(iemg[start_index:end_index,:])
            zc = featureZC(iemg[start_index:end_index,:])
            ssc = featureSSC(iemg[start_index:end_index,:])
            featureStack = np.hstack((rms, mav, wl, zc, ssc))
 
            featureData.append(featureStack)
            featureLabel.append(i)
 
featureData = np.array(featureData)
featureLabel = np.array(featureLabel)
 
print(featureData.shape)
print(featureLabel.shape)
 

对特征的定义可以在GITHUB中找到

def spectrum(signal, fs):
    m = len(signal)
    n = next_power_of_2(m)
    y = np.fft.fft(signal, n)
    yh = y[0:int(n / 2 )]
    fh = (fs / n) * np.arange(0, n / 2, 1)#
    power = np.real(yh * np.conj(yh) / n)#
    #print(fh.shape)
    #print(power.shape)
    return fh, power
def next_power_of_2(x):
    return 1 if x == 0 else 2 ** (x - 1).bit_length()
 
def frequency_ratio(frequency, power):
    power_low = power[(frequency >= 30) & (frequency <= 250)]
    power_high = power[(frequency > 250) & (frequency <= 500)]
    ULC = np.sum(power_low)
    UHC = np.sum(power_high)
 
    return ULC / UHC
def mean_freq(frequency, power):#频率和功率
    num = 0
    den = 0
    for i in range(int(len(power) / 2)):
        num += frequency[i] * power[i]
        den += power[i]      
        return num / den
 
 
def median_freq(frequency, power):
    power_total = np.sum(power) / 2
    temp = 0
    tol = 0.01
    errel = 1
    i = 0
 
    while abs(errel.all()) > tol:
        temp += power[i]
        errel = (power_total - temp) / power_total
        i += 1
        if errel.all() < 0:
            errel = 0
            i -= 1
 
    return frequency[i]
 
for i in range(classes):
    index = []
    for j in range(label.shape[0]):
        if np.all(label[j, :] == i):
            index.append(j)
    iemg = emg[index, :]
    length = math.floor((iemg.shape[0] - timeWindow) / strideWindow)
    print("class ", i, ", number of sample: ", iemg.shape[0], length)
 
sampling_frequency = 2e3
classes = 49
timeWindow = 200
strideWindow = 200
emg=normalized_filtered_emg
featureData = []
featureLabel = []
 
for j in range(length):
    start_index = strideWindow * j
    end_index = start_index + timeWindow
    if end_index <= iemg.shape[0]:
            # 提取频率特征
        fh, power = spectrum(iemg[start_index:end_index, :], fs)  # 使用您的频率特征提取函数
        fr = (frequency_ratio(fh, power))  # 使用您的频率特征提取函数
        mnp = (np.sum(power) / len(power)) # 使用您的频率特征提取函数
        tot = (np.sum(power))  # 使用您的频率特征提取函数
        mnf = (mean_freq(fh, power))  # 使用您的频率特征提取函数
        pkf = (fh[np.argmax(power)])  # 使用您的频率特征提取函数
 
 
        featureStack = np.hstack((fr, mnp, tot, mnf, pkf))
 
        featureData.append(featureStack)
        featureLabel.append(i)
 
featureData = np.array(featureData)
featureLabel = np.array(featureLabel)
 
print(featureData.shape)
print(featureLabel.shape)

对频域特征提取，采用与时域相同代码可以保证特征提取时，有对应标签。