ecg_data.cs

这个是用沃尔什变换对测井曲线的处理程序

using System;
using System.IO;
using System.Drawing;
using System.Collections;
using System.ComponentModel;
using System.Windows.Forms;
using System.Data;


namespace ECG_test
{
	/// <summary>
	/// ECG_DATA 的摘要说明。
	/// </summary>
	public class ECG_DATA
	{
		public bool DataInputed;         //防止用户在数据没有输入之前处理数据
		public bool WTProcessed;         //小波变换是否处理，测试用变量
		public bool WalshTtansformProcessed;  //沃尔什变换是否处理，测试用变量
		public bool QRSProcessed;        //防止用户在数据没有处理之前调用心率等参数
		public bool PTProcessed;         //P波T波是否处理，测试用变量

		public double[] SourceData;       //存储原始信号数据
		public double[,] WTcoefficients;  //不同尺度下小波变换系数数组，数组0尺度存贮原始信号
		//public double[,] SmoothedSignal;  //不同尺度下平滑近似信号数组
		public int[,] MaxMoudle;          //模极大值序列,当对应的值为1时为极大值。其他的点为零
		public int TotalHeartBeatNumber;  //总心跳数即总QRS波群数
		public int SampleRate;            //采样率
		private int SampleInterval;       //采样的时间间隔（单位为毫秒ms）
		public int[] Flag;                //对于Flag标记数组。Q波起点设为1,Q波顶点设为2,R峰顶点设为3,S波顶点设为4,S波终点设为5,T波起点设为6,T波顶点设为7,T波终点设为8
		private double []Slope;           //存储数据斜率
		private int length;               //数据长度
		public double Threshold;          //阈值（针对测井数据筛选含水层时使用，只有视电阻率高于该值被确定为含水层，需要归一化）

		#region   数据长度只读访问器
		public int Length                 //定义只读访问器，防止用户修改数据长度
		{
			get 
			{
				if(DataInputed) return length;
				else            
					return 0;
			}
		}
		#endregion

		#region 心率只读访问器
		public double HeartRate           //心率,只读访问器
		{
			get 
			{
				if(QRSProcessed)return (1000.0/MeanRRInterval)*60;
				else             return  0;
			}
		}
		#endregion

		#region 平均RR间期只读访问器
		public double MeanRRInterval      //平均RR间期,单位为毫秒，只读访问器
		{
			get
			{
				if(QRSProcessed)
				{
					int index1=0,index2=0;
					while(Flag[index1]!=3)index1++;
					while(Flag[length-index2-1]!=3)index2++;
					return (length-1-index1-index2)/TotalHeartBeatNumber*1000/SampleRate;
				}
				else 
					return 0;
			}
		}
		#endregion

		#region 类构造函数
		public ECG_DATA()                 //构造函数
		{
			DataInputed =false;
			WTProcessed=false;
			QRSProcessed=false;
			PTProcessed=false;
			Threshold=0.3;
			SampleRate=0;
			length=0;
			TotalHeartBeatNumber=0;
  
		}
		#endregion

		#region 打开文件函数
		public string OpenFile()            //读文件，输入12导联心电数据
		{
			OpenFileDialog file=new OpenFileDialog();
			file.Title="选取ECG测试文件";
			file.InitialDirectory = "C:\\Documents and Settings\\Mg\\桌面\\测井\\数据" ;
			file.Filter = "ecg files (*.txt)|*.txt|All files (*.*)|*.*" ;
			file.FilterIndex = 1 ;
			file.RestoreDirectory = true ;
            
			if(file.ShowDialog() == DialogResult.OK)
			{
				DataInputed=true;
				WTProcessed=false;
				QRSProcessed=false;
				PTProcessed=false;
				WalshTtansformProcessed=false;
				TotalHeartBeatNumber=0;
				StreamReader sr = new StreamReader(file.FileName);
				ArrayList values = new ArrayList();
				string line = sr.ReadLine();
				while (line != null)
				{
					values.Add(Convert.ToDouble(line));
					line = sr.ReadLine();
				}
				length=values.Count;

				SampleRate=360;         //MIT/BIH数据库心电信号统一采样率
				SampleInterval=1000/360;
				SourceData=new double[length];
				for(int index=0;index<length;index++)
					SourceData[index]=(double)values[index];
				return file.FileName;
			}
			else 
			{
				DataInputed=false;
				return null;
			}//end else

		}
		#endregion

