摘要

最近在做直流电机的毕设中，由于需要采集转速，电流，电压，温度等参数，常规的采集容易受到干扰，所以特意复习了一下关于数字滤波有关的知识，并作出相应的整理。本文对数字滤波进行简单介绍，讲解七种常用的滤波算法并用C语言实现，并比较其优缺点。由于篇幅原因未能在嵌入式系统中实验，读者可以自行验证。本篇篇幅较长，建议收藏。

所用工具：

1、测试软件：Visual C++6.0

知识概括：

通过本篇文章您将学到：

1、数字滤波算法基础知识

一、数字滤波简介

很多嵌入式系统都需要通过AD转换方式采集模拟信号，当干扰作用于模拟信号之后，AD转换结果就会偏离真实值。如果仅采样一次，是无法确定该结果是否可信的，必须多次采样，得到一个AD转换的数据序列，通过某种处理后，才能得到一个可信度较高的结果。这种用软件算法从采样数据序列中提取逼近真值数据的方法，通常称为数字滤波。数字滤波有硬件滤波的功效，却节省了硬件投资。实现数字滤波功能的软件算法称为数字滤波算法，由于数字滤波算法的灵活性，其效果往往是硬件滤波电路达不到的。
数字滤波的不足之处就是需要消耗一定的CPU机时，在进行实时信号处理时必须充分考虑到这一点。当有用信号为高频信号，干扰信号为低频信号时，需要采用高通滤波，而实时高通滤波算法收到CPU速度的限制，往往力不从心，而硬件高通滤波器却很容易处理这种实时信号。当有用信号为低频信号（如温度，湿度，重量，水位等），其频率通常很低（甚至接近0.001Hz），而干扰信号频率先对较高（如各种电磁干扰），需要采用低通滤波。由于硬件低通滤波体积往往较大，与当前电子产品微型化趋势不相容。而当前CPU的速度在执行实时低通滤波算法时完全可以胜任，故数字滤波主要应用领域为实时低通滤波。在某些情况下，数字滤波也用来进行带通滤波。
对于非实时信息处理，由于不受处理速度的制约，数字滤波算法广泛应用于各种音频数据，图片数据和说数据的加工处理。
如果用硬件手段（FPGA/CPLD）来实现数字滤波，其处理速度将有质的飞越，滤波算法采用数学工具软件MATLAB来设计，滤波功能也将更加丰富多彩，应用频率范围也大大扩展。以下介绍在嵌入式系统中广泛使用的几种数字滤波算法，这些算法都是用软件进行实时低频信号处理，达到提高系统抗干扰能力的目的。

二、常用数字滤波算法

1.限幅滤波

生活经验告诉我们，很多物理量的变化是需要一定时间的，相邻两次采样值之间的变化有一个合理的限度（即有用信号高频分量的最高频率有限）。例如在热处理车间的大型电炉里，工作的温度变化不可能在短时间（采样间隔）内发生剧烈变化，如果相邻两次采样值之间的变化未超过预定范围，说明该采样值未受到明显干扰，可以采用。如果相邻两次采样值之间的变化超过预定范围，说明该采样值受到了明显干扰（毛刺型突发干扰），不能采用。对于不可信得采样数据，必须输出一个合理的替代数据，以保证采样数据序列的连续性和完整性。在要求不高的场合，可以用前一个采样数据来替代这个受干扰的采样数据。
限幅滤波C语言代码如下：

#include "stdafx.h"

int adcData[10]={10,12,15,13,14,25,16,16,20,9};            //手动定义ADC采集的值

#define D 10                                               //相邻两次采样值之间最大允许变化值
int history;                                               //上次采样值

