常用的滤波算法

大傻子的文渊阁 2019-07-17

展开全文

滤波是传感器处理中的重要算法，经常接触底层常常用到，以下总结了一些滤波算法，供以后参考调用。

一、低通滤波

1.1RC滤波的数字低通滤波

　　指在截止频率fc的时候，增益为-3db(Aup=0.707)的滤波器，也是模电书中出现的第一种硬件滤波器，以下是对应的软件形式的1阶RC滤波器的数字形式(本断程序节选自匿名4轴)

　　一阶形式：Y(n)=(1-a)*Y(n-1)+a*X(n)

　　下式中 oldData表示上一次的输出Y(n-1) newData表示新的输入X(n)

 1 　　float LopPassFilter_RC_1st(float oldData, float newData, float a)
 2 　　{
 3 　　　　return oldData * (1 - a) + newData * a;
 4 　　}
 5 
 6 　　计算比例系数a:
 7 
 8 　　float LopPassFilter_RC_1st_Factor_Cal(float deltaT, float Fcut)
 9 　　{
10 　　　　return deltaT / (deltaT + 1 / (2 * M_PI * Fcut));
11 　　}

1.2均值滤波：

　　把一段时间内的数据累加后求平均值，达到平滑的作用，适用性广泛，元素越多滤波效果越好时延越高。

 1 　　uint16_t LowPassFilter_Average(uint16_t data[],uint16_t length)
 2 
 3 　　{
 4 
 5 　　　　uint32_t add=0;
 6 　　　　uint16_t result;
 7 　　　　int i;
 8 
 9 　　　　for(i=0;i<length;i++)
10 　　　　{
11 　　　　　　add += data[i];
12 　　　　}
13 　　　　result=add/length;
14 　　　　return result;
15 　　}
16 
17 　　//data[]放入一段时间里的数值，length：data数组的长度

1.3滑动滤波

　　在均值滤波的基础上，加上比例系数，最新的数据具有更大的比例，增加时效性。

　　以下列举了一个三个元素的滑动滤波

 1 　　double LowPassFilter_Silding(double dataNewest,double dataMiddle,double dataLast)
 2 　　{
 3 
 4 　　　　#define PROPORTIONNEW 0.55
 5 
 6 　　　　#define PROPORTIONMID 0.35
 7 
 8 　　　　#define PROPORTIONLAST (1- PROPORTIONNEW -PROPORTIONMID )
 9 
10 　　　　double result;
11 　　　　result = PROPORTIONNEW *dataNewest+ PROPORTIONMID *dataMiddle+PROPORTIONLAST *dataLast;
12 　　　　return result;
13 　　}

　　三个输入参数为连续3个读出的数据，其中 dataNewest 为最新数据，对应比例系数PROPORTIONNEW ，另外两个同理。比例系数根据需要改变

1.4中值滤波

　　当传感器采集的数据存在毛刺时，为了提取其中有效的数据，采用中值滤波的算法，只保存数据大小在中间的数值。

 1 　　uint16_t LowPassFilter_GetMiddleValue( uint16_t *pusBuffer, uint8_t ucLength )
 2 　　{
 3 　　　　uint8_t i;
 4 　　　　uint16_t usMiddle,usMax,usCount = 0,usPoint;
 5 
 6 　　　　while( usCount <= (ucLength/2) )
 7 　　　　{
 8 　　　　　　usPoint = 0;
 9 　　　　　　usMax = pusBuffer[0];
10 　　　　　　for( i=1;i<ucLength;i++ ) 
11 　　　　　　{ 
12 　　　　　　　　if( pusBuffer[i] > usMax ) 
13 　　　　　　　　{
14 　　　　　　　　　　usMax = pusBuffer[i];
15 　　　　　　　　　　usPoint = i;
16 　　　　　　　　} 
17 　　　　　　}
18 　　　　　　usMiddle = usMax;
19 　　　　　　pusBuffer[usPoint] = 0; 
20 　　　　　　usCount++; 
21 　　　　}
22 　　　　return usMiddle;
23 　　}

　　*pusBuffer:待改变的起始地址 ucLength:长度

　　特别注意：为了运行效率，这个方法中没有对数组进行备份，会直接操作到源 *data并改变元素的值。

　　这个快速中值算法原理是只排序最大或者最小值的一半，一半中的最后一个即为中间值

1.5去除极值的XX滤波

　　把中值滤波的优点结合给其他滤波，使其他滤波器也具有够剔除少量毛刺的能力。

