数字滤波器设计之一：巴特沃斯（Butterworth）滤波器

zhwzhwei

67518人浏览 · 2017-04-30 23:25:01

zhwzhwei · 2017-04-30 23:25:01 发布

我们常说的经典滤波器是根据傅里叶分析和变换设计出来的，只允许一定频率范围内的信号成分正常通过，而阻止另一部分频率成分通过。按照最佳逼近特性或者滤波通带特性分类，主要为巴特沃斯滤波器（Butterworth）、切比雪夫滤波器（Chebyshev）、贝塞尔滤波器（Bessel）和椭圆滤波器（Elliptic）四种。每种MATLAB都有相应的函数，用起来也比较方便，但是却缺少C/C++的程序，于是自己仔细研究了每种滤波器的特性和原理，并且部分滤波器实现了C语言的代码化，接下来的时间会对这些滤波器的原理和C语言的实现进行介绍。

该系列均以低通滤波器为原型来介绍，其他类型的滤波器可以由低通滤波器通过频率变换转换得到，这里不过多介绍。低通滤波器的主要性能指标有以下几个：通带截止频率fp、阻带截止频率fs、通带衰减 ( Ap)、阻带衰减 ( As) 以及归一化频率时需要用到的-3dB的转折频率fc。

1. Butterworth滤波器原理

Butterworth滤波器因其在通带内的幅值特性具有最大平坦的特性而闻名，是四种经典滤波器中最简单的，巴特沃斯滤波器只需要两个参数表征，滤波器的阶数N和-3dB处的截止频率 $\Omega _c$ 。其幅度平方函数为：

$\bg_white \left | \mathrm{H}_a\left ( j\Omega \right ) \right |^{2}=\frac{1}{1+\left ( \frac{\Omega }{\Omega _c} \right )^{2N}}$

$\bg_white A_{p}=10lg\left ( 1+\varepsilon _{p}^{2} \right ),A_{s}=10lg\left ( 1+\varepsilon _{s}^{2} \right ),A\left ( \Omega \right )=10lg\left ( 1+\left ( \frac{\Omega }{\Omega _c} \right )^{2N} \right )$

N是滤波器的阶数，从幅度平方函数可以看出，N阶滤波器有2N个极点，而且这2N个极点均布在一个圆上，圆的半径为，称之为Butterworth圆，Butterworth滤波器系统是一个线性系统，要使其稳定，其极点必须位于S平面的左半平面，所以取左半平面内的N个极点作为滤波器的极点，滤波器就是稳定的了，求出极点之后，计算模拟滤波器的系数as、bs，然后通过双线性变换（不懂得自行查书）由模拟域到数字域，求出系数az和bz 。最后通过差分方程就可以计算滤波结果了。

$\bg_white y(n)= \sum_{k=0}^{\propto }h(k)x(n-k)=-\sum_{k=1}^{M}a_{k}y(n-k)+\sum_{k=0}^{N}b_{k}x(n-k)$

2. C语言实现

A．求阶数

公式为：

$N=\frac{1}{2}\frac{lg\left ( \frac{10^{0.1A_s}-1}{10^{0.1A_p}-1} \right )}{lg\left ( \frac{\Omega _s}{\Omega _p} \right )}$

代码：

N = ceil(0.5*( log10 (( pow (10, Stopband_attenuation/10) - 1)/
	   ( pow (10, Passband_attenuation/10) - 1)) / log10 (Stopband/Passband) ));

B．求极点

公式：

代码：

for(k = 0;k <= ((2*N)-1) ; k++)
    {
		if(Cutoff*cos((2*k+N-1)*(pi/(2*N))) < 0.0)
         {	   
            poles[count].x = -Cutoff*cos((2*k+N-1)*(pi/(2*N)));
			poles[count].y = -Cutoff*sin((2*k+N-1)*(pi/(2*N)));	  
			 count++;
	        if (count == N) break;
         }
    }

C．计算模拟滤波器系数

公式：