int main(int argc, char* argv[])
{
	int neo;                                               //定义本次采样值
	history=adcData[0];                                    //定义上次采样值初值
	for(int i=0;i<10;i++)                                  //循环访问数组adcData
	{
		neo=adcData[i];                                    //将数组里的元素作为本次采样值
		if(i>0)                                            //i=0时，无上次采样值，故判断
		{
			if((neo-history>D)||(history-neo>D))           //限幅判断
				neo=history;                               //如果超过最大变化值，则新值变为上次采样值
			history=neo;                                   //本次采样值变为上次采样值，开始下一轮循环
		}
		printf("%d\n",neo);                                //输出滤波后的值
 	}
	return 0;
}

//循环neo输出结果为：10,12,15,13,14,14,16,16,20,20
//可见第六位和第10位由于相差超过10，故被前一个值替换

当被检测物理量处于明显变化阶段时，相邻两个采样值之间的差别也比较明显，这时用前一个采样值来替代当前受干扰的采样值存在较大的误差。因此，可以利用被检测物理量的变化趋势来进行预测，选择一个更加合理的数据来替代当前受干扰的采样值，使得替代误差大大减小。
为了掌握被检测物理量当前变化趋势，至少需要两个以上的采样值历史记录。我们采用最简单的线性预测算法，只有两个历史采样值就可以了。线性预测算法的含义为“当前采样值的变化趋势与前一次采样值的变化趋势相同”，即采样值在短期内保持相同的上升趋势和下降趋势。由于代码类似，只是加一些限制条件，故这里省略。
限幅滤波的关键是合理确定采样周期和相邻两次采样值的最大允许变化值。采样周期可以通过检测精度（硬件是AD模块分辨率，软件是采样频率）和物理量最大变化速率（比如温度每秒钟对多上升多少摄氏度）来计算。另外，采样周期不能太短，必须大于干扰的持续时间，一面一次干扰造成两次以上的采样值不准确。
限幅滤波本质上是低通滤波，由于毛刺型突发干扰为高频干扰，故可以被很好地滤除。另一方面，限幅滤波采用限幅手段来达到滤波目的，对于幅度比较小的噪声干扰无能为力，即不能滤除随机干扰，因而只能用于对精度要求不是很高的场合。

2.中值滤波

限幅滤波有一个潜在的隐患，在连续两个以上采样周期受到强干扰后，系统可能会不稳定。因此，限幅滤波只能在基本没有干扰的场合下采样，他只能滤除极个别偶发的毛刺型干扰。在环境恶化，干扰频繁的情况下，可以采用中值滤波算法来处理。该算法描述是：连续进行奇数次采样，然后将采样得到的数据样本进行排序，取中间数据样本作为有效采样值。例如连续采样5次，得到5个采样数据样本（例如17、17、29、27、16），然后进行排序（16、17、17、27、29），最后取中间的（第三个，即17）采样数据样本作为有效采样值输出。由于受到干扰的采样值偏离有效采样值，排序后必然处于两端的位置，只要受到干扰的采样数据样本个数小于总采样数据样本数目的一半，就可以确保中值采样数据样本的有效性。
在采样中值滤波算法时，需要确定两个采样周期：一个是主采样周期（T0）和子采样周期（T1），每个主采样周期输出一个有效采样值，具体的时间长度由系统精度和有效信号的高频分量决定。每个子采样周期进行一次采样操作，它是由干扰频繁程度和连续采样次数N决定。

采用中值滤波算法必须满足以下条件：
被采样的物理量本身在连续N次子采样周期期间是基本稳定的，其变化量小于系统精度要求，可以忽略不计。
每次干扰的最长持续时间已知，相邻两次干扰的间隔时间虽然不定，但远大于一次干扰的持续时间。这样，可以使得连续采样的数据样本中被干扰的样本数目不会超过总样本数目N的一半，确保中值样本的有效性。
中值滤波C语言代码如下：

#include "stdafx.h"

int adcData[20]={10,12,15,13,14,25,16,16,20,9,               //手动定义ADC采集的值
				13,14,25,16,16,20,9,10,12,15}; 

#define N 5                                                  //定义N个数取中值

int value_buff[N];                                           //定义中值计算缓冲区

