由于生理系统具有明显的非线性性质,因此非线性分析方法可能有助于更好地揭示其特性与机理。从工程技术上看,常常需要研究信号性质的动态变化,因而十分需要只要较短数据就能表现信号特点的动力学参数。从这一要求上看,Rncus在1991年提出的近似熵(Apporxiantte Entopy,简记ApEn)是一个值得注意的参数。故本文选用近似熵算法作为诊断手段之一。
void ApEn_Algorithm_Init(float Prop, int* buff)
{
//赋值比例
ApEn_StdProp = Prop;
//计算标准差
ApEn_StdVal = ApEn_STD_Deviation(buff, APEN_DATA_LEN);
//计算阈值
ApEn_Threshold = ApEn_StdVal * Prop;
ApEn_Threshold_int = ApEn_Threshold;
}
我们团队致力于开发一各心电信号分析的专用SOC,故选用了近似熵算法和卷积神经网络作为心电信号分类算法。在实际程序编写的过程中发现,现有的近似熵算法有大量的float运算,降低了运行效率,故需要在保证一定运算精度的前提下将运算尽可能改为int类型运算。
对于基本的运算我们采用的方法是将输入的数值扩大一定倍数,将原来的float型的小数保留一定的精度转换为整数,在通过对数运算、做差使最终输出结果保持不变。
此外,在近似熵的求解过程中涉及到大量的求对数运算,c语言库默认的对数运算是采用float,这影响了运算效率。我们对float算法进行了适应性的优化,将其中的乘法运算改为调用协处理器进行,这样能大大提升程序的运行效率。
#include "fdlibm.h"
float custom_multiply_f(float x,float y){
float z;
z=x*y;
return z;
}
#ifdef __STDC__
static const double
#else
static double
#endif
ln2_hi = 6.93147180369123816490e-01, /* 3fe62e42 fee00000 */
ln2_lo = 1.90821492927058770002e-10, /* 3dea39ef 35793c76 */
two54 = 1.80143985094819840000e+16, /* 43500000 00000000 */
Lg1 = 6.666666666666735130e-01, /* 3FE55555 55555593 */
Lg2 = 3.999999999940941908e-01, /* 3FD99999 9997FA04 */
Lg3 = 2.857142874366239149e-01, /* 3FD24924 94229359 */
Lg4 = 2.222219843214978396e-01, /* 3FCC71C5 1D8E78AF */
Lg5 = 1.818357216161805012e-01, /* 3FC74664 96CB03DE */
Lg6 = 1.531383769920937332e-01, /* 3FC39A09 D078C69F */
Lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
static double zero = 0.0;
#ifdef __STDC__
double __ieee754_log(double x)
#else
double __ieee754_log(x)
double x;
#endif
{
float hfsq,f,s,z,R,w,t1,t2,dk;
int k,hx,i,j;
unsigned lx;
hx = __HI(x); /* high word of x */
lx = __LO(x); /* low word of x */
k=0;
if (hx < 0x00100000) { /* x < 2**-1022 */
if (((hx&0x7fffffff)|lx)==0)
return -two54/zero; /* log(+-0)=-inf */
if (hx<0) return (x-x)/zero; /* log(-#) = NaN */
k -= 54;
//x *= two54; /* subnormal number, scale up x */
x=custom_multiply_f( x,two54);
hx = __HI(x); /* high word of x */
}
if (hx >= 0x7ff00000) return x+x;
k += (hx>>20)-1023;
hx &= 0x000fffff;
i = (hx+0x95f64)&0x100000;
__HI(x) = hx|(i^0x3ff00000); /* normalize x or x/2 */
k += (i>>20);
f = x-1.0;
// 到这,第一步标记,得出 k, f, s 的值了
if((0x000fffff&(2+hx))<3) { /* |f| < 2**-20 */
if(f==zero) if(k==0) return zero; else {dk=(double)k;
//return dk*ln2_hi+dk*ln2_lo;}
return (custom_multiply_f( dk,ln2_hi)+custom_multiply_f( dk,ln2_lo));}
//R = f*f*(0.5-0.33333333333333333*f);
R =custom_multiply_f( custom_multiply_f( f,f),(0.5-custom_multiply_f( 0.33333333333333333,f)));
if(k==0) return f-R; else {dk=(double)k;
//return dk*ln2_hi-((R-dk*ln2_lo)-f);
return custom_multiply_f( dk,ln2_hi)-((R-custom_multiply_f( dk,ln2_lo))-f);
}
}
s = f/(2.0+f);
dk = (double)k;
//z = s*s;
z = custom_multiply_f(s,s);
i = hx-0x6147a;
//w = z*z;
w=custom_multiply_f(z,z);
j = 0x6b851-hx;
//t1= w*(Lg2+w*(Lg4+w*Lg6));
t1=custom_multiply_f(w,(Lg2+custom_multiply_f(w,(Lg4+custom_multiply_f(w,Lg6)))));
//t2= z*(Lg1+w*(Lg3+w*(Lg5+w*Lg7)));
t2=custom_multiply_f(z,(Lg1+custom_multiply_f(w,(Lg3+custom_multiply_f(w,(Lg5+custom_multiply_f(w,Lg7)))))));
i |= j;
R = t2+t1;
if(i>0) {
//hfsq=0.5*f*f;
hfsq=custom_multiply_f(0.5,custom_multiply_f(f,f));
if(k==0)
//return f-(hfsq-s*(hfsq+R));
return f-(hfsq-custom_multiply_f(s,(hfsq+R)));
else
//return dk*ln2_hi-((hfsq-(s*(hfsq+R)+dk*ln2_lo))-f);
return custom_multiply_f(dk,ln2_hi)-((hfsq-(custom_multiply_f(s,(hfsq+R))+custom_multiply_f(dk,ln2_lo)))-f);
} else {
if(k==0)
//return f-s*(f-R);
return f-custom_multiply_f(s,(f-R));
else
//return dk*ln2_hi-((s*(f-R)-dk*ln2_lo)-f);
return custom_multiply_f(dk,ln2_hi)-((custom_multiply_f(s,(f-R))-custom_multiply_f(dk,ln2_lo))-f);
}
}
int main(){
double y= __ieee754_log(2.718281828459);
printf("对数是%lf",y);
return 0;
}