`
几种常用的滤波算法
在嵌入式开发中经常会用到一些滤波算法，我整理了一些资料把这些算法封装成可以直接调用的函数，方便以后的开发，包括限幅滤波、中位值滤波法、 算术平均滤波法、一阶滞后滤波法、加权递推平均滤波法、消抖滤波法 ，其他的一些包括卡尔曼滤波什么的之后有机会再整理
#include "filter.h"



/*
限幅滤波

　　A方法： 根据经验判断，确定两次采样允许的最大偏差值（设为A），每次检测到新值时判断： 如果本次值与上次值之差<=A，则本次值有效，如果本次值与上次值之差>A，则本次值无效，放弃本次值，用上次值代替本次值。

　　B优点： 能有效克服因偶然因素引起的脉冲干扰。

　　C缺点： 无法抑制那种周期性的干扰，平滑度差。

		A值可根据实际情况调整

		value为有效值，new_value为当前采样值

		滤波程序返回有效的实际值
*/

#define DEVIATION  10

float limit_filter(float new_value)
{
    static int num = 0;
    static float value = 0;  //需要赋一个初值
    num ++;
    if(num == 1)
    	value = new_value;
    else
    {
	    if ( ( new_value - value > DEVIATION ) || ( value - new_value > DEVIATION ))
	        return value;
	}
    return new_value;
}


/*
中位值滤波法

　　A方法： 取之前采样的N次（N取奇数），把N次采样值按大小排列，取中间值为本次有效值。

　　B优点： 能有效克服因偶然因素引起的波动干扰，对温度、液位的变化缓慢的被测参数有良好的滤波效果。

　　C缺点： 对流量、速度等快速变化的参数不宜。

    排序采用冒泡法 只需要移动最后一个元素即可

*/
#define MIDDLE_FILTER_N  11

float middle_filter( float new_value)
{
    static float value_buf[MIDDLE_FILTER_N];
    float temp ;
    uint8_t count, i;
    for ( count = 0; count < MIDDLE_FILTER_N - 1; count++)
    {
        value_buf[count] = value_buf[count + 1] ;
    }
    value_buf[MIDDLE_FILTER_N - 1] = new_value;

    for (i = MIDDLE_FILTER_N; i > 0 ; i --)
    {
        if ( value_buf[i] < value_buf[i - 1] )
        {
            temp = value_buf[i];
            value_buf[i] = value_buf[i - 1];
            value_buf[i - 1] = temp;
        }
    }
    return value_buf[(MIDDLE_FILTER_N - 1) / 2];
}


/*
   算术平均滤波法

		A方法： 连续取N个采样值进行算术平均运算，N值较大时：信号平滑度较高，但灵敏度较低；N值较小时：信号平滑度较低，但灵敏度较高。N值的选取：一般流量，N=12；压力：N=4。

　　B优点： 适用于对一般具有随机干扰的信号进行滤波，这样信号的特点是有一个平均值，信号在某一数值范围附近上下波动。

　　C缺点： 对于测量速度较慢或要求数据计算速度较快的实时控制不适用，比较浪费RAM 。

   无需每次求一编所有的和，减去第一个数据加上新数据

*/

#define AVERAGE_N 12

float average_filter(float new_value)
{
    static float average_value_buf[AVERAGE_N];
    static  float average_sum = 0;
    uint8_t count;
    average_sum -= average_value_buf[0];
    for ( count = 0; count < AVERAGE_N - 1; count++)
    {
        average_value_buf[count] = average_value_buf[count + 1] ;
    }
    average_value_buf[AVERAGE_N - 1] = new_value;
    average_sum += average_value_buf[11];

    return (average_sum /(AVERAGE_N * 1.0) );
}



/*
一阶滞后滤波法
		A方法： 取a=0~1，本次滤波结果=（1-a）*本次采样值+a*上次滤波结果。

　　B优点： 对周期性干扰具有良好的抑制作用，适用于波动频率较高的场合。

　　C缺点：相位滞后，灵敏度低，滞后程度取决于a值大小，不能消除滤波频率高于采样频率的1/2的干扰信号。

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

#define FIRST_LAG_PROPORTION 0.4

float first_order_lag_filter( float new_value)
{
    static float first_order_value , first_order_last_value;
    first_order_value = first_order_last_value;
    first_order_last_value = new_value;
    return (1 - FIRST_LAG_PROPORTION) * first_order_value + FIRST_LAG_PROPORTION * new_value;
}


/*
加权递推平均滤波法

　　A方法： 是对递推平均滤波法的改进，即不同时刻的数据加以不同的权，通常是，越接近现时刻的资料，权取得越大，给予新采样值的权系数越大，则灵敏度越高，但信号平滑度越低。

　　B优点： 适用于有较大纯滞后时间常数的对象和采样周期较短的系统。

　　C缺点： 对于纯滞后时间常数较小，采样周期较长，变化缓慢的信号，不能迅速反应系统当前所受干扰的严重程度，滤波效果差。

		coe数组为加权系数表，存在程序存储区。
*/

#define WEIGHT_AVERAGE_N 12



uint8_t  coe[WEIGHT_AVERAGE_N] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
uint8_t  sum_coe = 1 + 2 + 3 + 4 + 5 + 6 + 7 + 8 + 9 + 10 + 11 + 12;

float weighted_filter(float new_value)
{
    static float weight_average_buf[WEIGHT_AVERAGE_N];
    uint8_t count;
    float  sum = 0;
    for ( count = 0; count < AVERAGE_N - 1  ; count++)
    {
        weight_average_buf[count] = weight_average_buf[count+ 1] ;
    }
    weight_average_buf[AVERAGE_N - 1] = new_value;

    for (count = 0 ; count < WEIGHT_AVERAGE_N; count++)
        sum += weight_average_buf[count] * coe[count];
    return (sum / (sum_coe * 1.0));
}

/*
消抖滤波法
　　A方法： 设置一个滤波计数器，将每次采样值与当前有效值比较： 如果采样值＝当前有效值，则计数器清零。如果采样值<>当前有效值，则计数器+1，并判断计数器是否>=上限N(溢出)，如果计数器溢出，则将本次值替换当前有效值，并清计数器。

　　B优点： 对于变化缓慢的被测参数有较好的滤波效果，可避免在临界值附近控制器的反复开/关跳动或显示器上数值抖动。

　　C缺点： 对于快速变化的参数不宜，如果在计数器溢出的那一次采样到的值恰好是干扰值，则会将干扰值当作有效值导入系统。
*/

#define SHAKE_N 12

float shake_filter( float new_value , float now_value)
{
    static uint8_t count = 0;
    if(now_value != new_value)
    {
        count++;
        if (count >= SHAKE_N)
        {
            count = 0;
            return new_value;
        }
    }
    return now_value;
}

常用空间滤波处理分为两大类：线性滤波和非线性滤波。OpenCV里有这些滤波的函数，使用起来非常方便，现在简单介绍其使用方法。
线性滤波:
1.方框滤波：模糊图像
2.均值滤波：模糊图像
3.高斯滤波：信号的平滑处理，去除符合正太分布的噪声
非线性滤波
1.中值滤波：去除椒盐噪声
2.双边滤波：保边去噪
OpenCV中的函数：

1、方框滤波——boxFilter()函数

2、均值滤波——blur()函数

3、高斯滤波——GaussianBlur()函数

4、中值滤波——medianBlur()函数

5、双边滤波——bilateralFilter()函数
参考博客

几种视差图后处理方法，包括滤波去噪（中值滤波或双边滤波）、连通域检测和左右一致性检测。

1）滤波去噪

滤波去躁主要用于去除视差图中由于误匹配造成的孤立噪点，视差图后处理中常用的两种滤波方法有中值滤波和双边滤波。

中值滤波是一种典型的非线性滤波技术，其基本思想是用像素点邻域灰度值的中值来代替该像素点的灰度值，该方法在去除脉冲噪声、椒盐噪声的同时又能保留图像的边缘细节。如下图所示：



 

双边滤波是结合图像的空间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，达到保边去噪的目的。双边滤波器顾名思义比高斯滤波多了一个高斯方差sigma－d，它是基于空间分布的高斯滤波函数，所以在边缘附近，离的较远的像素不会太多影响到边缘上的像素值，这样就保证了边缘附近像素值的保存。如下图所示：



2）连通域检测

连通域检测主要用于去除视差图中由于误匹配造成的小团块。其算法大概流程如下：

0.首先设定speckleRange和speckleWindowSize，speckleWindowSize是指设置检测出的连通域中像素点个数，也就是连通域的大小。speckleRange是指设置判断两个点是否属于同一个连通域的阈值条件。

1.判断当前像素点四邻域的邻域点与当前像素点的差值diff，如果diff<speckleRange，则表示该邻域点与当前像素点是一个连通域，设置一个标记。然后再以该邻域点为中心判断其四邻域点，步骤同上。直至某一像素点四邻域的点均不满足条件，则停止。

2.步骤1完成后，判断被标记的像素点个数count，如果像素点个数count<=speckleWindowSize，则说明该连通域是一个小团块(blob)，则将当前像素点值设置为newValue(表示错误的视差值，newValue一般设置为负数或者0值)。否则，表示该连通域是个大团块，不做处理。同时建立标记值与是否为小团块的关系表rtype[label]，rtype[label]为0，表示label值对应的像素点属于小团块，为1则不属于小团块。

3.处理下一个像素点时，先判断其是否已经被标记：如果已经被标记，则根据关系表rtype[label]判断是否为小团块（blob），如果是，则直接将该像素值设置为newValue；如果不是，则不做处理。继续处理下一个像素。如果没有被标记，则按照步骤1处理。

4.所有像素点处理后，满足条件的区域会被设置为newValue值。

代码：

typedef cv::Point_<short> Point2s;
void myFilterSpeckles(cv::Mat &img, int newVal, int maxSpeckleSize, int maxDiff)
{
	int width = img.cols;
	int height = img.rows;
	int imgSize = width * height;
	int *pLabelBuf = (int*)malloc(sizeof(int)*imgSize);//标记值buffer
	Point2s *pPointBuf = (Point2s*)malloc(sizeof(short)*imgSize);//点坐标buffer
	uchar *pTypeBuf = (uchar*)malloc(sizeof(uchar)*imgSize);//blob判断标记buffer
	//初始化Labelbuffer
	int currentLabel = 0;
	memset(pLabelBuf, 0, sizeof(int)*imgSize);

	for (int i = 0; i < height; i++)
	{
		float *pData = img.ptr<float>(i);
		int *pLabel = pLabelBuf + width * i;
		for (int j = 0; j < width; j++)
		{
			if ((pData[j] - newVal) > 1e-1)
			{
				if (pLabel[j])
				{
					if (pTypeBuf[pLabel[j]])
					{
						pData[j] = newVal;
					}
				}
				else
				{
					Point2s *pWave = pPointBuf;
					Point2s curPoint(j, i);
					currentLabel++;
					int count = 0;
					pLabel[j] = currentLabel;
					while (pWave >= pPointBuf)
					{
						count++;
						float *pCurPos = &img.at<float>(curPoint.y, curPoint.x);
						float curValue = *pCurPos;
						int *pCurLabel = pLabelBuf + width * curPoint.y + curPoint.x;
						//bot
						if (curPoint.y < height - 1 && !pCurLabel[+width] && pCurPos[+width] != newVal && abs(curValue - pCurPos[+width]) <= maxDiff)
						{
							pCurLabel[+width] = currentLabel;
							*pWave++ = Point2s(curPoint.x, curPoint.y + 1);
						}
						//top
						if (curPoint.y > 0 && !pCurLabel[-width] && pCurPos[-width] != newVal && abs(curValue - pCurPos[-width]) <= maxDiff)
						{
							pCurLabel[-width] = currentLabel;
							*pWave++ = Point2s(curPoint.x, curPoint.y - 1);
						}
						//right
						if (curPoint.x < width - 1 && !pCurLabel[+1] && pCurPos[+1] != newVal && abs(curValue - pCurPos[+1]) <= maxDiff)
						{
							pCurLabel[+1] = currentLabel;
							*pWave++ = Point2s(curPoint.x + 1, curPoint.y);
						}
						//left
						if (curPoint.x > 0 && !pCurLabel[-1] && pCurPos[-1] != newVal && abs(curValue - pCurPos[-1]) <= maxDiff)
						{
							pCurLabel[-1] = currentLabel;
							*pWave++ = Point2s(curPoint.x - 1, curPoint.y);
						}

						--pWave;
						curPoint = *pWave;
					}

					if (count <= maxSpeckleSize)
					{
						pTypeBuf[pLabel[j]] = 1;
						pData[j] = (float)newVal;
					}
					else
					{
						pTypeBuf[pLabel[j]] = 0;
					}
				}
			}
		}
	}

	free(pLabelBuf);
	free(pPointBuf);
	free(pTypeBuf);
}


3)左右一致性检测

左右一致性检测是用于遮挡检测的，具体做法：

根据左右两幅输入图像，分别得到左右两幅视差图。对于左图中的一个点p，求得的视差值是d1，那么p在右图里的对应点应该是(p-d1)，(p-d1)的视差值记作d2。若|d1-d2|>threshold，p标记为遮挡点(occluded point)。

在工业过程控制系统中,由于被控对象的环境比较恶劣,干扰源比较多,仪器、仪表采集的信息常会受到干扰,所以在模拟系统中,为了消除干扰,常采用RC滤波电路,而在由工业控制计算机组成的自动检测系统中,为了提高采样的可靠性,减少虚假信息的影响,常常采用数字滤波的方法。

数字滤波的方法有很多种,可以根据不同的测量参数进行选择。下面给出几种常用的数字滤波方法的C语言函数,这些函数有一定的通用性,适用于PC机及其兼容机。 

1.程序判数滤波 

采样的信号,如因常受到随机干扰传感器不稳定而引起严重失真时,可以采用此方法。方法是:根据生产经验确定两交采样允许的最大偏差△×,若先后两次采样的信号相减数值大于△×,表明输入的是干扰信号,应该去掉;用上次采样值作为本次采样值,若小于、等于△×表明没有受到干扰,本次采样值效。该方法适用于慢变化的物理参数的采样,如温度、物理位置等测量系统。

程序判断滤波的C程序函数如下:

float program_detect_filter(float old_new_value[], float X)
{
   float sample_value;
   if (fabs(old_new_value[1] - old_new_value[0])>X)
      sample_value=old_new_value[0];
   else
      sample_value=old_new_value[1];
   retrun(sample_value);
}

 函数调用需一个一维的两个元素的数组(old_new_value[2],用于存放上次采样值(old_new_value[0],)和本次采样值(old_new_value[1],),函数中sample_value表示有效采样值,X表示根据根据经验确定的两次采样允许的最大偏差△x。

2、中值滤波

中值滤汉是对某一参数连续输入N次(一般N取奇数),从中选择一个中间值作为本次采样值,若变量变化比较缓慢,采用此方法效果比较好,但对快速变化过程的参数,如流量、自然伽玛等,则不宜采用。

中值滤波的C程序函数如下:

unsigned char GetMedianNum(int * bArray, int iFilterLen)  
{  
   int i,j;// 循环变量  
   unsigned char bTemp;        
    // 用冒泡法对数组进行排序  
   for (j = 0; j < iFilterLen - 1; j ++)  
   {  
      for (i = 0; i < iFilterLen - j - 1; i ++)  
      {  
          if (bArray[i] > bArray[i + 1])  
         {  
             // 互换  
            bTemp = bArray[i];  
            bArray[i] = bArray[i + 1];  
            bArray[i + 1] = bTemp;  
         }  
      }  
  }  
  // 计算中值  
  if ((iFilterLen & 1) > 0)  
  {  
     // 数组有奇数个元素，返回中间一个元素  
     bTemp = bArray[(iFilterLen + 1) / 2];  
  }  
 else  
  {  
     // 数组有偶数个元素，返回中间两个元素平均值  
     bTemp = (bArray[iFilterLen / 2] + bArray[iFilterLen / 2 + 1]) / 2;  
  }  
  return bTemp;  
}  

函数假设对某一参数连续采样3次,若多次采样,可 对该函数稍作修改即可。3次采样值存储在数组mddle_value[3],其中Sample-value表示有效采样值,count表示连续采样次数。

3、滑动算术平均值滤波

滑动算术平均值滤波是设一循环队列,依顺序存放N次采样数据,每次数据采集时,先将放在队列中第一个最早采集的数据丢掉,再把新数据放入队尾,然后求包括新数据在内的N个数据的算术平均值,便得到该次采样的有效数据。该方法主要用于对压力、流量等周期脉动的采样值进行平滑加工处理。

滑动算术平均值滤波C程序函数如下:

float move_average_filtaer(float data_buf[], int count)
{
   float sample_vaue,data=0;
   int i;
   for (i=0;i<cout;i++)
      data+=data_buf[i];
   sample_value=data/count;
   return(sample_value);
}

函数假设顺序存放5次采样数据的数据缓冲区data_buf[5],对于多于5次的滑动算术平均滤波,只需对该函数稍作修改即可,其中sample_value表示本次采样的有效数据,count表示数据有样次数。

4.滑动加权平均值滤波

滑动加权平均滤波是设一个数据缓冲区依顺序存放N次采样数据,每采进一个新数据,就将最先采集的数据丢掉,而后求包括新数据在内的N个数据的加权平均值,便得到该次采样的有效数据。该方法对脉冲性干扰的平滑作用尚不理想,不适用于脉冲性干扰比较严重的场合。

滑动加权平均值滤波的C程序函数如下:

floa move_times_filter(float data _buf [])
{
   float sample_value;
   float filter_k[3]={0.3 , 0.2 , 0.1 , 0.3};
   sample_value=filter_k[0]*data_buf[1]+filter_k[1]*(data_buf[2]
       +filter_k[2]*data_buf[3]+filter_k[3]*data_buf[4];
   sample_value = sample_value/4;
   return(sample_value);
}

    函数假设依次存放5次采样数据的数据缓冲区data_buf[5],对于多于5次的滑动加权平均滤波,只需对该函数稍作修改即可,其中数据组filter_k[3]表示加权系数,这三个系数的关系为filter_k[0]+2*filter_k[1]+2*filter_k[2]=1,本次采样的有效数据用sample_value表示。

5、防脉冲干扰平均值滤波

防脉冲干扰平均值滤波是连续进行N次采样,去掉其中最大值和最小值,然后求剩下的N-2个数据的平均值,作为本次采样的有效值。该方法适用于变量跳变比较严重的场合。这种滤波也应用边采样边计算的方法。

float max_min_chioce(float x_buffer[],int number)
{
   int max_value, min_value;
   float sample_value=0.0;
   int i;
   max_value = x_buffer[0];
   min_value = x_buffer[0];
   for(i=1;i<number;i++)
   {
      if(x_buffer[i]>max_value)
         max_value=x_buffer[i];
      if(x_buffer[i]]<min_value) 
         min_value = x_buffer[0];
   }
   for(i=0; i<number;i++) 
      sample_value =sample_value + x_buffer[i];
   sample_value = (sample_value - max_value - min_value)/(number-2);
   return(sample_value);
}

函数假设存放连续进行5次采样的数据缓冲区data_buf[5],对于多于5次的防脉冲干扰平均值滤波,只需对该函数稍作修改即可,其中Sample_value表示本次采样的有效数据,number表示连续进行的采样次数。

6、低通数字滤波

低通滤波也称一阶滞后滤波,方法是第N次采样后滤波结果输出值是(1-a)乘第N次采样值加a乘上次滤波结果输出值。可见a<<1。该方法适用于变化过程比较慢的参数的滤波的C程序函数如下:

float low_filter(float low_buf[])
{
   float sample_value;
   float X=0.01;
   sample_value=(1-X)*low_buf[1]+X*low_buf[0];
   retrun(sample_value);
}