int main(int argc, char* argv[])
{
	int neo;                                                 //定义滤波计算完输出的值
	char temp;                                               //定义冒泡排序法临时变量
	for(int major_col=0;major_col<20;major_col+=N)           //定义主采样周期，每个包括N个子采样周期
	{
		for(int minor_col=0;minor_col<N;minor_col++)         //采集每个子采样周期的值
		{
			value_buff[minor_col]=adcData[minor_col+major_col];
		}
		//冒泡排序法
		for(int i=0;i<N-1;i++)
		{
			for(int j=0;j<N-i-1;j++)
			{
				if(value_buff[j]>value_buff[j+1])
				{
					temp=value_buff[j];
					value_buff[j]=value_buff[j+1];
					value_buff[j+1]=temp;
				}
			}
 		}
		neo=value_buff[(N-1)/2];                            //输出的值取缓冲区内的中值
		printf("%d\n",neo);
  	}
	return 0;
}

//循环neo输出结果为：13,16,16,12
//可见四个主采样周期的值都是子采样周期数据的中值

中值滤波本质上也是低通滤波，由于毛刺型突发干扰为高频干扰，故可以被很好地滤除。由于中值滤波进行N次采样才输出一次有效值（即采样输出比N：1），抗突发干扰能力比限幅滤波要提高很多，能够在干扰比较频繁的场合正常工作。另外，在N个数据样本中，中值样本通常受到的噪声干扰也比较小，故中值滤波算法也具备一定的抗随机干扰能力（其能力与N正相关）。
中值滤波的效果比限幅滤波提高很多，但必须满足主采样周期远大于子采样周期的条件。而子采样周期受到毛刺型突发干扰持续时间的限制不能太短，故主采样周期必然较长，即被检测物理量的基波频率一定比较低（慢变信号），其应用场合受到一定的限制。

3.算术平均滤波

由于随机干扰（噪声型干扰）随着数据样本的增加其统计平均值区域零，故对被检测物理量进行连续多次采样，然后求其算数平均值作为有效采样值，就可以达到抑制随机干扰的效果。连续采样次数越多，抑制随机干扰的效果越好。这种滤波算法就是算术平均滤波算法。
算术平均滤波C语言代码如下：

#include "stdafx.h"

float adcData[20]={10,12,15,13,14,25,16,16,20,9,               //手动定义ADC采集的值
				13,14,25,16,16,20,9,10,12,15}; 

#define N 5												   	//定义N个数取算术平均

float value_buff[N];                                           //定义算术平均计算缓冲区

int main(int argc, char* argv[])
{
	float neo=0.0;                                                 //定义滤波计算完输出的值
	for(int major_col=0;major_col<20;major_col+=N)                 //定义主采样周期，每个包括N个子采样周期
	{
		for(int minor_col=0;minor_col<N;minor_col++)               //采集每个子采样周期的值
		{
			value_buff[minor_col]=adcData[minor_col+major_col];
			neo+=value_buff[minor_col];
		}
		neo=neo/N;                                                 //输出算数平均
		printf("%0.2f\n",neo);
	}
	return 0;
}

//循环neo输出结果为：12.80 19.76 20.75 17.35
//可见四个主采样周期的值都是子采样周期数据的算数平均

算术平均滤波对随机干扰的抑制效果与连续采样次数N密切相关，当采样次数N增加到一定程度后，被测物理量本身产生的变化量就会超过允许范围，再也不能忽略不计了，因此连续采样次数N是不能任意增加的。其应用场合与中值滤波的应用场合类似，必须满足主采样周期远大于子采样周期的条件。而子采样周期受到毛刺型突发干扰持续时间的限制不能太短，故主采样周期必然较长，即被检测物理量的基波频率一定比较低（慢变信号），其应用场合受到一定限制。
虽然算术平均滤波对随机干扰的抑制效果比较好，但却不能消除毛刺型突发干扰。只要有一个采样数据样本受到毛刺型突发干扰，算术平均值将明显偏离真实值，所以一般需要结合其他限制条件去对算术平均滤波算法进行优化。