　　以下例子结合的是均值滤波，即去除极值的均值滤波

 1 　　int16_t LowPassFilter_RemoveExtremumAverage(int16_t data[],int16_t length)
 2 　　{
 3 　　　　#define EXTERMUM_NUMBER 1 //number of extermum
 4 　　　　int32_t sum=0;
 5 　　　　int16_t i;
 6 
 7 　　　　if(length<=2*EXTERMUM_NUMBER)
 8 　　　　{
 9 　　　　　　return 0;//length shorter than 2*EXTERMUM_NUMBER, return err
10 　　　　}
11 
12 　　　　//sort
13 　　　　BubbleSort_int16(data,length);
14 　　　　
15 
16 　　　　//average filter
17 　　　　for(i=EXTERMUM_NUMBER;i<length-EXTERMUM_NUMBER;i++)
18 　　　　{
19 　　　　　　sum+=data[i];
20 　　　　}
21 
22 　　　　return sum/(length-2*EXTERMUM_NUMBER);
23 　　}

　　输入长度为length的数组data[]。EXTERMUM_NUMBER定义去除极值的个数（1表示去掉1个最大值和1个最小值）。先用排序法（本例使用下面给出的冒泡排序）对data数组进行排序，去除1个最大值和1个最小值后剩余求平均值。

　　这段代码可能出现两个错误

　　1、length比2倍的EXTERMUM_NUMBER短，本例用return 0一笔带过，实际需要在程序return 0处作处理

　　2、如果输入length过长会使32位的sum在累加过程溢出，如果length很长建议改变均值滤波的顺序在累加之前线进行除法

　　下面给出子程序冒泡排序

 1 　　void BubbleSort_int16(int16_t *a,int len)
 2 　　{
 3 　　　　int i;
 4 　　　　int j;
 5 　　　　int mid;
 6 　　　　for(i=0;i<len;i++)
 7 　　　　{
 8 　　　　　　for(j=0;j<len-i-1;j++)
 9 　　　　　　{
10 　　　　　　　　if(*(a+j)>*(a+j+1))
11 　　　　　　　　{
12 　　　　　　　　　　mid=*(a+j);
13 　　　　　　　　　　*(a+j)=*(a+j+1);
14 　　　　　　　　　　*(a+j+1)=mid;
15 　　　　　　　　}
16 　　　　　　}
17 　　　　}
18 　　}

1.6巴特沃斯低通滤波

　　借助MATLAB使用数字滤波的方法来模拟巴特沃斯滤波器：

　　IIR滤波的基本结构 y(n)+a1*y(n-1)+a2*y(n-2)....ai*y(n-i)=b0*x(n)+b1*x(n-1)+b2*x(n-2)....+bi*x(n-i);

　　这里简单把IIR滤波理解为滑动滤波的进阶形式，其中i为阶数，x(n)为本次输入，x(n-1)为上一次输入，如此类推x(n-i)为i次前的输入。y(n)为输出，y(n-1)为上一次的输出。y(n-i)为i次前的输出。 a和b分别为比例。

　　系数a和b改变达到模拟不同采样频率、不同截止频率的滤波效果。

　　以下使用matlab计算一个1000HZ采样频率，100HZ截止频率，2阶的巴特沃斯滤波器的IIR系数：

　　>>fdatool

1.7切比雪夫低通滤波

1.8卡尔曼滤波

卡尔曼滤波是一种预测型滤波，这里把卡尔曼滤波归纳到低通滤波中，是因为在适用过程中卡尔曼增益Kg会慢慢收敛，变成跟低通滤波一样的效果

 1 float LowPassFilter_kalman(float data)
 2 { 
 3 　　static float Xlast=0.01;//初值 Xlast&Xpre
 4 　　static float P=0.1;//Plast&Pnow
 5 　　const float Q=0.44;//自己感觉的方差
 6 　　const float R=0.54;//测量器件的方差 
 7 　　float Kg,PP;//
 8 　　float Xnow;//经过卡尔曼滤波的值 Xnow
 9 
10 　　//1式A=1 无输入
11 　　PP=P+Q; //Ppre=Plast+Q （2）
12 　　Kg=PP/(PP+R); //更新Kg Ppre/Ppre+R（4）
13 　　Xnow=Xlast+Kg*(data-Xlast); // Xnow = Xpre+Kg*(Z(k)-H*Xpre)(3)
14 　　P=(1-Kg)*PP; //Pnow=(I-Kg)*Ppre(5) 由于Pnow不会再使用，所以更新Plast=Pnow
15 　　Xlast=Xnow; 
16 
17 　　return Xnow;
18 }