代码（这部分需要自己推导，多算几步就找到规律了）：

     Res[0].x = poles[0].x; 
     Res[0].y = poles[0].y;

     Res[1].x = 1; 
     Res[1].y= 0;

    for(count_1 = 0;count_1 < N-1;count_1++)//N个极点相乘次数
    {
	     for(count = 0;count <= count_1 + 2;count++)
	     {
	          if(0 == count)
			  {
 	              Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] );  
			  }

	          else if((count_1 + 2) == count)
	          {
	                Res_Save[count].x += Res[count - 1].x;
					Res_Save[count].y += Res[count - 1].y;	
	          }		  
		    else 
		    {
				Res_Save[count] = ComplexMul( Res[count], poles[count_1+1] );	               				
			    Res_Save[count].x += Res[count - 1].x;
			    Res_Save[count].y += Res[count - 1].y;	
		    }
	     }

	     for(count = 0;count <= N;count++)//Res[i]=a(i),i越大次数越高
	     {
			Res[count].x = Res_Save[count].x; 
			Res[count].y = Res_Save[count].y;
				 
		    *(a + N - count) = Res[count].x;
	     }				 
	}

     *(b+N) = *(a+N);

D．双线性变换

用下式进行替换 $H_a\left ( s \right )$ 中的变量，得到 $H\left ( \textrm{z} \right )$ ，然后类似上面的计算方法计算Z域的系数az和bz，其中T为采样周期，但是因为在计算中会被约去，所以简化计算，这里取1。

公式：

代码：

int Count = 0,Count_1 = 0,Count_2 = 0,Count_Z = 0;
      double *Res, *Res_Save;
	  Res = new double[N+1]();
	  Res_Save = new double[N+1](); 
	  memset(Res, 0, sizeof(double)*(N+1));
	  memset(Res_Save, 0, sizeof(double)*(N+1));

      	for(Count_Z = 0;Count_Z <= N;Count_Z++)
		{
            *(az+Count_Z)  = 0;
		    *(bz+Count_Z)  = 0;
		}

	
	for(Count = 0;Count<=N;Count++)
	{    	
	      for(Count_Z = 0;Count_Z <= N;Count_Z++)
	      {
	      	    Res[Count_Z] = 0;
				Res_Save[Count_Z] = 0;	 
	      }
          Res_Save [0] = 1;
	
	      for(Count_1 = 0; Count_1 < N-Count;Count_1++)//计算（1-Z^-1）^N-Count的系数,
	      {												//Res_Save[]=Z^-1多项式的系数，从常数项开始
				for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
				{
					 if(Count_2 == 0)  
					 {
					   Res[Count_2] += Res_Save[Count_2];
					 } 	

					 else if((Count_2 == (Count_1+1))&&(Count_1 != 0))  
					 {
						 Res[Count_2] += -Res_Save[Count_2 - 1];   
					 } 

					 else  
					 {
						 Res[Count_2] += Res_Save[Count_2] - Res_Save[Count_2 - 1];
					 }				 
				}

		      for(Count_Z = 0;Count_Z<= N;Count_Z++)
		      {
 				   Res_Save[Count_Z]  =  Res[Count_Z] ;
				   Res[Count_Z]  = 0;   
		      }	
	      }

		for(Count_1 = (N-Count); Count_1 < N;Count_1++)//计算(1-Z^-1)^N-Count*（1+Z^-1）^Count的系数,
	    {												//Res_Save[]=Z^-1多项式的系数，从常数项开始
                for(Count_2 = 0; Count_2 <= Count_1+1;Count_2++)
				{
	      			 if(Count_2 == 0)  
					 {
					   Res[Count_2] += Res_Save[Count_2];
					 } 	

					 else if((Count_2 == (Count_1+1))&&(Count_1 != 0))  
					 {
						 Res[Count_2] += Res_Save[Count_2 - 1];
					  } 

					 else  
					 {
						  Res[Count_2] += Res_Save[Count_2] + Res_Save[Count_2 - 1];
					 }				 
				}

		      for(Count_Z = 0;Count_Z<= N;Count_Z++)
		      {
 		          Res_Save[Count_Z]  =  Res[Count_Z] ;
			      Res[Count_Z]  = 0;    
		      }
	     }

		for(Count_Z = 0;Count_Z<= N;Count_Z++)
		{
            *(az+Count_Z) +=  pow(2,N-Count)  *  (*(as+Count)) * Res_Save[Count_Z];
			*(bz+Count_Z) +=  (*(bs+Count)) * Res_Save[Count_Z];		       
		 }	

	}//最外层for循环
   
     for(Count_Z = N;Count_Z >= 0;Count_Z--)
     {
          *(bz+Count_Z) =  (*(bz+Count_Z))/(*(az+0));
          *(az+Count_Z) =  (*(az+Count_Z))/(*(az+0));
     }

