例说信号处理与滤波器设计

例说信号处理与滤波器设计

许多公式在转换时成了乱码，相应的word版本请点这里

目录

数字时代    2

数字信号处理的应用    3

频率——信号的指纹    5

卷积可以不卷    8

向量运算的启示    11

滤波器设计征程    16

最后一击——滤波的实现方法    22

纵览全局    27

数字时代

信号处理是对原始信号进行改变，以提取有用信息的过程，它是对信号进行变换、滤波、分析、综合等处理过程的统称。数字信号处理是将信号以数字方式表示并处理的理论和技术；模拟信号处理是指用模拟系统对模拟信号进行处理的方法或过程。

    数字信号处理课程的主要内容包括信号分析与处理。两者并不是孤立的，不同的信号处理方法往往需要选择不同的信号表示形式。两者的区别主要表现在，信号处理是用系统改变输入信号，以得到所期望的输出信号，如信号去噪；而信号分析往往是通过变换(傅里叶变换、小波变换等)，或其它手段提取信号的某些特征，如语音信号的基本频率，图像的直方图等。

早期的信号处理局限于模拟信号，随着数字计算机的飞速发展，信号处理的理论和方法得以飞速发展，出现了不受物理制约的纯数学的加工，即算法，并确立了数字信号处理的领域。现在，对于信号的处理，人们通常是先把模拟信号变成数字信号，然后利用高效的数字信号处理器(DSP：Digital Signal Processor)或计算机对其进行数字形式的信号处理。

　　一般地讲，数字信号处理涉及三个步骤：

1. 模数转换(A/D转换)：把模拟信号变成数字信号，是一个对自变量和幅值同时进行离散化的过程，基本的理论保证是采样定理。
2. 数字信号处理(DSP)：包括变换域分析(如频域变换)、数字滤波、识别、合成等。
3. 数模转换(D/A转换)：把经过处理的数字信号还原为模拟信号。通常，这一步并不是必须的。

图 1数字信号处理基本步骤

数字信号处理的应用

图 2家庭影院(视+听)

高分辨率图像、高保真音质

图 3语音识别

噪声环境下高识别率

图 4图像增强

更清晰、美观

图 5无人驾驶

高度智能、安全

图 6医学成像

更高效、更精确的成像结果

狭义地讲，信号处理可以统称为滤波，根据不同的要求，选用不同性能的滤波器。在数字信号处理应用中，设计合适的滤波器至关重要。什么是数字滤波器呢？数字滤波器就是数字信号处理(Digital Signal Processing)算法，或者是数字信号处理器件(Digital Signal Processor)。什么是数字信号处理算法呢？我们需要借助基本的数学工具和方法。

首先回到信号处理的根本目的：用滤波器改变输入信号，以期得到理想的输出信号，如果输入信号用$x\left[ n \right]$表示，输出信号用$y\left[ n \right]$表示，则数字信号处理或滤波可用如下框图表示：

图 7数字信号滤波图示

上述框图并没有给出数字信号处理算法的具体实现方法，但提供了数字信号处理的基本数学思想：通过系统将$x\left[ n \right]$和$y\left[ n \right]$进行关联。在数学上描述$x$和$y$之间关系的手段有很多，如函数$y = f\left( x \right)$，方程组\[Ax = y\]等，从数学角度来理解，如果$f\left( \cdot \right)$为不同映射，即使同样的自变量$x$，也会得到不同的因变量$y$；如果系数矩阵\[A\]不同，同样的\[x\]经变换之后，将得到不同的\[y\]。在数字信号处理课程中，描述输入输出信号关系的基本数学工具是差分方程：

    \[\sum\limits_{k = 0}^Q {{a_k}y\left[ {n - k} \right]} = \sum\limits_{k = 0}^P {{b_k}x\left[ {n - k} \right]} \]    

令$N = \max \left( {P,Q} \right)$，上述差分方程可写为

    \[\sum\limits_{k = 0}^N {{a_k}y\left[ {n - k} \right]} = \sum\limits_{k = 0}^N {{b_k}x\left[ {n - k} \right]} \]    