4.去极值平均滤波

算术平均滤波不能消除毛刺型突发干扰，只是通过平均操作将其影响削弱。因毛刺型突发干扰使采样值远离真实值，针对毛刺型突发干扰这一特点，可以比较容易地将其剔除，不参加平均值计算，从而使平均滤波的输出值更接近真实值。算法原理如下：连续采集N次，将其累加求和，同时找出其中的最大值和最小值，再从累加和中减去最大值和最小值，将剩余的N-2个采样值求平均，即得有效采样值。这种滤波算法就是去极值平均滤波算法。
算术平均滤波C语言代码如下：

#include "stdafx.h"

float adcData[20]={10,12,15,13,14,25,16,16,20,9,               //手动定义ADC采集的值
				13,14,25,16,16,20,9,10,12,15}; 

#define N 5												   	//定义N个数取算术平均

float value_buff[N];                                           //定义算术平均计算缓冲区
float nvalue_buff[N-2];                                        //定义去极值之后的缓冲区
int main(int argc, char* argv[])
{
	float neo=0.0;                                                 //定义滤波计算完输出的值
	float temp;
	for(int major_col=0;major_col<20;major_col+=N)           //定义主采样周期，每个包括N个子采样周期
	{
		for(int minor_col=0;minor_col<N;minor_col++)         //采集每个子采样周期的值
		{
			value_buff[minor_col]=adcData[minor_col+major_col];
		}
		//冒泡排序法
		for(int i=0;i<N-1;i++)
		{
			for(int j=0;j<N-i-1;j++)
			{
				if(value_buff[j]>value_buff[j+1])
				{
					temp=value_buff[j];
					value_buff[j]=value_buff[j+1];
					value_buff[j+1]=temp;
				}
			}
		}
		for(int k=0;k<N-2;k++)                                //去极值并且放入缓冲区
		{
			nvalue_buff[k]=value_buff[k+1];
			neo+=nvalue_buff[k];
		}
		neo=neo/(N-2);                                        //输出值为去极值之后的算数平均值
		printf("%0.2f\n",neo);
	}
	return 0;
}

//循环neo输出结果为：13.00 21.67 22.56 19.85
//可见四个主采样周期的值都是子采样周期数据去极值后的算数平均

去极值平均滤波同时具有消除毛刺型突发干扰和噪声型随机干扰的能力，在低频信号采集系统中被广泛运用。比如在各种文艺表演大赛的评奖计分中经常听到主持人所说的“去掉一个最高分，去掉一个最低分，最终的得分是xxx”，这实际上就是采用了去极值平均滤波算法。
如果连续N次采样中没有受到毛刺型干扰的数据样本，被剔除的将是两个随机误差最大的数据样本，滤波效果仍然很好。去极值平均滤波算法也有效果不佳的时候，如果连续N次采样中有两个以上的正毛刺数据样本或者两个以上的负毛刺数据样本，将至少有一个毛刺数据样本不能剔除，最终平均值必然与真实值存在明显偏差。
去极值平均滤波的应用条件与算术平均滤波的应用条件相同，主要用来对低频信号进行软件滤波。

5.滑动平均滤波

前面介绍的两种平均滤波算法有一个共同点，即每个主采样周期取得的一个有效采样值是由连续若干个子采样值计算出来的，为此，这些算法的主采样周期都是比较长的。当被检测物理量变化快时，较长的主采样周期不能及时采集到被检测物理量中的高频信息，致使系统反应迟钝，实时性难以保证。要保证系统的实时性，必须要按技术指标的要求相应缩短主采样周期，而子采样周期不能随意缩短（受干扰持续时间的制约），只好减少连续采样次数，连续采样次数的减少又直接降低了滤波效果。既然如此，干脆直接按技术指标的要求确定一个合理的采样周期，每个采样周期只进行一次采样操作（即采样输出比1:1），然后将当前的采样值和之前连续的若干次采样值（最近历史采样值）一起进行平均，将得到的平均值作为当前有效采样值投入使用。由于参与平均运算的历史采样值个数固定且内容不断更新，相当于一个滑动的时间窗口，故将这种平均滤波方式称为“滑动平均滤波”，窗口的大小用采样数据样本个数m来表示。按这种方式进行采样就不再存在主采样周期和子采样周期的问题。

