📄 ecg_data.cs
字号:
if(Flag[index]==3)
{
// ************************ index值为当前R峰位置 ****************************************
// 检测中使用的极限位置为当前和下一个R峰，Q峰，S峰位置，不使用Q波，S波起点和终点是因为
// 之前的算法中R峰，Q峰，S峰位置识别较为精确，Q波起点和S波终点可能会出现较大偏差
// ****************************************************************************************
int NextRcrest=index+1;
while(Flag[NextRcrest]!=3 && NextRcrest<length-1) NextRcrest++; // 下一个R峰波峰
double max=0;
double maxY=0; //用来记录T波波峰近似位置（最大斜率处）
int Tonset=0;
int Toffset=0;
int Tcrest=index+120/SampleInterval; //从当前R峰后120ms位置开始搜索T波
if(index+120/SampleInterval>length-1)break; //如果R峰后120ms内就到达了数据的结尾.说明最后的R峰后没有T波;跳出最外层的for循环
else max=Y[index+120/SampleInterval];
if(NextRcrest==length-1) //没有下一个R峰了，就从当前R峰后120ms开始寻找最大值，其位置定为T波峰值位置
{
max=0;
for(int search=index+120/SampleInterval;search<=length-1;search++)
{
if(max<=Math.Abs(MeanAtPosition(search)))
{max=Math.Abs(MeanAtPosition(search));Tcrest=search;}
}
Flag[Tcrest]=7;
index=length-1;
break; //结束整个检测函数
}
int Screst=index+1;
while(Flag[Screst]!=4) Screst++; // 当前S波波峰
int NextQcrest=index+1;
while(Flag[NextQcrest]!=2)NextQcrest++; // 下一个Q波波峰
int boundary=450/SampleInterval;
do
{
boundary -=50/SampleInterval;
for(int search=index+120/SampleInterval;search<index+boundary && search< NextQcrest-3;search++) //搜索当前峰值位置的后150ms到400ms之间，限制搜索位置不超过下一个Q波峰值位置，寻找最大斜率点,这是比较接近T波峰顶的位置,
{
if(maxY<Y[search])
{maxY=Y[search];Tcrest=search;}
}
}while(Y[Tcrest]>0.4*Y[NextRcrest-3] && boundary>0);
max=Math.Abs(MeanAtPosition(Tcrest)); //校正T波顶点位置
for(int search=Tcrest-60/SampleInterval;search<=Tcrest+80/SampleInterval && search> Screst+5 && search< NextQcrest-5;search++) //从准T波波峰位置的前60ms后80ms区间内搜索真正峰顶位置
{
if(max<=Math.Abs(MeanAtPosition(search)))
{max=Math.Abs(MeanAtPosition(search));Tcrest=search;}
}
Flag[Tcrest]=7; //标记T波顶点为7
int k=0;
double dt=0;
while(k<10)
{
dt+=0.025;
for( Toffset=Tcrest+2;Toffset<=Tcrest+150/SampleInterval && Toffset<NextQcrest-1;Toffset++) //从T波峰顶位置开始的后150ms区间内搜索T波终点
if(Y[Toffset]<=maxY*dt && Y[Toffset+1]<=maxY*dt)
{
Flag[Toffset]=8; //标记T波终点为8
k=10;
break;
}
k++;
}
k=0;
dt=0.05;
while(k<5)
{
dt+=0.05;
for( Tonset=Tcrest-100/SampleInterval;Tonset<=Tcrest-2 ;Tonset++) //从T波峰顶位置开始的前100ms区间内搜索T波起点
if(Y[Tonset]<=maxY*dt && Y[Tonset-1]<=maxY*dt)
{
Flag[Tonset]=6; //标记T波起点为6
k=10;
break;
}
k++;
}
}//end if(Flag[index]==3)
SingleLeadProcessed=true;
}//end if(QRSProcessed)
}//end PTDetect()
#endregion
#region Marr小波变换
public void WaveletTransform()
{
if(DataInputed)
{
WTcoefficients=new double[6,length]; //不同尺度下小波变换系数数组
SmoothedSignal=new double[6,length]; //不同尺度下平滑近似信号数组
for(int i=0;i<length;i++)
SmoothedSignal[0,i]=SourceData[i]; //尺度为1(2的零次幂)时存贮原始信号
double[] h={0.0032 ,-0.0132, -0.0393, 0.0450 ,0.2864 ,0.4347, 0.2864 ,0.0450 ,-0.0393 ,-0.0132 ,0.0032}; //等价低通滤波器系数
double[] g={0.0039, -0.0062, -0.0226 ,-0.1120 ,-0.2309 ,0.7118 ,-0.2309 ,-0.1120 ,-0.0226 ,-0.0062 ,0.0039}; //等价高通滤波器系数
for(int scale=1;scale<=6;scale++)
{
for(int position=160;position<length-160;position++)
{
for(int k=-5;k<=5;k++)
{
SmoothedSignal[scale,position] += h[k+5]*SmoothedSignal[scale-1,position-(int)Math.Pow(2,scale-1)*k];
WTcoefficients[scale,position] += g[k+5]*SmoothedSignal[scale-1,position-(int)Math.Pow(2,scale-1)*k];
}
}
}
MaxMoudle=new double [6,length];
WTProcessed=true;
}
*/
}
#endregion
#region 数据归一化处理
public void Normalzie ()
{
if(DataInputed)
{
double maxdata=SourceData[0];
for(int index=0;index<length;index++)
if(maxdata<SourceData[index])maxdata=SourceData[index]; //找到最大数据
double mindata=SourceData[0];
for(int index=0;index<length;index++)
if(mindata>SourceData[index])mindata=SourceData[index]; //找到最小数据
for(int index=0;index<length;index++) //规一化
SourceData[index]=(SourceData[index]-mindata)/(maxdata-mindata);
}
}
#endregion
#region 二次样条小波变换
public bool SplineWT()
{
if(DataInputed)
{
WTcoefficients=new double [6,length];
double [,] SmoothedSignal=new double [6,length];
for(int index=0;index<length;index++)
{
SmoothedSignal[0,index]=SourceData[index]; //2的零次幂尺度上的平滑近似信号为原始信号
WTcoefficients[0,index]=SourceData[index];
}
double [] h={0 ,0.125 ,0.375, 0.375 ,0.125 }; //低通滤波器系数
double [] g={0 ,0 ,-2.0 , 2.0 ,0 }; //高通滤波器系数
for(int Scale=1;Scale<5;Scale++) //完成从2的一次幂到2的5次幂尺度上的小波变换
{
for(int n=0;n<length;n++)
{
int offset=0;
for(int k=-2;k<=2;k++)
{
offset=n-(int)Math.Pow(2,Scale-1)*k;
if(offset>=0 && offset<length)
{
SmoothedSignal[Scale,n] += h[k+2]*SmoothedSignal[Scale-1,offset] ;
WTcoefficients[Scale+1,n] += g[k+2]*SmoothedSignal[Scale-1,offset] ;
}
}
}
}
MaxMoudle=new int[5,length];
Flag=new int[length];
WTProcessed=true; //已完成小波变换
return true;
}
else
return false;
}
#endregion
#region 拉德梅克函数
private double rad(int i,double t)
{
return Math.Sign(Math.Sin( Math.Pow(2,i)* Math.PI *t /length));
}
#endregion
#region 沃尔什函数
private double wal(int n,double t)
{
if(n==0)return 1;
if(n==1)return rad(1,t);
if(n==2)return rad(1,t)*rad(2,t);
if(n==3)return rad(2,t);
if(n==4)return rad(2,t) *rad(3,t);
if(n==5)return wal(4,t)*wal(1,t);
if(n==6)return rad(1,t) *rad(3,t);
if(n==7)return rad(3,t);
if(n==8)return rad(3,t)*rad(4,t);
if(n==9)return rad(1,t)*wal(8,t);
if(n==10)return wal(2,t) *wal(8,t);
if(n==11)return wal(3,t) *wal(8,t);
if(n==12)return wal(3,t) *rad(4,t);
if(n==13)return wal(2,t) *rad(4,t);
if(n==14)return rad(1,t) *rad(4,t);
if(n==15)return rad(4,t);
if(n==16)return rad(4,t)*rad(5,t);
if(n==17)return rad(1,t)*rad(4,t)*rad(5,t);
if(n==18)return wal(17,t)*rad(2,t);
if(n==19)return rad(2,t)*rad(4,t)*rad(5,t);
if(n==20)return wal(19,t)*rad(3,t);
if(n==21)return rad(1,t) *rad(2,t) *rad(3,t) *rad(4,t) *rad(5,t);
if(n==22)return rad(1,t) *rad(3,t) *rad(4,t) *rad(5,t);
if(n==23)return rad(3,t) *rad(4,t) *rad(5,t);
if(n==24)return rad(3,t) *rad(5,t);
if(n==25)return rad(3,t) *rad(1,t) *rad(5,t);
if(n==26)return wal(25,t) *rad(2,t);
if(n==27)return rad(3,t) *rad(2,t) *rad(5,t);
if(n==28)return rad(2,t) *rad(5,t);
if(n==29)return rad(1,t) *rad(2,t) *rad(5,t);
if(n==30)return rad(1,t) *rad(5,t);
if(n==31)return rad(5,t);
if(n==32)return rad(5,t) *rad(6,t);
if(n==33)return rad(1,t)*rad(5,t) *rad(6,t);
else return 0;
}
#endregion
#region 离散沃尔什变换
public bool WalshTransform()
{
if(DataInputed)
{
double []Fn=new double [34]; //生成34个列率的离散沃尔什变换列率谱系数
for(int n=0;n<34;n++)
Fn[n]=F(n);
for(int index=0;index<length;index++)
{
WTcoefficients[3,index]=0; //把沃尔什变换的系数存在了WTcoefficients[3, ]对应位置上。这样做的理由是
//最终图形输出的时候输出的是WTcoefficients数组不同尺度上数据的图形
//所以用WTcoefficients[0,]存贮原始信号；WTcoefficients[1,]WTcoefficients[2,]
//存贮1，2尺度上小波变换的结果；WTcoefficients[3,]存储沃尔什变换系数
for(int n=0;n<=33;n++) //截至列率到34，运算量很大，截至列率越高，变换系数的步长越短。
WTcoefficients[3,index]+=Fn[n]*wal(n,index);
}
WalshTtansformProcessed=true;
return true;
}
else
return false;
}
#endregion
#region 离散沃尔什变换列率谱系数
private double F(int n)
{
double Fn=0;
for(int i=0;i<length;i++)
Fn += SourceData[i]*wal(n,i);
Fn/=length;
return Fn;
}
#endregion
#region 由沃尔什变换系数判定层届点
public void Delaminate ()
{
if(WalshTtansformProcessed)
{
for(int index=0;index<length-1;index++)
{
int onset=0,offset=0; //含水层起点，终点
double OnsetThresh=0.08; //起点阈值
while( index<=length-2 && (WTcoefficients[3,index+1]-WTcoefficients[3,index])<OnsetThresh )index++; //寻找突然上升的点（上升高度要足够高）
if(index<=length-1)onset=index; //标记为含水层起点（要求不超过数组的长度）
while( index<length-2 && (WTcoefficients[3,index+1]-WTcoefficients[3,index])>=0 )index++; //寻找到上升的尽头（开始下降的位置）
bool stop=false;
offset=index;
for(;!stop;index++)
{
while( index<length-2 && (WTcoefficients[3,index+1]-WTcoefficients[3,index])<0 )index++;//追踪到下降的尽头（遇到平台了）记下平台起点（可能为含水层终点）
offset=index;
while( index<length-2 && (WTcoefficients[3,index+1]-WTcoefficients[3,index])==0 )index++;//考察平台的另一头是上升还是下降
if(index<length-2 && (WTcoefficients[3,index+1]-WTcoefficients[3,index])>0 && WTcoefficients[3,index]<0.47) {stop =true;index-=2;} //上升的话。平台的起点就是含水层终点
else
if(index>=length-2)stop=true;
}//end for(;!stop;index++)
if(onset<=length-1 && offset<=length-1 && offset-onset <30 &&offset -onset>=2)
{
double max=WTcoefficients[0,onset];
for(int i=onset;i<=offset;i++)
{
if(max<WTcoefficients[0,i])max=WTcoefficients[0,i];
}
if(max>Threshold)
{Flag[onset]=1;Flag[offset]=2;}
}//end if(onset<=length-1 &&
}//end for(int index=0;index<length-1;index++
}//end if(WalshTtansformProcessed)
}
#endregion
}//end class
}//end namespace
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -