单片机中常用的一些滤波算法

0 引言

在单片机的数据采集系统中，测量通道串入随机干扰是难免的，从而使A/D转换送入单片机的数据存在误差，这种因随机干扰产生的误差称为随机误差。随机误差虽然无法预测，但多次测量结果是符合统计规律的。为克服随机干扰产生的误差，可以采用硬件抗干扰的方法，也可以采用按统计规律用软件方法实现，即采用数字滤波的方法来抑制有效信号中的干扰成分，从而消除随机误差。

1 常用的数字滤波算法及实现

在单片机的数据采集系统中，常用的数字滤波算法有限幅滤波、中位值滤波、算术平均滤波、移动平均滤波、加权平均滤波、低通滤波、中位值平均滤波等。以下就分别介绍其原理及滤波程序。

1.1 限幅滤波

限幅滤波的基本原理是把两次相邻时刻（n和n-1）的采样值Yn和Yn-1相减，求出其差值，以绝对值表示，然后将这个差值与两次采样允许的最大偏差值ΔY比较，如果两次采样值的差值超过了允许的最大偏差值ΔY，则认为发生了随机干扰，并认为最后一次采样值Yn非法，应予剔除。剔除Yn后，可用Yn-1代替Yn；若未超过允许的最大偏差值范围，则认为本次采样值有效。可用如下公式表示：

|Yn-Yn-1|≤ΔY；则Yn有效

|Yn-Yn-1|＞ΔY；则Yn-1有效

此算法的样例子程序如下：

#define A 10   //A值可根据实际情况调整

char data;    //上一次的数据

char filter_1()

{

char datanew;  //新数据变量

datanew=get_data();   //获得新数据

//滤波算法

if ((datanew-data>A)||(data-datanew>A) 

return data;

return datanew;

}

    该算法主要用于处理变化比较缓慢的数据，如温度、物体的位置等。使用时关键在于最大偏差值的Δy的选择，通常可根据经验获得，也可按照输出参数可能的最大变化速度Vmax及采样周期T来决定ΔY的值[1]，即ΔY=VmaxT。

1.2 中位值滤波

中位值滤波是先对某一参数连续采样N次（一般N取奇数），然后把N次采样值按从小到大排列，取中间值为本次采样值。

该滤波方法实际上是一种排序方法，文献[2，7，9]采用的是冒泡法排序。由于在冒泡法排序中，每出现一次前者数据大于后者数据，就要进行二者数据的交换，效率很低，故在验证时改用选择排序法。

该算法的样例子程序如下：

#define N 11  //N值可根据实际情况调整

char filter_2()

{

char value_buf[];

char count,i,j,k,temp;

for(count=0;count<N;count++) //获取数据

{

value_buf[count]=get_data();

delay();

}

for(i=0;i<N-1;i++) //选择排序

 {k=i; 

for(j=i+1;j<N;j++)

      if(value_buf[j]<value_buf[k]) k=j;

      temp=value_buf[k];

      value_buf[k]=value_buf[i];

      value_buf[i]=temp;

      }

return value_buf[(N-1)/2];

}

中位值滤波能有效地克服偶然因素引起的波动或采样器不稳定引起的误码等脉冲干扰。对温度、液位等缓慢变化的被测参数采用此算法能收到良好的滤波效果，但对于流量、压力等快速变化的数据，不宜采用中位值滤波。

1.3 算术平均滤波

算术平均滤波法是要按输入的N个采样数据Xi（i从1到N），寻找这样一个Y，使Y与各采样值之间的偏差的平方和最小，即

                     （1）

由一元函数求极值的原理，可得

                    （2）

此即为算术平均滤波的基本算式。具体实现此算法的子程序如下：

#define N 12

char filter_3()

{

int count;

int sum=0;

for(count=0;count<N;count++)

 {

  sum+=get_ad();

  delay();

  }

 return (char)(sum/N);

}

算术平均滤波适用于对一般具有随机干扰的信号进行滤波。这种信号的特点是有一个平均值，信号在某一数值范围附近做上下波动，在这种情况下仅取一个采样值做依据显然是不准确的。算术平均滤波对信号的平滑程序完全取决于N，当N较大时，平滑度高，但灵敏度低；当N较小时，平滑度低，但灵敏度高，应视具体情况选取N，以便既少占用计算时间，又达到最好的效果。

1.4 移动平均滤波

算术平均滤波需要连续采样若干次后，才能进行运算而获得一个有效的数据，因而速度较慢。为了克服这一缺点，可采用移动平均滤波。即先在RAM中建立一个数据缓冲区，按顺序存放N次采样数据，然后每采进一个新数据，就将最早采集的数据去掉，最后再求出当前RAM缓冲区中的N个数据的算术平均值。这样，每进行一次采样，就可计算出一个新的平均值，即测量数据取一丢一[3]，测量一次便计算一次平均值，大大加快了数据处理的能力。具体的滤波程序如下：

#define N 12

char value_buf[N];

char i=0;

char filter_4()

{

char  count;

int sum=0;

value_buf[i++]=get_ad();

if (i==N) i=0;

for (count=0;count<N;count++)

sum=value_buf[count];

return(char)(sum/N);

}

移动平均滤波对周期性干扰有良好的抑制作用，平滑度高，灵敏度低，但对偶然出现的脉冲性干扰抑制作用差，因此它不适用于脉冲干扰比较严重的场合，而适用于高频振荡的系统。通过观察不同N值下移动平均的输出响应来选取N值，以便既少占用单片机的时间，又能达到最好的滤波效果。

1.5 加权平均滤波

在算术平均滤波和移动平均滤波中，N次采样值在输出结果中的权重是均等的，取1/N。用这样的滤波算法，对于时变信号会引入滞后，N值越大，滞后越严重[1,4]。为了增加新采样数据在移动平均中的权重，以提高系统对当前采样值中所受干扰的灵敏度，可采用加权平均滤波，它是移动平均滤波算法的改进。

加权平均滤波是对连续N次采样值分别乘上不同的加权系统之后再求累加和，加权系统一般先小后大，以突出后面若干采样的效果，加强系统对参数变化趋势的辨识。各个加权系统均为小于1的小数，且满足总和等于1的约束条件。这样，加权运算之后的累加和即为有效采样值。

为方便计算，可取各加权系数均为整数，且总和为256，加权运算后的累加和除以256，即舍去低字节后便是有效采样值。具体的样例子程序如下：

//code数组为加权系统表，存在ROM区。

#define N 12

char code jq[N]={1,2,3,4,5,6,7,8,9,10,11

,12}; 

char code sum_jp=1+2+3+4+5+6+7+8+9+10+11

+12;

char filter_5()

{

char count;

char value_buf[N];

int sum=0;

for(count=0;count<N;count++)

{

value_buf[count]=get_data();

delay();

}

for(count=0;count<N;count++)

sum+=value_buf[count]*jq[count];

return (char)(sum/sum_jq);

}

1.6 中位值平均滤波

该算法的基本原理是，对连续的n个数据进行排序，去掉其中最大和最小的2个数据，将剩余数据求平均值。

为了加快数据处理速度，n可以取值6，而对于具有较快速度的处理器，n值可以适当取大一些。但最好是n=2k+2[5]，k为整数,因为这样可以在求平均值时用移位的方法，加快处理速度。

上述算法存在的一个不足之处就是每采集一个数据就要进行一次排序，这样会大量占用系统宝贵的时间，文献[8]的算法程序即是如此。

可以通过存储当前数据中的最大值和最小值来改进，具体做法是：

系统中用两个变量来存储当前n个数据中最大值和最小值在数组中的偏移量。这样只有在当前输入的数据将要覆盖的数据正好是当前的最大值或最小值时，才在下个数组中查找最大值或最小值，而其他情况下，则只要将输入的数据与最大值和最小值比较就可以修改最大值和最小值了，而且不用进行数据排序。

改进后的算法程序如下：

#define dtype unsigned int 

#define uint8 char

#define LEN  6 

#define SHIFT 2

uint8 pdata;

uint8 pmax,pmin; 

dtype datas[LEN];

dtype szlb(dtype _data)

{ 

uint8 i;

dtype average=0;

pdata=(pdata+1)%LEN; 

datas[pdata]=_data; 

for(i=0;i<LEN;i++)

average+=datas[i]; 

if(_data>datas[pmax])

pmax=pdata;

else if(_data<datas[pmin])

pmin=pdata; 

if(pdata==pmax) 

 { for(i=0;i<LEN;i++)

   if(datas[i]>datas[pmax])

   pmax=i;

 }

 else if(pdata==pmin) 

{

 for(i=0;i<LEN;i++)

 if(datas[i]<datas[pmin])

 pmin=i;

 }

average=average-datas[pmax]-datas[pmin]; 

 return (average>>SHIFT);  

}

这种滤波方法兼容了移动平均滤波算法和中位值滤波算法的优点，所以无论对缓慢变化的信号，还是对快速变化的信号，都能取得较好的滤波效果。

1.7低通滤波

将描述模拟RC低通滤波特性的微分方程用差分方程来表示，便可以用软件算法来实现模拟滤波的功能[4]。

最简单的一阶RC低通滤波器描述其输入x(t)与输出y(t)的微分方程为

                （3）

以采样周期T对x(t)、y(t)进行采样得

yn=y(Nt)

xn=x(Nt)

如果T<<RC，则由微分方程可得差分方程

               （4）

令                    （5）

可得yn=axn+(1-a)yn-1,其中：

xn：当前的数据

yn-1：上次的滤波输出值

a：滤波系统

yn：本次的滤波输出值

由上式可以看出，本次滤波的输出值主要取决于上次滤波的输出值，这和加权平均滤波是有本质区别的，本次采样值对滤波输出的贡献是比较小的，但多少有些修正作用，这种算法便模拟了具有较大惯性[6]的低能滤波功能。当目标参数为变化很慢的物理量时，这是很有效的。

低通滤算法程序和加权平均滤波相似，但加权系统只有两个：a和1-a，因为只有两项，a和1-a均以立即数的形式编入程序中，不另设表格。虽然采样值为单元字节，为保证运算精度，滤波输出值用双字节表示，其中一字节整数，一字节小数，否则有可能因为每次舍去尾数而使输出不会变化。

低通滤波的程序如下：

//为加快程序处理速度，假定基数为100，//a=0～100

#define a 50

char value;

char filter_6()

{

char new_value;

new_value=get_data();

return(100-a)*value+a*new_value;

}

一阶低通滤波算法对周期性干扰具有良好的抑制作用，适用于波动频繁的参数滤波，其不足之处是带来了相位滞后，灵敏度低。滞后的程度取决于a值的大小。同时，它不能滤除频率高于采样频率1/2（奈奎斯特频率）的干扰信号[7]，对于高于奈奎斯特频率的干扰信号，应该采用模拟滤波器。