滑动平均滤波C语言代码如下：

#include "stdafx.h"

float adcData[20]={10,12,15,13,14,25,16,16,20,9,               //手动定义ADC采集的值
				13,14,25,16,16,20,9,10,12,15}; 

#define N 5												   	//定义N个数取算术平均

float value_buff[N];                                           //定义算术平均计算缓冲区

int main(int argc, char* argv[])
{
	float neo=0.0;                                           //定义滤波计算完输出的值
	for(int step=0;step<20-N;step++)                         //缓冲区的循环
	{
		for(int count=0;count<N;count++)                     //缓冲区内N个数的赋值
		{
			value_buff[count]=adcData[count+step];
			neo+=value_buff[count];
		}
		neo=neo/N;                                           //输出滑动平均后的值，一共20-N个
		printf("%0.2f\n",neo);
	}

	return 0;
}

//循环neo输出结果为：12.80 18.36 20.27 20.85 22.37 21.67 19.13 18.23
//                   19.85 19.37 20.67 22.33 21.67 18.53 17.11
//可见添加滑动窗之后会更加平滑
//理应输入输出是1:1，这里为了简便代码，故这样输出，只要知道思想就可以了，具体算法具体使用

6.滑动加权滤波

在滑动平均滤波算法中，窗口中的m个采样数据样本以平等身份参与运算，这对抑制随机干扰比较有利，但将有超时严重的滞后效果，降低系统反应速度。为了提高系统的反应速度，滤波算法的输出应该及时反映当前采样值中包含的有效信息，即增加即时数据样本和近期数据样本的权重，降低早期数据样本的权重。为此，为滑动窗口中的各个数据样本分配不同的“加权系数”，进行加权平均运算。这种滤波算法称为“滑动加权滤波”。

滑动加权滤波C语言代码如下：

#include "stdafx.h"

float adcData[20]={10,12,15,13,14,25,16,16,20,9,               //手动定义ADC采集的值
				13,14,25,16,16,20,9,10,12,15}; 

#define N 5												   	//定义N个数取算术平均

float value_buff[N];                                           //定义算术平均计算缓冲区
float perc_buff[N]={0.1,0.2,0.1,0.3,0.3};                  //定义缓冲区权重


int main(int argc, char* argv[])
{
	for(int step=0;step<20-N;step++)                         //缓冲区的循环
	{
		float neo=0.0;                                                 //定义滤波计算完输出的值
		for(int count=0;count<N;count++)                     //缓冲区内N个数的赋值
		{
			value_buff[count]=adcData[count+step];    
			value_buff[count]*=perc_buff[count];              //每个值乘对应的权重
			neo+=value_buff[count];							//累加
		}
		printf("%0.2f\n",neo);
	}
	return 0;
}

//循环neo输出结果为：13.00 17.20 17.80 16.20 18.80 16.00 13.40 14.60
//                   16.80 17.20 16.20 18.80 16.00 12.50 13.10
//可见添加滑动窗之后会更加平滑

7.一阶滞后滤波

将一阶滞后滤波器的微分方程用差分方程来表示，便可以用软件算法来模拟硬件一阶滞后滤波器（RC低通滤波器）的功能。

由上式可以看出，本次滤波的输出值主要取决于上次滤波的输出值，本次采样值对本次滤波输出的贡献是比较小的，但多少有些修正作用。这种算法便模拟了具有较大惯性的一阶滞后滤波器的功能。
当目标参数为变化很慢的物理量时（如大型储水池的水位信号），采用一阶滞后滤波算法是很有效的。
另一方面，他不能滤除高于二分之一采样频率的干扰信号。一阶滞后滤波算法程序流程与加权平均滤波相似，而加权系数只有两个：a和（1-a）。
一阶滞后滤波C语言代码如下：