二、高通滤波

三、融合滤波

3.1互补滤波

当有两种同一类型参数的时候，虽然不同它们的特性又刚好互有长短，为了各取所长，互补其短，所以叫互补滤波。以下为这种情况的举例：

在姿态角测量中加速度计（偏移小、精度低）和陀螺仪（偏移大、发散、精度高）

在高度测量中超声波传感器（偏移小、采样率低）和气压计（偏移大、采样率高）

在位置测量中 GPS（偏移小、采样率低）和惯性导航系统（偏移大、采样率高）

以下这段程序以姿态角计算为例说明

提取Pitch轴说明：陀螺仪读出角速度积分A，利用重力加速度读出加速度计可以计算出夹角B

制作工艺使得陀螺仪积分角度A随时间的推移发生零点偏移（角度产生严重的偏差），积分后产生越来越严重的累计误差。

另一方面，加速度计算夹角B没有积分过程不会产生如上问题，当物体运动时，物体本身的加速度没有被消除导致加速度计的值在物体加减速运动瞬间偏差。

故使用互补滤波让陀螺仪A为精确测量部件占一个大的比例，加速度B占一个小的比例作为修正

 1 float FusionFilter_Complementary_1p(float AccAngle,float GyrRate,double PRYAngleLast,float dtC)
 2 {
 3 　　float tau=3;
 4 　　float porprtion;
 5 　　float RPYAngleNew;
 6 
 7 　　porprtion=tau/(tau+dtC); 
 8 　　RPYAngleNew=porprtion*(PRYAngleLast+GyrRate*dtC) + (1-porprtion)*(AccAngle);
 9 　　return RPYAngleNew;
10 }

AccAngle:输入加速度计算夹角后的值，GyrRate：角速度，PRYAngleLast：上一次的姿态角（Pitch或Roll适用），dtC：两次调用的时间间隔用于积分

3.2卡尔曼滤波

与上面低通滤波中的卡尔曼滤波不同，这次的输入参数包括两种，以一种对另一种进行预测，以下是一段修改自互联网的程序，使用加速度计和陀螺仪来计算姿态角

 1 void FusionFilter__Kalman_Filter_IMU(float Accel,float GyroRate,double *Angle ,double dt) 
 2 {
 3 
 4 static const float Q_angle=0.001; 
 5 static const float Q_gyro=0.003;
 6 static const float R_angle=0.5;
 7 //    static const float dt=0.01;    
 8 static const char C_0 = 1;
 9 static float Q_bias, Angle_err;
10 static float PCt_0, PCt_1, E;
11 static float K_0, K_1, t_0, t_1;
12 static float Pdot[4] ={0,0,0,0};
13 static float PP[2][2] = { { 1, 0 },{ 0, 1 } };
14 
15 *Angle +=(GyroRate - Q_bias) * dt; //先验估计
16 
17 Pdot[0]=Q_angle - PP[0][1] - PP[1][0]; // pk-先验估计误差协方差的微分
18 
19 Pdot[1]= -PP[1][1];
20 Pdot[2]= -PP[1][1];
21 Pdot[3]=Q_gyro;
22 
23 PP[0][0] += Pdot[0] * dt; //pk-先验估计误差协方差微分的积分
24 PP[0][1] += Pdot[1] * dt; //=先验估计误差协方差
25 PP[1][0] += Pdot[2] * dt;
26 PP[1][1] += Pdot[3] * dt;
27 
28 Angle_err = Accel - *Angle;    //zk-先验估计
29 
30 PCt_0 = C_0 * PP[0][0];
31 PCt_1 = C_0 * PP[1][0];
32 
33 E = R_angle + C_0 * PCt_0;
34 
35 K_0 = PCt_0 / E;
36 K_1 = PCt_1 / E;
37 
38 t_0 = PCt_0;
39 t_1 = C_0 * PP[0][1];
40 
41 PP[0][0] -= K_0 * t_0;    //后验估计误差协方差
42 PP[0][1] -= K_0 * t_1;
43 PP[1][0] -= K_1 * t_0;
44 PP[1][1] -= K_1 * t_1;
45 
46 *Angle    += K_0 * Angle_err;    //后验算估计
47 Q_bias    += K_1 * Angle_err;    //后验估计
48 //    Gyro_y = GyroRate - Q_bias;    //输出值（后验估计）的微分 = 实际角速度
49 }