使用了GPU的基于Python的MSE(多尺度熵)
这边用了cupy库来进行GPU加速,相对于CPU版本有了明显的加速
测试环境是Linux Mint 20.2 Cinnamon,GPU为GT730
在Windows上由于未知原因不能使用cupy
cupy具体内容请百度
代码
MSE部分
import numpy as np
import cupy as cp
from matplotlib import pyplot as plt
import time
from SampEn import sampEn
import math
def MSE(signal , max_scale:int = 20):
result = cp.zeros(max_scale)
length = len(signal)
std = cp.std(signal)
for scale in range(1 , max_scale + 1):
# 确定截取的长度
length = int(len(signal) / scale) - 1
# 分段取平均
scale_i = cp.array(signal[ : len(signal) : scale][:length])
for i in range(1,scale):
scale_i = scale_i + signal[i: len(signal) : scale][:length]
scale_i = scale_i / scale
#计算样本熵
result[scale - 1] = sampEn(scale_i, std)
print("scale:" , scale, 'SampEn' , result[scale - 1] )
return result
white_noise = np.loadtxt("white_noise.csv")
# RR_Data = np.loadtxt("RR_Inrerival_healthy_000.txt")
white_noise = cp.array(white_noise)
begin = time.time()
entropy = MSE(white_noise[1:50000] , 20)
# entropy = MSE(RR_Data)
end = time.time()
样本熵部分(SampEn)
import numpy as np
import cupy as cp
import time
def sampEn(L:cp.array, std : float = 1,m: int= 2, r: float = 0.15):
"""
计算时间序列的样本熵
Input:
L: 时间序列
std: 原始序列的标准差
m: 1或2
r: 阈值
Output:
SampEn
"""
N = len(L)
B = 0.0
A = 0.0
# Split time series and save all templates of length m
xmi = cp.array([L[i:i+m] for i in range(N-m)])
xmj = cp.array([L[i:i+m] for i in range(N-m+1)])
# Save all matches minus the self-match, compute B
B = cp.sum(cp.array([cp.sum(cp.abs(xmii-xmj).max(axis=1) <= r * std)-1 for xmii in xmi]))
# Similar for computing A
m += 1
xm = cp.array([L[i:i+m] for i in range(N-m+1)])
A = cp.sum(cp.array([cp.sum(cp.abs(xmi-xm).max(axis=1) <= r * std)-1 for xmi in xm]))
# Return SampEn
return -cp.log(A/B)
if __name__ == "__main__":
white_noise = np.loadtxt("white_noise.csv")
white_noise = cp.array(white_noise)
# RR_Data = np.loadtxt("RR_Inrerival_healthy_000.txt")
begin = time.time()
entropy = sampEn(white_noise[1:50000])
# entropy = MSE(RR_Data)
end = time.time()
print("用时:" , int((end - begin)/60) ,'min', (end - begin) % 60 , 's')
print(entropy)