如果仅局限于基础(经典)的数字信号处理算法研究，全部内容都是围绕这个常系数线性差分方程(Constant Coefficient Linear Difference Equation，CCLDE)展开的。那么，我们需要研究CCLDE什么呢？很显然，我们可以将CCLDE看成一个系统，它描述了输入$x\left[ n \right]$和输出$y\left[ n \right]$之间的关系。哪些参数决定这个系统的性能呢？方程中除了$x\left[ n \right]$和$y\left[ n \right]$，还有三个参数：方程阶数$N$，系数${a_k}$和${b_k}$——数字滤波器设计的根本任务就是确定这三个参数——目标非常明确!!!

    从何下手呢？如何判定所设计的滤波器符合预期呢？迷茫中ing…

频率——信号的指纹

    烈日当空，窗外的知了热得撕心裂肺地吼叫

声音特点：又大又尖，怎叫人不心烦意乱？!

大：能量大，或者信号幅度大(声压级超过120dB即达到痛阈)；

尖：频率高，即知了的叫声中包含了丰富的高频成分，因此听起来很"刺耳"。

图 8知了声的频谱图(横坐标：频率(Hz)；纵坐标：幅度(dB))

图 9 听觉曲线

图 10 声音频率范围

可见，幅度和频率是声音信号的主要参数。诸如"低音炮"、"高音喇叭"、"男低音"和"女高音"都是从幅度和频率去描述声音信号的，因此，分析信号的频率特性是信号处理领域的重要内容。

分析信号的频率特性无非就是想知道信号包含了哪些频率成分，各频率成分的大小是多少。傅里叶分析完美地解决了这个问题。在信号与系统课程中，我们通常将傅里叶分析分为傅里叶级数和傅里叶变换，前者用来对付周期信号，后者用来处理非周期信号，但无论是傅里叶级数，还是傅里叶变换(当然，傅里叶级数也可以纳入到傅里叶变换体系中)，它们的原理或宗旨都一样：将一般的信号分解为基本信号的线性组合

(连续周期信号)    \[\tilde x\left( t \right) = \sum\limits_{k = - \infty }^{ + \infty } {{a_k}{e^{jk{\omega _0}t}}} \]    

(离散非周期序列)    \[x\left[ n \right] = \frac{1}{{2\pi }}\int_{\left\langle {2\pi } \right\rangle } {X\left( {{e^{j\omega }}} \right)} {e^{j\omega n}}d\omega \]    

以式为例，等式左边是一般的(满足收敛条件的)周期信号，而右边则是频率为$k{\omega _0}$、幅度为${a_k}$的虚指数信号\[{e^{jk{\omega _0}t}}\]的线性组合。

现在可以谈论一下滤波的概念了。假设式中\[\tilde x\left( t \right)\]就是知了发出的"嗞哇嗞哇……"没完没了的周期信号，现在我们想把烦人的高频成分去掉，最直观的想法就是需要设计一个滤波器，这个滤波器可以让较低频率的信号顺利通过，同时又能阻止较高频率的信号，这就是所谓的低通滤波器。