		#region Marr小波变换
		public bool MarrWT()
		{
			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.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++)                      //完成从2的一次幂到2的5次幂尺度上的小波变换
				{
					for(int n=0;n<length;n++)
					{
						int offset=0;
						for(int k=-5;k<=5;k++)
						{
							offset=n-(int)Math.Pow(2,Scale-1)*k;
							if(offset>=0 && offset<length)
							{
								SmoothedSignal[Scale,n] += h[k+5]*SmoothedSignal[Scale-1,offset] ;
								WTcoefficients[Scale,n] += g[k+5]*SmoothedSignal[Scale-1,offset] ; 
							}
						}
					}
				}
                MaxMoudle=new int[5,length];
				Flag=new int[length];
				WTProcessed=true;     //已完成小波变换
				return true;
			}
			else
				return false;
		}
		#endregion

		#region 生成指定尺度上的模极大值序列
		private bool MaxModuleAtScale(int Scale)
		{
			if(Scale<1 ||Scale>5) return false;
			
			Slope=new double [length]; 
			Slope[1]=Slope[0]=Slope[length-1]=Slope[length-2]=0;//WTcoefficients[Scale,index]-WTcoefficients[Scale,index-1];//
			for(int index=2;index<length-3;index++)	
				Slope[index] =-2*WTcoefficients[Scale,index-2]-WTcoefficients[Scale,index-1]+
					WTcoefficients[Scale,index+1]+2*WTcoefficients[Scale,index+2];
			double max=0;      //以下的模块初始化max的值。
		
		  {
			 double [] maxk=new double[20];
			 for(int k=0;k<20;k++)
				for(int index=k*180+1;index<=180*(k+1);index++)
					if(maxk[k]<Math.Abs(Slope[index]))
						maxk[k]=Math.Abs(Slope[index]);
			 double temp=0;
			 for(int i=0;i<19;i++)           //20个最大值排序
			 {
				temp=0;
				for(int j=0;j<=i;j++)
					if(maxk[j]>maxk[j+1])
					{
						temp=maxk[j];
						maxk[j]=maxk[j+1];
						maxk[j+1]=temp;
					}
			 }
			 for(int i=19;i>=0;i--)
			 {
				double mean=0;
				for(int j=0;j<=i;j++)
					mean +=maxk[j];
				mean /=i;
				//max=mean;
			   if(maxk[i]<=2.0*mean)
					max=maxk[i];
				if(maxk[i]<=2.0*mean)
					break;
			 }
		   }
           
			
			int pre =0;
			for(int index=0;index<length-1;index++)
			{
			 double thresh=0.5*max;
			  if(Math.Abs(Slope[index])>=thresh && Math.Abs(Slope[index+1])>=thresh)
			  {
				 double temp=WTcoefficients[Scale,index];
				 int start=index-5;
				 int end =index+15;
				 int position=index;
				 if(start<0)start=0;
				 if(end>length-1)end=length-1;
				 double NewMax=Math.Abs(Slope[index]);
				 for(int k=start;k<=end;k++)   //在斜率最大值附近寻找系数极大值和斜率最大值
				 {
					 if(temp<WTcoefficients[Scale,k])
					 {temp=WTcoefficients[Scale,k];position =k;}
					 if(NewMax<Math.Abs(Slope[k]))
						 NewMax=Math.Abs(Slope[k]);
					
				 }
				 
				  pre=position;
				 if(Math.Abs(WTcoefficients[0,position])>=0.5*Math.Abs(WTcoefficients[0,pre]) && Math.Abs(WTcoefficients[0,position])<=2.0*Math.Abs(WTcoefficients[0,pre]))
				  { 
					 MaxMoudle[Scale-1,position]=1;          //MaxMoudle 仅为5行
					  Flag[position]=3;           //标记
					  TotalHeartBeatNumber++;
					  if(NewMax<1.3*max && NewMax>0.5*max) //index>20*180 && 
						  max=max*0.9+0.1*NewMax;
					  thresh=max*0.5;
					  index=position+75;
				  }

			 }//end if(Math.Abs(Slope[index])>=thresh && Math.Abs(Slope[index+1])>=thresh)
				if(index-pre>=2000/SampleInterval)
				{index=pre+75;max=0.5*max;}

			}//end for(int index=0;index<length-1;index++)
			return true;
		}
		#endregion

		#region QRS波群特征点检测
		public int  QRSDetect()         //返回值为QRS波群的个数  
		{
			if(DataInputed && WTProcessed)
			{
				//Flag=new int[length];//对于Flag数组。Q波起点设为1,Q波顶点设为2,R峰顶点设为3,S波顶点设为4,S波终点设为5
                MaxModuleAtScale(4);
				/*MaxModuleAtScale(2);
				MaxModuleAtScale(4);
				for(int index=0;index<length;index++)
				{
					if(MaxMoudle[3,index]==1)
					{
						int start=index-5;
						int end =index+5;
						int position=index;
						if(start<0)start=0;
						if(end>length-1)end=length-1;
						bool condition1=false;
						bool condition2=false;
						for(int i=start;i<end;i++)
						{
							if(MaxMoudle[2,i]==1)condition1=true;
                            if(MaxMoudle[4,i]==1)condition2=true;
						}
						if(condition1 && condition2)
						{
							Flag[index]=3;TotalHeartBeatNumber++;
						}
						else
							MaxMoudle[3,index]=0;
					}
				}*/
				QRSProcessed=true;
				return TotalHeartBeatNumber;
				
			}//end if(DataInputed && WTProcessed)
			else        //数据没有输入就不错处理
			{
				QRSProcessed=false;
				return 0;
			}

		}// end QRSDetect()
		#endregion

		#region T波特征点检测
		public void PTDetect()            
		{
			/*if(QRSProcessed)
			{
				PTProcessed=true;
				for(int index=0;index<length-1;index++)