数字滤波器的设计和分析方法

在使用Matlab设计了一个滤波器后，可以通过两个矢量来得到离散的差分方程：B（分子系数）和A（分母系数）。在输入这两个矢量的情况下，Matlab可以用工具箱的很多函数快速地分析和画出滤波器的性能。通过help signal可以查阅Matlab信号处理工具箱的帮助说明。输入这个命令后的结果如下：

我们特别感兴趣的是频率响应、脉冲响应、零极点图和群延时。

创建一个单位脉冲响应图

和FIR滤波器不同，IIR滤波器的滤波器系数不是组成单位脉冲响应的每一项。IIR的单位脉冲响应必须迭代地进行计算。Matlab的impz命令可以大大简化这个过程。例如，如果我们使用butter命令（对于巴特沃斯滤波器）来设计一个滤波器并希望查看它的单位脉冲响应，我们需要如下命令：

这样就可以创建一个单位脉冲响应的图形。在第1行代码中，我们设计了一个4阶的巴特沃斯低通滤波器，它的截止频率为0.25F_S/2。在第2行代码中，我们使用impz命令来确定单位脉冲响应。其中，“B”和“A”分别为分子和分母的系数矢量，10表示单位脉冲响应计算和作图涉及的点数，48000是我们需要的滤波器采样频率。因为我们需要包含采样频率（第4个输入），结果图形中水平轴的坐标单位是时间。

除了使用butter命令，我们还可以使用SPTool（Signal Processing Tool，信号处理工具）命令。它调出一个用于设计数字滤波器的图形用户界面（Graphical User Interface，GUI）和一些其他的输出。当你从SPTool输出一个文件到Matlab工作环境中时，用来保存滤波器设计结果默认的结构的文件叫做filt1。在这种情况下，不是用上面的“B”和“A”矢量来表示滤波器的分子和分母的系数，而是使用“filt1.tf.num”和“filt1.tf.den”。图4.7所示为filt1的数据结构。

图4.7 输出到Matlab工作空间的Filt1的数据结构

另外，impz命令的第三个输入如下（假设使用SPTool）：

因为没有说明采样频率，结果的绘图中横坐标没有显示采样的单位（采样数n）。

最后，这里有一个二输入的impz命令的版本，它显示如下。如果你使用的是这个版本的impz命令格式如下：

这里的横坐标是采样数，但是这个算法自动估计采样的数目并画图。和这个简单版本的impz命令相关的图形如图4.8所示。

创建一个频率响应图：

通过使用一个Matlab命令freqz来画出频率响应图。使用这个命令的例子如下：

图4.9所示为结果频率响应曲线。如果在freqz命令中没有使用额外的参数来说明采样频率，图中横轴的频率归一化为×π的单位弧度频率。这等效归一化到F_S/2的频率。如果归一化的坐标等于0.25，它对应的频率为f=0.25F_S/2。这是截止频率（滤波器的幅频响应下降3dB的频率数值），我们在设计滤波器的时候就指定过这个数值。要注意IIR的非线性的相位响应。

图4.8 与一个4阶巴特沃斯低通滤波器相关的单位脉冲响应曲线，其截止频率为0.25F_S/2

图4.9 与一个4阶巴特沃斯低通滤波器相关的频率响应曲线， 它的截止频率为0.25F_S/2

创建一个零极点图：

在z平面画出零点和极点的位置要使用Matlab的zplane命令，这个命令的举例如下：

图4.10所示为这个零极点图。Matlab命令zplane调用其他的函数如zplaneplot。这个函数实际负责建一个包含单位圆的z平面的零极点图。不幸的是，尽管作者反馈意见到美国MathWorks公司，Matlab软件中仍然只使用70个点来形成单位圆。因为Matlab常常使用直线来连接各个点，从而导致形成69个多边形。所以在你的零极点非常接近于单位圆的时候，这样的情况就非常受限。有两个方法来解决这个问题。(www.xing528.com)

图4.10 与一个4阶巴特沃斯低通滤波器相关的零极点图，它的截止频率为0.25F_S/2

●推荐的方法：使用CD中包含的ucf函数来修正这个问题。一般地，ucf用白色线重写的方法擦除了原来的单位圆，并画出更为精确的单位圆。ucf的命令代表单位圆修正器，可以在没有影响的前提下对任何Matlab工具箱的函数进行改进。

●不推荐的方法：编辑Matlab的m文件zplaneplot的第91行，如下：

linspace命令创建了一个变量theta，它由均匀分布于0～2π的70个元素组成。进行变换，将第70元素作为更大的整数如100，将导致使用更多点来绘制单位圆。由于它修改了美国MathWorks公司提供的Matlab的工具箱函数所以不推荐。修改花钱购买的代码会给其他人造成开发和维护的困难，因如下四个原因而不推荐。

（1）如果你认为这个代码有错误，可以提出修正代码的申请，这样也能使得其他人从中获益。

（2）如果你修改了Matlab的代码，下一次升级这个工具箱的时候，安装的过程中会覆盖这个你已经修改的函数，这样你的修改就失效了。

（3）很久以后你还记得所做的修改的可能性是很小的，你可能在下次安装工具箱时再次发生单位圆绘制的问题。

（4）如果你破坏了工具箱函数，它不能正常工作或者表现更坏，但你仍然认为它能很好地工作，这实际上将返回不正确的结果。

你在别人的软件中会不可避免地遇到这些问题，这时要记得上面四点，只有从长远考虑才能更好地避免遇到问题。

创建群延时图

使用Matlab的grpdelay命令可以创建一个群延时的图形。

这样会产生如图4.11所示的曲线。主要群延时不是一个常数，这是由这个滤波器的非线性相位响应所决定的。

图4.11 与一个4阶巴特沃斯低通滤波器相关的群延时，它的截止频率为0.25F_S/2

使用FDATool和FVtool

Matlab软件环境提供了更多的设计数字滤波器的工具。例如，你可以使用Matlab的FDATool（Filter Design and Analysis Tool，滤波器设计和分析工具箱）来设计你的滤波器。如图4.12所示，软件的一些按钮可以用于设计你的滤波器，也可以分析并对滤波器性能作图。FDATool输出的滤波器的系数是矢量（默认情况下）而不是像SPTool那样嵌入到filt结构中。

图4.12 FDATool的滤波器设计和分析程序的图形用户程序

最后，不管如何设计你的滤波器，你可能通过使用Matlab的FVTool（Filter Visualization Tool，滤波器观测工具）来分析你的滤波器，格式如下：

如果你的滤波器参数是如下的结构形式的，这个Matlab命令为

图4.13所示为一些FVTool中软件的按钮，它允许你观察滤波器分析的图形，通过图中标签和箭头可以明确它们的含义。在这个图中选择了滤波器的群延时。

虽然所有的Matlab的命令和工具是就IIR滤波器而说的，但是可以非常容易地用于FIR滤波器设计。其中，主要的不同是FIR滤波器的“A”矢量的系数等于数值1。