E．由由差分方程计算滤波结果

采用滤波器直接II型结构，可以减少一半的中间缓存内存，具体代码如下：

double FiltButter(double *pdAz,	//滤波器参数表1
                  double *pdBz,	//滤波器参数表2
				  int	nABLen,	//参数序列的长度
				  double dDataIn,//输入数据
				  double *pdBuf)	//数据缓冲区
{
	int	i;
	int	nALen;
	int nBLen;
	int	nBufLen;
	double	dOut;

	if( nABLen<1 )return 0.0;
	//根据参数,自动求取序列有效长度
	nALen = nABLen;
	for( i=nABLen-1; i; --i )
	{
		if( *(pdAz+i) != 0.0 )//从最后一个系数判断是否为0
		{
			nALen = i+1;	
			break;
		}
	}
	//printf("%lf ", nALen);
	if( i==0 ) nALen = 0;

	nBLen = nABLen;
	for( i=nABLen-1; i; --i )
	{
		if( *(pdBz+i) != 0.0 )
		{
			nBLen = i+1;
			
			break;
		}
	}
	//printf("%lf ", nBLen);
	if( i==0 ) nBLen = 0;
	//计算缓冲区有效长度
	nBufLen = nALen;
	if( nALen < nBLen)
		nBufLen = nBLen;

	//滤波: 与系数a卷乘
	dOut = ( *pdAz ) * dDataIn;  // a(0) * x(i)   

	for( i=1; i<nALen; i++)	// a(i) * w(n-i),i=1toN
	{
		dOut -= *(pdAz+i) * *(pdBuf + (nBufLen - 1) - i);
	} 

	//卷乘结果保存为缓冲序列的最后一个
	*(pdBuf + nBufLen - 1) = dOut;
	
	//滤波: 与系数b卷乘
	dOut = 0.0;
	for( i=0; i<nBLen; i++)	// b(i) * w(n-i)
	{    	
		dOut += *(pdBz+i) * *(pdBuf + (nBufLen - 1) - i);
	}

	//丢弃缓冲序列中最早的一个数, 最后一个数清零
	for( i=0; i<nBufLen-1; i++)
	{
		*(pdBuf + i) = *(pdBuf + i + 1);
	}
	*(pdBuf + nBufLen - 1) = 0;
	
	//返回输出值
	return dOut; 
}

VC6.0开发环境滤波结果如下：

至此就完成了Butterworth滤波器的设计过程，完整代码请自行下载。

参考资料：http://blog.csdn.net/zhoufan900428/article/details/9069475

源码地址：http://download.csdn.net/detail/zhwzhaowei/9830114

GitCode 开源社区

旨在为数千万中国开发者提供一个无缝且高效的云端环境，以支持学习、使用和贡献开源项目。

更多推荐

[转载]在Windows环境下安装GNU Radio

转自：在Windows环境下安装GNURadio_恐弱智_新浪博客GNU Radio是用Python开发的，大部分开源的工程能够在Linux环境下运行良好，而Windows下却运行的很勉强，而且安装配置都很复杂。GNU Radio算是个例外了，不光提供了Windows的二进制安装，还有比较详细的说明。我是Python小白，所以折腾了好久才弄好，特意记录下来，免得以后再装还折腾。GNU Radio的

GitCode 开源社区

centOS 8 使用dnf安装Docker

DNF是什么？CentOS 8使用YUM软件包管理器版本v4.0.4。现在，该版本使用DNF(已删除YUM)。DNF是软件包管理器。它会在Linux发行版上安装，执行更新并删除软件包。使用DNF安装Docker跳过具有损坏依赖性的程序包一个有效的解决方案是使您的CentOS 8系统使用以下--nobest命令安装最符合条件的版本：sudo dnf install docker...

GitCode 开源社区

定时同步数据库表(mysql+linux+crontab)

sync.sh里面的参数需要改变，ip/username/password/database/tablesync.sh#!/bin/sh# Please change the IP and password of the data source db.# Then change the table name.filename=/home/nington/db/$(date +%Y-%m