卷积可以不卷

    再回到常系数线性差分方程，参数$N$，${a_k}$和${b_k}$完全决定了方程所描述的系统的所有特性。什么？难道与输入$x\left[ n \right]$和输出$y\left[ n \right]$没关系？对，没有关系!电阻是一个系统，其阻值仅与自身的材料及结构有关，虽然有关系式$R = v\left( t \right)/i\left( t \right)$，事实上，即使两端没有电压，电阻依然存在；再比如理想的线性放大电器，它的增益仅取决于内部结构，而与输入输出无关，但为了测量放大器的增益(放大倍数)，我们可以在输入端接入幅度为$i$的信号，然后测量输出信号的幅度$o$，这样就可以得到放大倍数$g = \frac{o}{i}$，更特殊地，如果取$i = 1$，此时输出信号的幅度就是放大器的增益，即$g = o$。类似的概念推广到CCLDE描述的LTI系统，我们如何获得这个系统的特性呢？输入——输出描述法，即令输入$x\left[ n \right]$取某种特殊值时，计算(或测量)系统的输出$y\left[ n \right]$。那么，什么样的$x\left[ n \right]$算特殊呢？单位脉冲!

    \[\delta \left[ n \right] = \left\{ \begin{array}{l}

1,\;\;\;n = 0\\

0,\;\;n \ne 0

\end{array} \right.\]    

也就是说，令$x\left[ n \right] = \delta \left[ n \right]$，此时系统的输出就是所谓的单位脉冲响应$h\left[ n \right]$。再与放大器的例子对比一下：输入信号幅度$i = 1$，输出信号的幅度$o$就是放大倍数(放大器的重要指标，从应用者角度而言，其实就是唯一关心的指标)。我们有理由相信，对于CCLDE描述的LTI系统，得到了单位脉冲响应$h\left[ n \right]$，就能够掌握系统的全部特性(因果性、稳定性、频率选择性等)。何以见得呢？之后将逐步解答。

    现在，我们假设已经知道了LTI系统的单位脉冲响应$h\left[ n \right]$，对于任意输入$x\left[ n \right]$，如何求系统的输出$y\left[ n \right]$呢？其实就是求滤波后的信号。首先，我们要建立已知和需求之间的联系。已知的是

，需求是

。

图 11 单位脉冲响应

事实上，任意序列$x\left[ n \right]$很容易用$\delta \left[ n \right]$的移位、加权的线性组合来表示：$x\left[ n \right] = \sum\limits_{k = - \infty }^{ + \infty } {x\left[ k \right]\delta \left[ {n - k} \right]} $。

图 12 序列的单位脉冲表示

结合系统的线性和时不变性，有

式中最后一个等式\[y\left[ n \right] = \sum\limits_{k = - \infty }^{ + \infty } {x\left[ k \right]h\left[ {n - k} \right]} \]就是卷积和，记作$y\left[ n \right] = x\left[ n \right]*h\left[ n \right]$，表明LTI系统的(零状态)响应$y\left[ n \right]$是单位脉冲响应$h\left[ n \right]$的移位、加权、求和。这句话也描述了求$y\left[ n \right]$的步骤。根据图 4，我们可以假想有这样的一个LTI系统——出钞机：如果你今天往出钞机投币口(输入端)投币一元，它将在接下来的三天，每一天都会从出币口(输出端)吐出一元。试问：如果某个月的1、3、4日，你分别投了2、2、3元，那么，出钞机每天的输出是多少元呢？不妨用图来说明一下。

图 13 卷积和的图解

图解结果告诉我们，出钞机在2-7日分别吐出2,2,4,5,5,3元。图解的过程包含了乘积与求和("积"与"和")，但并有体现出"卷"(翻转)这一操作。如果只想求某一天(比如5日)出钞口吐出多少钱，此时就要用另一种方法，即许多教材中描述的步骤：

1. 变量替换：$x\left[ n \right] \to x\left[ k \right],h\left[ n \right] \to h\left[ k \right]$
2. 将$x\left[ k \right]$或$y\left[ k \right]$翻转
3. 移位、相乘、相加

以上步骤包含了"卷积和"所有的操作(卷——翻转；积——相乘；和——相加)。

图解过程是根据线性时不变系统的定义导出的一种结果，是系统特性的直接反映。卷积和是信号与系统、数字信号处理中最重要的公式之一，它描述了LTI系统输入$x\left[ n \right]$、输出$y\left[ n \right]$以及单位脉冲响应$h\left[ n \right]$之间的关系，已知三者中的任何两个，就可以确定第三者，于是就有以下应用场合：

1. 已知输入$x\left[ n \right]$和系统单位脉冲响应$h\left[ n \right]$，求输出$y\left[ n \right]$，即用确定的系统对输入信号滤波处理；
2. 已知输出$y\left[ n \right]$和系统单位脉冲响应$h\left[ n \right]$，求输入$x\left[ n \right]$；

3. 已知输入$x\left[ n \right]$和输出$y\left[ n \right]$，求系统单位脉冲响应$h\left[ n \right]$。

   现在，我们将注意力集中到卷积和公式上来，公式重写如下

    \[y\left[ n \right] = \sum\limits_{k = - \infty }^{ + \infty } {x\left[ k \right]h\left[ {n - k} \right]} \]    

卷积和只包含了移位、数乘和相加运算，看来，数字信号处理的计算十分简单。如果再上升一个层次，卷积和能很直观解释系统如何改变输入信号的吗？从图 6的结果很难看出系统是如何将$x\left[ n \right]$变成$y\left[ n \right]$的，而且也很难说明这种变化意味着什么，因此，卷积和在解释系统行为特性时很不直观，这是卷积和运算本身决定的。

向量运算的启示

了解运算背后的物理或几何含义对理解问题的本质至关重要。例如向量的乘法和加法：设两个向量\[\vec A = {x_1} + i{y_1}\]，\[\vec B = {x_2} + i{y_2}\]，向量的加法符合平行四边形法则，从图 7很容易看出两个向量相加的几何意义，而相应的代数运算也很简单——对应分量相加即可，即有

    \[\vec C = \vec A + \vec B = \left( {{x_1} + {x_2}} \right) + i\left( {{y_1} + {y_2}} \right)\]    

但如果要求两个向量的乘法呢？代数运算的结果是

相比两个向量的加法，向量的乘法运算要复杂一些，更重要的是，根据运算结果，我们很难看清向量乘法的几何意义。如果将向量$\vec A$和$\vec B$用极坐标表示，情况又会如何呢？设

    \[\begin{array}{l}

\vec A = \left( {{\rho _1},{\theta _1}} \right) = {\rho _1}{e^{j{\theta _1}}}\\

\vec B = \left( {{\rho _2},{\theta _2}} \right) = {\rho _2}{e^{j{\theta _2}}}

\end{array}\]    

则两个向量的乘法

计算简单，几何意义明确——模相乘，相位相加。

图 14 向量加的直角坐标表示

图 15 向量乘法的极坐标表示

    向量运算的例子告诉我们，选择不同的坐标系，不仅影响运算的复杂度，在解释运算的几何意义时也各有千秋。上例告诉我们，如果是向量的加法，直角坐标系比极坐标系方便；如果是向量的乘法，很显然，极坐标系比直角坐标系方便。

    回到数字信号处理的话题，卷积和是让无数人困惑的公式，然而，它又是经典数字信号处理算法的根源。既然卷积和运算就是滤波，我们如何评判滤波的效果呢？从什么角度理解系统的滤波行为更好呢？

    再想想知了的叫声，之所以烦人，是因为包含了大量的高频分量。看来，用频率这一参数来解释信号特性很符合人的直观感觉。我们有更好的方式来洞悉卷积和的物理意义吗？能不能站在频(率)域的角度来解释卷积和呢？当然可以，因为我们有卷积定理(性质)：

通过傅里叶变换，将时域中的卷积和运算换成了频域中的相乘运算。那么频域中的相乘运算有什么好处呢？不要忘了，我们主要是为了直观解释滤波的行为和特性。

    图 9是一段滤波前知了声$x\left[ n \right]$的时域波形，单从图形上看，似乎很难说出这一段信号的特点，运用傅里叶变换，得到$x\left[ n \right]$的幅度谱\[\left| {X\left( {{e^{j\omega }}} \right)} \right|\]如图 10所示，很明显，该图说明信号$x\left[ n \right]$包含了丰富的高频成分。现在，我们希望设计一个滤波器，可以滤除2500Hz以上的频率分量，目的是仅保留红框内的频率分量，因此要求滤波器的频率响应(幅频响应)是

滤波器的幅频特性如图 11所示，截止频率${f_c} = 2500$Hz。由卷积定理可知，输出信号的幅度谱为

    \[\left| {Y\left( {{e^{j\omega }}} \right)} \right| = \left| {X\left( {{e^{j\omega }}} \right)} \right|\left| {H\left( {{e^{j\omega }}} \right)} \right|\]    

从图形上来看就更直观了：图 10与图 11相乘，得到了滤波后信号的幅度谱，滤波前、后信号频谱的变化一目了然，很显然，它是按照滤波器频率响应\[H\left( {{e^{j\omega }}} \right)\]所规定的形状去改变的。而对比滤波前、后的时域波形，很难解释信号的什么特性被改变了。

    在上述的讨论中，我们只关注幅度谱，实际中，相位谱也必须考虑，更完整的表达式为

对比极坐标系下两个向量的乘法，能找到相似之处吗？——模相乘、相位相加。如果将卷积和$y\left[ n \right] = x\left[ n \right]*h\left[ n \right]$看作直角坐标系下两个向量的乘法，那么\[Y\left( {{e^{j\omega }}} \right) = X\left( {{e^{j\omega }}} \right)H\left( {{e^{j\omega }}} \right)\]就类似于极坐标系下两个向量的乘法——运算简单、物理意义直观。[通过向量乘法运算的例子，我们可以看到，向量的表示方式起到了关键性的作用；同样地，信号也有不同的描述方式，如时域、频域、复频域，不同的表示方法都是为了使分析问题更简单、物理意义更明确]

图 16滤波前知了声$x\left[ n \right]$时域波形

图 17滤波前知了声的幅度谱\[\left| {X\left( {{e^{j\omega }}} \right)} \right|\] (红框是希望保留的低频分量)

图 18低通滤波器幅频特性\[\left| {H\left( {{e^{j\omega }}} \right)} \right|\] (蓝色线)

图 19滤波后知了声的幅度谱\[\left| {Y\left( {{e^{j\omega }}} \right)} \right|\]

图 20滤波后知了声$y\left[ n \right]$时域波形

滤波器设计征程

    借助系统的频率响应\[H\left( {{e^{j\omega }}} \right)\]，可以很方便地解释滤波器的特性(主要是指频率选择性)。有什么定量指标来描述频率选择性吗？通常有四个指标——横坐标两个，纵坐标两个。以低通滤波器为例来说明。如图 14所示，横坐标两个指标为通带截止频率\[{\omega _p}\]和阻带截止频率\[{\omega _s}\]，纵坐标两个指标为通带衰减\[{A_p}\]和阻带衰减\[{A_s}\]。

图 21数字低通滤波器指标

    通常而言，滤波器设计就是指根据应用的需求，提出合理的滤波器指标(通带截止频率和衰减，阻带截止频率和衰减)，然后借助滤波器设计工具得到滤波器参数(系数)。滤波器系数？没听说过啊，是什么东西？其实我们早见过面了，就是差分方程中的${a_k}$和${b_k}$!早就说过，数字滤波器设计要不忘初心，根本任务就是确定这三个参数——方程的阶数$N$，系数${a_k}$和${b_k}$。

    如何建立起滤波器指标和滤波器系数之间的联系呢？很显然，滤波器的指标约束着其频率响应\[H\left( {{e^{j\omega }}} \right)\]，而\[H\left( {{e^{j\omega }}} \right)\]与CCLDE的阶数$N$，系数${a_k}$和${b_k}$直接相关。傅里叶变换又该出场了。对差分方程两边取傅里叶变换得

    \[\sum\limits_{k = 0}^N {{a_k}{e^{ - jk\omega }}Y\left( {{e^{j\omega }}} \right)} = \sum\limits_{k = 0}^N {{b_k}{e^{ - jk\omega }}X\left( {{e^{j\omega }}} \right)} \]    

从而可得频率响应

    \[H\left( {{e^{j\omega }}} \right) = \frac{{Y\left( {{e^{j\omega }}} \right)}}{{X\left( {{e^{j\omega }}} \right)}} = \frac{{\sum\limits_{k = 0}^N {{b_k}{e^{ - jk\omega }}} }}{{\sum\limits_{k = 0}^N {{a_k}{e^{ - jk\omega }}} }}\]    

易见，频率响应\[H\left( {{e^{j\omega }}} \right)\]完全由$N$，${a_k}$和${b_k}$确定。综上可知，滤波器设计就是如何根据滤波器指标求解滤波器系数的过程，从数学角度来看，就是函数逼近和数值优化的过程。经典滤波器设计理论发展较完备(巴特沃斯、切比雪夫I、切比雪夫II、椭圆……，各自的特点是什么？更多内容，请点这里)，许多教科书都有详细讨论，此处不再赘述，我们的重点是要借助现代化的设计工具，轻松完成设计任务并实施滤波处理。此处仅以MATLAB工具为例，并通过两种方式来实现滤波器设计——脚本代码和可视化方法。

    需求：设计一个低通滤波器来处理知了声，通带截止频率\[{F_p} = 2000\]Hz，通带衰减${A_p} = 0.1$dB，阻带截止频率\[{F_s} = 2500\]Hz，通带衰减${A_p} = 50$dB。

乍一看，这些指标与图 14中的指标完全不同，实际需求的指标往往同本例，频率指标为模拟频率，衰减指标以dB为单位，而设计数字滤波器时，往往需要数字频率指标。如何进行指标转换呢？直接用公式计算${\omega _p} = 2\pi {F_p} = 4000\pi $，${\omega _s} = 2\pi {F_s} = 5000\pi $吗？图 14告诉我们，无论通带截止频率，还是阻带截止频率都不会超过$\pi $。难道最大的数字角频率就是$\pi $吗？为什么？这个问题同样令许多人困惑。如何衡量数字频率的大小呢？先考虑一个简单的数字信号：$x\left[ n \right] = \cos \left( {\omega n} \right)$，针对$\omega $取不同的值，分别画出对应的波形如下：

图 22$\cos \left( {\omega n} \right)$波形图

振荡得越厉害，意味着频率越高，这是符合常理的。上图中的9个波形，哪个振荡得最厉害呢？第二排中间那个——+1和-1交替，而此时$\omega = \pi $。而且还发现，$\omega = 0$和$\omega = 2\pi $的波形一样，并不是$\omega $的值越大，振荡得越快。如果是连续时间信号$x\left( t \right) = \cos \left( {2\pi ft} \right) = \cos \left( {\Omega t} \right)$，大家都知道，随着$\Omega $的增大，波形会振荡得越快。从连续时间信号到离散时间信号，是什么原因导致两者出现如此大的差异呢？通过对采样过程的探索，我们便可理解其中的奥妙。设以时间间隔${T_s}$对$x\left( t \right) = \cos \left( {\Omega t} \right)$采样，这样便得到

    $x\left( {n{T_s}} \right) = \cos \left( {\Omega n{T_s}} \right) = \cos \left( {\omega n} \right)$    

因此有

    $\omega = \Omega {T_s} = \Omega /{f_s}$    

其中\[{f_s} = 1/{T_s}\]为采样频率，$\Omega $为模拟角频率，$\omega $为数字角频率。结论是：数字频率是模拟频率关于采样频率的归一化。

    根据Nyquist采样定理，当采样频率大于信号最高频率2倍时，可以由样点完全恢复原来的信号。我们不妨假设信号的最高频率正好是采样频率的一半，即\[{f_{\max }} = {f_s}/2\]，对应的最高模拟角频率${\Omega _{\max }} = 2\pi {f_{\max }} = \pi {f_s}$，代入式得

    ${\omega _{\max }} = {\Omega _{\max }}/{f_s} = \pi {f_s}/{f_s} = \pi $    

这也证明了为什么最高数字角频率为$\pi $。(如果不满足Nyquist采样条件又会怎样呢？事实上并不影响结论成立)。关于各种频率之间的关系的讨论，可以参考这里(注意文中"fdtool工具中归一化频率为什么最大只到0.5的原因"，应改为"最大只到1的原因"。)

    再回到滤波器设计题目，第一步，需要将模拟频率指标转换为归一化数字频率指标，转换公式是

    \[\omega = 2f/{f_s}\;\;\;\;(\pi )\]    

很显然，我们需要知道采样频率${f_s}$。知了声的采样频率${f_s} = 11025$Hz。将模拟频率指标代入式计算出对应的归一化数字频率指标

    \[\begin{array}{l}

{\omega _p} = 2{F_p}/{f_s} = 2 \times 2000/11025 \approx 0.36\\

{\omega _s} = 2{F_s}/{f_s} = 2 \times 2500/11025 \approx 0.45

\end{array}\]    

这样，我们就得到了数字域的四个指标：${\omega _p},{A_p},{\omega _s},{A_s}$。如何从这四个指标去计算滤波器系数呢？有现成的公式套用，即许多教科书上所说的各类滤波器原型。我们选择椭圆型滤波器为例。MATLAB代码如下：

clc;

close all

fs=11025;                 %采样频率(Hz)

Fp = 2000;             %通带截止频率(Hz)

Fs = 2500;             %阻带截止频率(Hz)

Ap = 0.1;             %通带衰减(dB)

As = 50;             %阻带衰减(dB)

wp=2*Fp/fs;             %归一化通带截止频率

ws=2*Fs/fs;             %归一化阻带截止频率 

[N,Wp]=ellipord(wp,ws,Ap,As);         %确定带通滤波器的阶数和截止频率 

[b,a]=ellip(N,Ap,As,Wp);     %确定滤波器系数 

[h,w]=freqz(b,a);     %求数字带通滤波器的频率响应 

%以下为绘图命令，绘制带通滤波器的幅频响应 

figure;

plot(w*fs/(2*pi),20*log10(abs(h)/max(abs(h))));

axis([0,fs/2,-100,0]);

title(\'数字低通滤波器的幅度响应\');

xlabel(\'频率(Hz)\');

ylabel(\'幅度(dB)\');

grid

滤波器系数

a = [1 -2.98876196757509 5.41154092244290 -6.18465137960809    4.94650037623018     -2.65349822687766 0.903727225595414 -0.150115115698030]

b = [0.0204237714181223    0.0185281305134523    0.0518202070683316    0.0515988082549051    0.0515988082549052    0.0518202070683315    0.0185281305134523    0.0204237714181222]

分别对应式中的${a_k}$和${b_k}$，得到了这两个参数，滤波器设计基本上算完成了。

    如果你掌握了滤波器设计的基本理论，但又不想写代码，你可以利用FDATool轻松完成滤波器设计的全过程。

在命令窗口输入：fdatool，回车，就可以看到下图所示的界面：

然后动动鼠标就可以完成滤波器设计了。针对本例，我们按下图红框所示设定指标，在完成指标设定之后，点"Design Filter"按钮即完成设计。关于滤波器的所有信息，都可以通过菜单或工具栏获得。注意"Current Filter Information"部分，它描述了当前所设计出的滤波器基本信息，包括结构、阶次、稳定性等。有一个问题值得思考：滤波器的结构对滤波器的工程实现有什么影响呢？成本、稳定性、抗噪声性能，这些都与结构有关，因此，掌握滤波器各种结构的优缺点也十分重要。

    很显然，用FDATool设计滤波器十分简单!无论采用什么方法，都不要忘记滤波器设计的根本目标——得到滤波器系数。关于FDATool更多内容，请点这里和这里。

    当滤波器设计完成之后，我们该考虑如何用所得到的系数对输入信号进行滤波了。

最后一击——滤波的实现方法

    当滤波器系数${a_k}$和${b_k}$确定之后，方程也就完全确定了，所谓的滤波，就是要对给定的输入$x\left[ n \right]$，求系统的响应$y\left[ n \right]$。MATLAB提供了几种滤波实现函数：

conv：线性卷积，可以用来实现滤波

filter：滤波，但是输出点数和输入点数相同，其余的存在内部状态中，以便处理后续输入

filtfilt：零相位滤波，输出通过滤波器后再反向通过滤波器，消除相位影响

fftfilt：利用基于FFT的重叠相加法对数据进行滤波，这种频域滤波技术只对FIR滤波器有效

完整例子：用MATLAB读入知了声的音频文件，画时域波形和频谱，设计滤波器，对音频文件进行滤波，最后画滤波后的信号时域波形和频谱。程序如下：

clc;

close all

[xn, fs] = wavread(\'cicada.wav\'); %读wav文件(知了声)，获得样点值xn和采样频率fs

%计算信号xn的频谱

xn_FFT = fft(xn);     %FFT频谱分析

xn_FFT_AMP = abs(xn_FFT);     %求xn的幅度谱

xn_FFT_AMP_db = 20*log10(xn_FFT_AMP/max(xn_FFT_AMP)); %用对数表示

%画xn时域波形

figure

n = 1:length(xn);         %生成样点横坐标向量

plot(n,xn);     %以n为横坐标，画xn

axis([1 length(xn) -1 1])     %设定坐标轴显示范围

xlabel(\'样点\');     %设定x坐标标记

ylabel(\'幅度\');     %设定y坐标标记

title(\'知了声时域波形\');     %设定标题

%画xn的幅度谱(对数形式)

figure

faxis = (0:length(xn)-1)*fs/length(xn);         %生成横坐标(频率轴)刻度，FFT的频率范围[0 fs]

plot(faxis,xn_FFT_AMP_db);             %以faxis为横轴画xn频谱图

axis([0 fs/2 -100 0]) %设定显示范围，由于FFT频谱具有对称性，因此只显示前半部分

xlabel(\'频率(Hz)\');

ylabel(\'幅度(dB)\');

title(\'知了声幅度谱\');

% 设计低通滤波器

Fp = 2000;     %通带截止频率(Hz)

Fs = 2500;     %阻带截止频率(Hz)

Ap = 0.1;     %通带衰减(dB)

As = 50;     %阻带衰减(dB)

wp=2*Fp/fs;     %归一化通带截止频率

ws=2*Fs/fs;     %归一化阻带截止频率

[N,Wp]=ellipord(wp,ws,Ap,As);             %确定带通滤波器的阶数和边缘频率

[b,a]=ellip(N,Ap,As,Wp);     %确定滤波器系数

[h,w]=freqz(b,a); %求数字带通滤波器的频率响应

%以下为绘图命令，绘制带通滤波器的幅频响应

figure;

plot(w*fs/(2*pi),20*log10(abs(h)/max(abs(h))));

axis([0,fs/2,-100,0]);

title(\'数字低通滤波器的幅度响应\');

xlabel(\'频率(Hz)\');

ylabel(\'幅度(dB)\');

grid

%利用设计的滤波器对xn进行滤波，得到滤波输出yn

yn = filter(b,a,xn);                     %滤波函数，b和a为滤波器系数，xn是输入信号

wavwrite(yn,fs,8,\'cicada_LPFiltered.wav\'); %保存滤波后的信号

%画滤波后的信号yn时域波形

figure

n = 1:length(yn);

plot(n,yn);

axis([1 length(yn) -1 1])

xlabel(\'样点\');

ylabel(\'幅度\');

title(\'低通滤波后知了声时域波形\');

%计算滤波后信号yn的频谱

yn_FFT = fft(yn);

yn_FFT_AMP = abs(yn_FFT);

yn_FFT_AMP_db = 20*log10(yn_FFT_AMP/max(yn_FFT_AMP));

figure

faxis = (0:length(yn)-1)*fs/length(yn);

plot(faxis,yn_FFT_AMP_db);

axis([0 fs/2 -100 0])

xlabel(\'频率(Hz)\');

ylabel(\'幅度(dB)\');

title(\'滤波后知了声幅度谱\');

结果：

图 23知了声时域波形

图 24知了声频谱(幅度频)

图 25低通滤波器频率响应(幅频响应)

图 26低通滤波后知了声时域波形

图 27滤波后知了声频谱图(幅度谱)

纵览全局

经典信号处理课程的主要内容：

1. 信号的表示方法(时域表示、频域表示、复频域表示)
2. 系统的描述方法(CCLDE、单位脉冲响应、频率响应、系统函数、零极点图)
3. 滤波器设计方法(代码和可视化方法)
4. 滤波的实施方法(滤波器类型的选择——IIR、FIR；滤波器结构的选择，量化对系数的影响)

欢迎加入数字信号处理群：152346662