#include "stdafx.h"

float adcData[10]={10,12,15,13,14,25,16,16,20,9};            //手动定义ADC采集的值

#define D 10                                               //相邻两次采样值之间最大允许变化值
#define A 0.85                                             //滤波系数
float history=0.0;                                               //上次采样值

int main(int argc, char* argv[])
{
	float neo=0.0;                                               //定义本次采样值
	history=adcData[0];                                    //定义上次采样值初值
	for(int i=0;i<10;i++)                                  //循环访问数组adcData
	{
		neo=adcData[i];                                    //将数组里的元素作为本次采样值
		neo=A*neo+(1-A)*history;                           //滤波加权
		printf("%0.2f\n",neo);                                //输出滤波后的值
 	}
	return 0;
}

//循环neo输出结果为：10.00 11.70 14.25 12.55 13.40 22.75 15.10 15.10 18.50 9.15
//可见滤波系数的取值对滤波效果尤为重要

在设置滤波系数时需要同时兼顾系统的平稳性和灵敏性，而两者在一阶滞后滤波算法中属于不能同事达到最优，只能取得相对最优解，所以找到一个平衡点很重要，这也需要具体工程具体应用。当数据变化较为快速时需要以灵敏性优先考虑，当数据趋于稳定时需要以平稳性优先考虑。

三、数字滤波小结

下面将从采样输出比，康毛刺干扰能力和抗噪声干扰能力三个方面进行对比，如下表所示：

每种数字滤波算法都有自己的特点和应用条件，应该根据实际应用场合来合理选择数字滤波算法，并根据技术参数要求来确定采样周期。以上这些也只算是比较基础的滤波方式，可以相互结合使用，但也要考虑采集数据的类型。读者掌握这些之后可以在学习一些更为升入的滤波方式，比如卡尔曼，维纳等，技多不压身。

PS：最最最重要的，码字不易，求各位大佬看完点个赞，感谢！！！

三种常见的模糊图像滤波算法：均值滤波、高斯滤波、中值滤波。

均值滤波

用滤波掩膜确定的领域内像素的平均灰度值去代替图像的每个像素点，目的是为了“减噪”。

缺点：存在不希望的边缘模糊效应。

高斯滤波

高斯滤波是一种线性平滑滤波，适用于消除高斯噪声，高斯滤波器是一类根据高斯函数的形状来选择权值的线性平滑滤波器。在对图像进行卷积的时候其实是对图像进行点操作，用所求像素点周围的像素通过模板的权值计算得到新的像素值。而图像形状的不同也就代表了不同的参数，从而得到不同的模板。

二维高斯分布：

A：表示高斯分布的峰值，（,）表示均值，表示方差，其中假设了，即在x轴与y轴方向的方差相同。

 当且仅当时，二维高斯函数为正太分布随机变量的概率密度函数，高斯函数的曲面下的积分为1：

为什么需要二维高斯函数的积分为1呢？

我的理解是，我们是希望我们基于高斯函数设计的高斯核不会使得卷积之后的像素值超过255，假如说一个3*3的卷积核用来处理像素值均为255的一个原图像的区域，那么得到的中心点的像素值为255*(卷积核的和)，只有当卷积核的和为1时，此时中心点的像素值不会超过255。

我们设计高斯核是基于高斯函数设计的，也就是说我们的高斯核的形状与高斯函数保持一致。因此高斯核有如下特点：

