心电滤波(STM32)小结

心电滤波(STM32)小结

 心电信号经过心电电路的放大滤波后，仍存在着50HZ工频干扰，基线漂移，高斯噪声等干扰信号。还可以在软件层面通过一系列算法滤除，本文总结了深圳暑期培训期间，通过陷波器滤除50Hz工频干扰，高通滤波去除基线漂移，以及FIR滤波，光滑滤波去除高斯噪声，这四种滤波算法

目录

心电滤波(STM32)小结

心电信号

陷波器

相关参数

c代码

高通滤波

相关参数

c代码：

FIR滤波

相关参数

c代码：

光滑滤波

c代码：

心率测量

---

---

原始信号（500uf工频干扰，1000uf基线漂移）

算法滤波后的信号

心电信号

心电信号由心电模拟仪或人体采集，通过心电滤波放大电路后，由STM32的ADC进行采样，并进行串口打印输出，采样的频率和周期如下：

采样频率Fs=500

周期T=1/Fs=0.002

陷波器

陷波器是一种特定类型的滤波器，其主要功能是从信号中去除或极大衰减特定频率的成分，也就是消除特定频率的“陷波”。陷波器在通信、信号处理和控制工程等领域具有广泛应用。

在本实验中，使用陷波器滤除信号中由电源引起的50Hz工频干扰。

相关参数

陷波频率Fc=50Hz

差分方程：

y[n] = a0 * x[n] + a1 * x[n-1] + a2 * x[n-2] - b1 * y[n-1] - b2 * y[n-2]

相关系数：

A[3]    = {1 , --1.618 ,    1};

B[3]    = {1 , --1.5533,    0.9216};

c代码

动态滤波

uint16_t Notch(uint16_t x1)

{

    uint16_t y = 0;

    static uint16_t y_last = 0;

    static uint16_t y_last_last = 0;

    static uint16_t x_last = 0;

    static uint16_t x_last_last = 0;

    y=x1 +x_last_last - 1.618*x_last + 1.5533*y_last - 0.9216*y_last_last;

    x_last_last=x_last;

    x_last=x1;

    y_last_last = y_last;

    y_last = y;

    return y;

}

静态滤波（L是数组的长度）

void Notch(uint16_t *Array,uint16_t* ecg1)

{

    uint16_t y = 0;

    uint16_t y_last = 0;

    uint16_t y_last_last = 0;

    uint16_t x_last = 0;

    uint16_t x_last_last = 0;

    for(m=0;m<L;m++)    

    {

        y=Array[m] +x_last_last - 1.618*x_last + 1.5533*y_last - 0.9216*y_last_last;

        y_last_last = y_last;

        y_last = y;

        x_last_last=x_last;

        x_last=Array[m];

        ecg1[m]=y;

    }

}

高通滤波

相关参数

截至频率Fc=0.5Hz

差分方程

y[n]=(b0*x[n]+b1*x[n-1]+b2*x[n-2])*Gain-a1*y[n-1]-a2*y[n-2]

相关系数

A[3]    = {1 , -1.98430 ,    0.98438156};

B[3] = {1 , -1.999993 , 1};

Gain = 0.99217 ;

c代码：

动态滤波

double HighFilter(double indata)
{
    double outdata = 0;
    

    double B[3] = {1 , -1.999993004365251403342540470475796610117 , 1};
    double A[3]    = {1 , -1.984302608109661525404021631402429193258 ,    0.984381559862382404801905977365095168352};
    double Gain = 0.992172777212600220941851603129180148244 ;

    static    double w_x[3] ={0 , 0 , 0};
    static    double w_y[3] ={0 , 0 , 0};
        
    w_x[0] = indata;
    w_y[0] = (B[0] * w_x[0] + B[1] * w_x[1] + B[2] * w_x[2]) * Gain - w_y[1] * A[1] - w_y[2] * A[2];
    
    outdata = w_y[0]/A[0];
    
    w_x[2]=w_x[1];w_x[1]=w_x[0];
    w_y[2]=w_y[1];w_y[1]=w_y[0];
 
    return outdata;
}

FIR滤波

相关参数

截至频率Fc=62.5Hz

滤波器阶数degree=4；

差分方程

y[n]=b0*x[n]+b1*x[n-1]+b2*x[n-2]+b3*x[n-3]+b4*x[n-4]

相关系数

b[5]={0.0246 ,  0.2344 , 0.4821 , 0.2344 ,  0.0246 }

c代码：

动态滤波

uint16_t FIR_ECG(uint16_t output)

{

    static uint16_t x[5];

    uint16_t newdata;

    x[0]=output;    

    for(m=4;m>0;m--)

    {

        x[m]=x[m-1];

    }

    newdata = b0*x[0]+b1*x[1]+b2*x[2]+b3*x[3]+b4*x[4];

    return newdata;

}

静态滤波（L是数组的长度）

void FIR_ECG(uint16_t *input,uint16_t* output)

{

    uint16_t x[5];

    uint16_t j ,m;

    for(m=0;m<L;m++)

    {

        x[0]=input[m];

        output[m]=b0*x[0]+b1*x[1]+b2*x[2]+b3*x[3]+b4*x[4];

        for(j=4;j>0;j--)

        {

            x[j]=x[j-1];

        }

    }

}

光滑滤波

光滑滤波通过移动窗口，取前10个数的均值来达到滤波效果。其主要作用是减小或去除信号中的高频噪声或不规则波动，从而使信号变得更加平稳和稳定。

c代码：

动态滤波

static u16 sliding_average_filter(u16 value)
{
  static int data_num = 0;
  static u16 output = 0;
    int i;
  float sum = 0;

  if (data_num < 10) //不满窗口，先填充
  {
    buffer[data_num++] = value;
    output = value; //返回相同的值
  }
  else
  {
    memcpy(&buffer[0], &buffer[1], (int)(10 - 1) * 4); //将1之后的数据移到0之后，即移除头部
    buffer[10 - 1] = value;                       //即添加尾部

    for (i = 0; i < 10; i++) //每一次都计算，可以避免累计浮点计算误差
    {    
      sum += buffer[i];
    }
    output = sum / 10;
  }

  return output;
}

心率测量

        作者通过静态心电信号测量来测量心率，即先在数组里面保存好心电信号，传递数组来寻找峰值，计算峰值之间的间距来测量心率。因为通过陷波器滤波之后的心电信号前小段有一些类似正弦波的无效值，所以一定要把这段去除后，再进行峰值判断才准确

//寻找峰值
uint16_t maxECG(uint16_t *Ecg)
{
	uint16_t max;
	max = Ecg[0];
	for(m=0;m<1000;m++)
	{
		if(max<Ecg[m])
		{
			max=Ecg[m];
		}
	}	
	return max;
}

//测量心率
uint8_t H_Rate(uint16_t *Ecg)
{
	uint16_t r_time;
	uint8_t rate;
	uint16_t point;
	uint8_t flagg = 0;
	uint16_t point_x1 = 0;
	uint16_t point_x2 = 0;
	point=maxECG(Ecg)-100;
	for(n=0;n<L;n++) 
	{
		if(Ecg[n]>point &&flagg==0 )
		{
			flagg=1;
			point_x2 = point_x1;
			point_x1=n;			
			if(point_x2 !=0 )
			{		
				r_time=(point_x1 - point_x2);
				rate = 30000/r_time;
				return rate;
			}
		}
		if(flagg==1  && Ecg[n]<point)
		{
			flagg=0;
		}
	}
	return 255;
}