1. （x,y）均值为卷积核的中心点坐标（kSize-1）/2，kSize为卷积核的大小。比如3*3的卷积核，其中心点坐标为（1,1），图像坐标从0开始。
2. 卷积核的方差根据卷积核的大小kSize计算，在Opencv中的GaussianBlur函数原型是这么设计方差的：σ=0.3((ksize−1)×0.5−1)+0.8。我们可以通过方差来计算卷积核大小，也可以通过卷积核大小来计算方差。
3. 需要归一化处理，使得高斯核的和为1。

参考博文：高斯滤波及其原理_声希Censh的博客-CSDN博客_高斯滤波

参考博文中给出了自己创建高斯卷积核的C的源代码：

代码中用的高斯函数是：，所以在归一化之前，卷积核中心点的值为1。卷积核采用的归一化方法为：高斯核各点的值/sum，sum为高斯核的和。

#include <stdio.h>
#include <math.h>
double** createGuassion(int, double);
double** createGuassion(int kSize, double sigma)
{
	if (kSize % 2 == 0 || kSize < 0 || (kSize == 0 && sigma == 0))
	{
		return 0;//当模板大小为偶数、或者小于0、或者模板大小与sigma同时为0时返回
	}
	else if (kSize == 0)
	{
		kSize = 2 * (sigma - 0.8) / 0.3 + 3;
		kSize = round(kSize);//
		if (kSize % 2 == 0)
		{
			kSize += 1;
		}
	}
	else if (sigma == 0)
	{
		sigma = 0.3 * ((kSize - 1) * 0.5 - 1) + 0.8;//sigma缺省时，利用该公式计算sigma
	}

	double** guass;//生成的高斯矩阵
	double sum = 0;//用于归一化求和
	double x2 = 0;
	double y2 = 0;
	int center = (kSize - 1) / 2;//高斯核中心点的行标和列标（从0开始计），所以必须保证kSize为奇数
	guass = new double* [kSize];//注意，double*[k]表示一个有10个元素的指针数组
	for (int i = 0; i < kSize; ++i)
	{
		//先创建一个kSize×kSize的模板矩阵
		guass[i] = new double[kSize];
	}
	for (int i = 0; i < kSize; i++)
	{
		x2 = pow(double(i - center), 2); //高斯函数中的x平方
		for (int j = 0; j < kSize; j++)
		{
			y2 = pow(double(j - center), 2);//高斯函数中的y平方
			sum += guass[i][j] = exp(-(x2 + y2) / (2 * sigma * sigma));//赋值并累加
		}
	}
	//归一化处理
	if (sum != 0)
	{
		for (int i = 0; i < kSize; i++)
		{
			for (int j = 0; j < kSize; j++)
			{
				guass[i][j] /= sum;
			}
		}
	}
	return guass;
}

int main()
{
	double** p = createGuassion(3, 0);
	for (int i = 0; i < 3; ++i)
	{
		for (int j = 0; j < 3; ++j)
		{
			printf("%-9f\x20", *(*(p + i) + j));
			//p+i表示a[i]的地址
			//"%-9f中, '-'表示左对齐(无'-'默认为右对齐)；9表示这个元素输出时占两个空格的空间
			//"\x20"表示以十六进制数20所代表的字符，即空格（space的ASCII为32）
		}
		printf("\n");
	}
	printf("\n");
	for (int i = 0; i < 3; ++i)
	{
		for (int j = 0; j < 3; ++j)
		{
			printf("%7.4f\x20", *(*(p + i) + j));
			//p+i表示a[i]的地址
			//"%-7.4f中, '-'表示左对齐(无'-'默认为右对齐)；
			//7.4表示这个元素输出时占7列，其中包括4个小数
		}
		printf("\n");
	}
	printf("\n");
	return 0;
}

运行结果如下：

0.057118  0.124758  0.057118
0.124758  0.272496  0.124758
0.057118  0.124758  0.057118

 0.0571  0.1248  0.0571
 0.1248  0.2725  0.1248
 0.0571  0.1248  0.0571

中值滤波 

中值滤波对于椒盐噪声的处理特别好，它算是统计排序滤波的一种，中值指的就是领域内像素的统计中值。

数字图像处理中对于中值滤波的描述是：用n*n的中值滤波器去除那些相对于其领域像素更亮或更暗，并且其区域小于（滤波器区域的一半）的孤立像素集。.....而对较大的像素影响较小。

这段话理解起来是这样：

1. 中值滤波只对其亮度明显大于或者小于其他周围像素点的噪声有作用。
2. 中值滤波器的大小必须是噪声点尺寸的2倍以上。因为中值定大于或小于一半的像素值。

而椒盐噪声非常符合第1点，所以中值滤波对于椒盐噪声的作用非常明显。

中值滤波与均值滤波对于椒盐噪声滤除作用的对比代码：

#include <iostream>
#include <opencv2/core.hpp>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>
#include <opencv2/opencv.hpp>

using namespace cv;
using namespace std;

//给图像加入椒盐噪声
void salt(Mat& image, int n)
{
	for (int k = 0; k < n; k++)
	{
		int i = rand() % image.cols;
		int j = rand() % image.rows;

		if (image.channels() == 1)
		{   //判断是一个通道
			image.at<uchar>(j, i) = 255;
		}
		else {
			image.at<cv::Vec3b>(j, i)[0] = 255;
			image.at<cv::Vec3b>(j, i)[1] = 255;
			image.at<cv::Vec3b>(j, i)[2] = 255;
		}
	}
}

int main()
{

	Mat src = Mat::zeros(Size(10, 10), CV_8U);
	salt(src, 20);
	cout << "椒盐噪声污染之后：\n" << src << endl;
	//imshow("椒盐噪声污染之后", src);

	Mat median_blur;
	medianBlur(src, median_blur, 3);
	cout << "中值滤波之后：\n" << median_blur << endl;
	//imshow("中值滤波之后", median_blur);

	Mat uniform_blur;
	blur(src, uniform_blur, Size(3, 3));
	cout << "均值滤波之后：\n" << uniform_blur << endl;
	//imshow("均值滤波之后", uniform_blur);

	//waitKey(0);

	return 0;
}

程序运行结果如下：

椒盐噪声污染之后：
[  0,   0,   0,   0, 255,   0,   0,   0,   0,   0;
   0, 255,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0, 255,   0,   0, 255,   0,   0,   0, 255;
   0,   0, 255,   0,   0,   0,   0,   0,   0,   0;
   0, 255, 255,   0,   0, 255,   0,   0,   0, 255;
   0,   0,   0,   0,   0, 255,   0,   0, 255,   0;
   0, 255,   0,   0,   0,   0,   0, 255,   0,   0;
   0, 255,   0,   0,   0,   0,   0,   0,   0,   0;
   0, 255,   0,   0,   0,   0,   0,   0, 255,   0;
   0,   0,   0,   0,   0,   0,   0, 255,   0,   0]
中值滤波之后：
[  0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0;
   0,   0,   0,   0,   0,   0,   0,   0,   0,   0]
均值滤波之后：
[113,  57,  57,  28,  28,  28,   0,   0,   0,   0;
  57,  57,  57,  57,  57,  57,  28,   0,  28,  28;
  57,  85,  85,  57,  28,  28,  28,   0,  28,  28;
  57, 113, 113,  85,  57,  57,  57,   0,  57,  57;
  57,  85,  85,  57,  57,  57,  57,  28,  57,  85;
 113,  85,  85,  28,  57,  57,  85,  57,  85,  85;
 113,  57,  57,   0,  28,  28,  57,  57,  57,  57;
 170,  85,  85,   0,   0,   0,  28,  57,  57,  57;
 113,  57,  57,   0,   0,   0,  28,  57,  57,  57;
 113,  57,  57,   0,   0,   0,  28,  85,  85, 113]

当然可以明显看到中值滤波比均值滤波的效果好不少。对于实际图像的过滤效果如下：

（a）原图 （b）椒盐噪声污染之后

（c）均值滤波效果

（d）中值滤波效果