信号分析与处理概述
博客对图片和公式的支持不是很好，相应的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)。什么是数字信号处理算法呢？我们需要借助基本的数学工具和方法。
首先回到信号处理的根本目的：用滤波器改变输入信号，以期得到理想的输出信号，如果输入信号用$xleft[ n ight]$表示，输出信号用$yleft[ n ight]$表示，则数字信号处理或滤波可用如下框图表示：
图 7数字信号滤波图示
上述框图并没有给出数字信号处理算法的具体实现方法，但提供了数字信号处理的基本数学思想：通过系统将$xleft[ n ight]$和$yleft[ n ight]$进行关联。在数学上描述$x$和$y$之间关系的手段有很多，如函数$y = fleft( x ight)$，方程组$Ax = y$等，从数学角度来理解，如果$fleft( cdot ight)$为不同映射，即使同样的自变量$x$，也会得到不同的因变量$y$；如果系数矩阵$A$不同，同样的$x$经变换之后，将得到不同的$y$。在数字信号处理课程中，描述输入输出信号关系的基本数学工具是差分方程：
[sumlimits_{k = 0}^Q {{a_k}yleft[ {n - k} ight]} = sumlimits_{k = 0}^P {{b_k}xleft[ {n - k} ight]} ]
令$N = max left( {P,Q} ight)$，上述差分方程可写为
[sumlimits_{k = 0}^N {{a_k}yleft[ {n - k} ight]} = sumlimits_{k = 0}^N {{b_k}xleft[ {n - k} ight]} ]
如果仅局限于基础(经典)的数字信号处理算法研究，全部内容都是围绕这个线性常系数差分方程(Linear Constant-Coefficient Difference Equation，LCCDE)展开的。那么，我们需要研究LCCDE什么呢？很显然，我们可以将LCCDE看成一个系统，它描述了输入$xleft[ n ight]$和输出$yleft[ n ight]$之间的关系。哪些参数决定这个系统的性能呢？方程中除了$xleft[ n ight]$和$yleft[ n ight]$，还有三个参数：方程阶数$N$，系数${a_k}$和${b_k}$——数字滤波器设计的根本任务就是确定这三个参数——目标非常明确!!!
从何下手呢？如何判定所设计的滤波器符合预期呢？迷茫中ing…
频率——信号的指纹
烈日当空，窗外的知了热得撕心裂肺地吼叫
声音特点：又大又尖，怎叫人不心烦意乱？!
大：能量大，或者信号幅度大(声压级超过120dB即达到痛阈)；
尖：频率高，即知了的叫声中包含了丰富的高频成分，因此听起来很“刺耳”。
图 8知了声的频谱图(横坐标：频率(Hz)；纵坐标：幅度(dB))
图 9 听觉曲线
图 10 声音频率范围
可见，幅度和频率是声音信号的主要参数。诸如“低音炮”、“高音喇叭”、“男低音”和“女高音”都是从幅度和频率去描述声音信号的，因此，分析信号的频率特性是信号处理领域的重要内容。
分析信号的频率特性无非就是想知道信号包含了哪些频率成分，各频率成分的大小是多少。傅里叶分析完美地解决了这个问题。在信号与系统课程中，我们通常将傅里叶分析分为傅里叶级数和傅里叶变换，前者用来对付周期信号，后者用来处理非周期信号，但无论是傅里叶级数，还是傅里叶变换(当然，傅里叶级数也可以纳入到傅里叶变换体系中)，它们的原理或宗旨都一样：将一般的信号分解为基本信号的线性组合
(连续周期信号) [ ilde xleft( t ight) = sumlimits_{k = - infty }^{ + infty } {{a_k}{e^{jk{omega _0}t}}} ]
(离散非周期序列) [xleft[ n ight] = frac{1}{{2pi }}int_{leftlangle {2pi } ight angle } {Xleft( {{e^{jomega }}} ight)} {e^{jomega n}}domega ]
以式为例，等式左边是一般的(满足收敛条件的)周期信号，而右边则是频率为$k{omega _0}$、幅度为${a_k}$的虚指数信号${e^{jk{omega _0}t}}$的线性组合。
现在可以谈论一下滤波的概念了。假设式中$ ilde xleft( t ight)$就是知了发出的“嗞哇嗞哇……”没完没了的周期信号，现在我们想把烦人的高频成分去掉，最直观的想法就是需要设计一个滤波器，这个滤波器可以让较低频率的信号顺利通过，同时又能阻止较高频率的信号，这就是所谓的低通滤波器。
卷积可以不卷
再回到线性常系数性差分方程，参数$N$，${a_k}$和${b_k}$完全决定了方程所描述的系统的所有特性。什么？难道与输入$xleft[ n ight]$和输出$yleft[ n ight]$没关系？对，没有关系!电阻是一个系统，其阻值仅与自身的材料及结构有关，虽然有关系式$R = vleft( t ight)/ileft( t ight)$，事实上，即使两端没有电压，电阻依然存在；再比如理想的线性放大电器，它的增益仅取决于内部结构，而与输入输出无关，但为了测量放大器的增益(放大倍数)，我们可以在输入端接入幅度为$i$的信号，然后测量输出信号的幅度$o$，这样就可以得到放大倍数$g = frac{o}{i}$，更特殊地，如果取$i = 1$，此时输出信号的幅度就是放大器的增益，即$g = o$。类似的概念推广到LCCDE描述的LTI系统，我们如何获得这个系统的特性呢？输入——输出描述法，即令输入$xleft[ n ight]$取某种特殊值时，计算(或测量)系统的输出$yleft[ n ight]$。那么，什么样的$xleft[ n ight]$算特殊呢？单位脉冲!
[delta left[ n
ight] = left{ egin{array}{l}
1,;;;n = 0\
0,;;n
e 0
end{array}
ight.]
也就是说，令$xleft[ n ight] = delta left[ n ight]$，此时系统的输出就是所谓的单位脉冲响应$hleft[ n ight]$。再与放大器的例子对比一下：输入信号幅度$i = 1$，输出信号的幅度$o$就是放大倍数(放大器的重要指标，从应用者角度而言，其实就是唯一关心的指标)。我们有理由相信，对于LCCDE描述的LTI系统，得到了单位脉冲响应$hleft[ n ight]$，就能够掌握系统的全部特性(因果性、稳定性、频率选择性等)。何以见得呢？之后将逐步解答。
现在，我们假设已经知道了LTI系统的单位脉冲响应$hleft[ n ight]$，对于任意输入$xleft[ n ight]$，如何求系统的输出$yleft[ n ight]$呢？其实就是求滤波后的信号。首先，我们要建立已知和需求之间的联系。已知的是 ，需求是 。
图 11 单位脉冲响应
事实上，任意序列$xleft[ n ight]$很容易用$delta left[ n ight]$的移位、加权的线性组合来表示：$xleft[ n ight] = sumlimits_{k = - infty }^{ + infty } {xleft[ k ight]delta left[ {n - k} ight]} $。
图 12 序列的单位脉冲表示
结合系统的线性和时不变性，有
式中最后一个等式[yleft[ n ight] = sumlimits_{k = - infty }^{ + infty } {xleft[ k ight]hleft[ {n - k} ight]} ]就是卷积和，记作$yleft[ n ight] = xleft[ n ight]*hleft[ n ight]$，表明LTI系统的(零状态)响应$yleft[ n ight]$是单位脉冲响应$hleft[ n ight]$的移位、加权、求和。这句话也描述了求$yleft[ n ight]$的步骤。根据图 11，我们可以假想有这样的一个LTI系统——出钞机：如果你今天往出钞机投币口(输入端)投币一元，它将在接下来的三天，每一天都会从出币口(输出端)吐出一元。试问：如果某个月的1、3、4日，你分别投了2、2、3元，那么，出钞机每天的输出是多少元呢？不妨用图来说明一下。
图 13 卷积和的图解
图解结果告诉我们，出钞机在2-7日分别吐出2,2,4,5,5,3元。图解的过程包含了乘积与求和(“积”与“和”)，但并有体现出“卷”(翻转)这一操作。如果只想求某一天(比如5日)出钞口吐出多少钱，此时就要用另一种方法，即许多教材中描述的步骤：
(1) 变量替换：$xleft[ n ight] o xleft[ k ight],hleft[ n ight] o hleft[ k ight]$
(2) 将$xleft[ k ight]$或$yleft[ k ight]$翻转
(3) 移位、相乘、相加
以上步骤包含了“卷积和”所有的操作(卷——翻转；积——相乘；和——相加)。
图解过程是根据线性时不变系统的定义导出的一种结果，是系统特性的直接反映。卷积和是信号与系统、数字信号处理中最重要的公式之一，它描述了LTI系统输入$xleft[ n ight]$、输出$yleft[ n ight]$以及单位脉冲响应$hleft[ n ight]$之间的关系，已知三者中的任何两个，就可以确定第三者，于是就有以下应用场合：
(1) 已知输入$xleft[ n ight]$和系统单位脉冲响应$hleft[ n ight]$，求输出$yleft[ n ight]$，即用确定的系统对输入信号滤波处理；
(2) 已知输出$yleft[ n ight]$和系统单位脉冲响应$hleft[ n ight]$，求输入$xleft[ n ight]$；
(3) 已知输入$xleft[ n ight]$和输出$yleft[ n ight]$，求系统单位脉冲响应$hleft[ n ight]$。
现在，我们将注意力集中到卷积和公式上来，公式重写如下
[yleft[ n ight] = sumlimits_{k = - infty }^{ + infty } {xleft[ k ight]hleft[ {n - k} ight]} ]
卷积和只包含了移位、数乘和相加运算，看来，数字信号处理的计算十分简单。如果再上升一个层次，卷积和能很直观解释系统如何改变输入信号的吗？从图 13的结果很难看出系统是如何将$xleft[ n ight]$变成$yleft[ n ight]$的，而且也很难说明这种变化意味着什么，因此，卷积和在解释系统行为特性时很不直观，这是卷积和运算本身决定的。
向量运算的启示
了解运算背后的物理或几何含义对理解问题的本质至关重要。例如向量的乘法和加法：设两个向量$vec A = {x_1} + i{y_1}$，$vec B = {x_2} + i{y_2}$，向量的加法符合平行四边形法则，从图 14很容易看出两个向量相加的几何意义，而相应的代数运算也很简单——对应分量相加即可，即有
$vec C = vec A + vec B = left( {{x_1} + {x_2}} ight) + ileft( {{y_1} + {y_2}} ight)$
但如果要求两个向量的乘法呢？代数运算的结果是
相比两个向量的加法，向量的乘法运算要复杂一些，更重要的是，根据运算结果，我们很难看清向量乘法的几何意义。如果将向量$vec A$和$vec B$用极坐标表示，情况又会如何呢？设
则两个向量的乘法
计算简单，几何意义明确——模相乘，相位相加。
图 14 向量加的直角坐标表示
图 15 向量乘法的极坐标表示
向量运算的例子告诉我们，选择不同的坐标系，不仅影响运算的复杂度，在解释运算的几何意义时也各有千秋。上例告诉我们，如果是向量的加法，直角坐标系比极坐标系方便；如果是向量的乘法，很显然，极坐标系比直角坐标系方便。
回到数字信号处理的话题，卷积和是让无数人困惑的公式，然而，它又是经典数字信号处理算法的根源。既然卷积和运算就是滤波，我们如何评判滤波的效果呢？从什么角度理解系统的滤波行为更好呢？
再想想知了的叫声，之所以烦人，是因为包含了大量的高频分量。看来，用频率这一参数来解释信号特性很符合人的直观感觉。我们有更好的方式来洞悉卷积和的物理意义吗？能不能站在频(率)域的角度来解释卷积和呢？当然可以，因为我们有卷积定理(性质)：
通过傅里叶变换，将时域中的卷积和运算换成了频域中的相乘运算。那么频域中的相乘运算有什么好处呢？不要忘了，我们主要是为了直观解释滤波的行为和特性。
图 16是一段滤波前知了声$xleft[ n ight]$的时域波形，单从图形上看，似乎很难说出这一段信号的特点，运用傅里叶变换，得到$xleft[ n ight]$的幅度谱$left| {Xleft( {{e^{jomega }}} ight)} ight|$如图 17所示，很明显，该图说明信号$xleft[ n ight]$包含了丰富的高频成分。现在，我们希望设计一个滤波器，可以滤除2500Hz以上的频率分量，目的是仅保留红框内的频率分量，因此要求滤波器的频率响应(幅频响应)是
滤波器的幅频特性如图 18所示，截止频率${f_c} = 2500$Hz。由卷积定理可知，输出信号的幅度谱为
[left| {Yleft( {{e^{jomega }}} ight)} ight| = left| {Xleft( {{e^{jomega }}} ight)} ight|left| {Hleft( {{e^{jomega }}} ight)} ight|]
从图形上来看就更直观了：图 17与图 18相乘，得到了滤波后信号的幅度谱，滤波前、后信号频谱的变化一目了然，很显然，它是按照滤波器频率响应$Hleft( {{e^{jomega }}} ight)$所规定的形状去改变的。而对比滤波前、后的时域波形，很难解释信号的什么特性被改变了。
在上述的讨论中，我们只关注幅度谱，实际中，相位谱也必须考虑，更完整的表达式为
对比极坐标系下两个向量的乘法，能找到相似之处吗？——模相乘、相位相加。如果将卷积和$yleft[ n ight] = xleft[ n ight]*hleft[ n ight]$看作直角坐标系下两个向量的乘法，那么[Yleft( {{e^{jomega }}} ight) = Xleft( {{e^{jomega }}} ight)Hleft( {{e^{jomega }}} ight)]就类似于极坐标系下两个向量的乘法——运算简单、物理意义直观。[通过向量乘法运算的例子，我们可以看到，向量的表示方式起到了关键性的作用；同样地，信号也有不同的描述方式，如时域、频域、复频域，不同的表示方法都是为了使分析问题更简单、物理意义更明确]
图 16滤波前知了声$xleft[ n ight]$时域波形
图 17滤波前知了声的幅度谱$left| {Xleft( {{e^{jomega }}} ight)} ight|$ (红框是希望保留的低频分量)
图 18低通滤波器幅频特性$left| {Hleft( {{e^{jomega }}} ight)} ight|$ (蓝色线)
图 19滤波后知了声的幅度谱$left| {Yleft( {{e^{jomega }}} ight)} ight|$
图 20滤波后知了声$yleft[ n ight]$时域波形
滤波器设计征程
借助系统的频率响应$Hleft( {{e^{jomega }}} ight)$，可以很方便地解释滤波器的特性(主要是指频率选择性)。有什么定量指标来描述频率选择性吗？通常有四个指标——横坐标两个，纵坐标两个。以低通滤波器为例来说明。如图 21所示，横坐标两个指标为通带截止频率${omega _p}$和阻带截止频率${omega _s}$，纵坐标两个指标为通带衰减${A_p}$和阻带衰减${A_s}$。
图 21数字低通滤波器指标
通常而言，滤波器设计就是指根据应用的需求，提出合理的滤波器指标(通带截止频率和衰减，阻带截止频率和衰减)，然后借助滤波器设计工具得到滤波器参数(系数)。滤波器系数？没听说过啊，是什么东西？其实我们早见过面了，就是差分方程中的${a_k}$和${b_k}$!早就说过，数字滤波器设计要不忘初心，根本任务就是确定这三个参数——方程的阶数$N$，系数${a_k}$和${b_k}$。
如何建立起滤波器指标和滤波器系数之间的联系呢？很显然，滤波器的指标约束着其频率响应$Hleft( {{e^{jomega }}} ight)$，而$Hleft( {{e^{jomega }}} ight)$与LCCDE的阶数$N$，系数${a_k}$和${b_k}$直接相关。傅里叶变换又该出场了。对差分方程两边取傅里叶变换得
[sumlimits_{k = 0}^N {{a_k}{e^{ - jkomega }}Yleft( {{e^{jomega }}} ight)} = sumlimits_{k = 0}^N {{b_k}{e^{ - jkomega }}Xleft( {{e^{jomega }}} ight)} ]
从而可得频率响应
[Hleft( {{e^{jomega }}} ight) = frac{{Yleft( {{e^{jomega }}} ight)}}{{Xleft( {{e^{jomega }}} ight)}} = frac{{sumlimits_{k = 0}^N {{b_k}{e^{ - jkomega }}} }}{{sumlimits_{k = 0}^N {{a_k}{e^{ - jkomega }}} }}]
易见，频率响应$Hleft( {{e^{jomega }}} ight)$完全由$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。
乍一看，这些指标与图 21中的指标完全不同，实际需求的指标往往同本例，频率指标为模拟频率，衰减指标以dB为单位，而设计数字滤波器时，往往需要数字频率指标。如何进行指标转换呢？直接用公式计算${omega _p} = 2pi {F_p} = 4000pi $，${omega _s} = 2pi {F_s} = 5000pi $吗？图 21告诉我们，无论通带截止频率，还是阻带截止频率都不会超过$pi $。难道最大的数字角频率就是$pi $吗？为什么？这个问题同样令许多人困惑。如何衡量数字频率的大小呢？先考虑一个简单的数字信号：$xleft[ n ight] = cos left( {omega n} ight)$，针对$omega $取不同的值，分别画出对应的波形如下：
图 22$cos left( {omega n} ight)$波形图
振荡得越厉害，意味着频率越高，这是符合常理的。上图中的9个波形，哪个振荡得最厉害呢？第二排中间那个——+1和-1交替，而此时$omega = pi $。而且还发现，$omega = 0$和$omega = 2pi $的波形一样，并不是$omega $的值越大，振荡得越快。如果是连续时间信号$xleft( t ight) = cos left( {2pi ft} ight) = cos left( {Omega t} ight)$，大家都知道，随着$Omega $的增大，波形会振荡得越快。从连续时间信号到离散时间信号，是什么原因导致两者出现如此大的差异呢？通过对采样过程的探索，我们便可理解其中的奥妙。设以时间间隔${T_s}$对$xleft( t ight) = cos left( {Omega t} ight)$采样，这样便得到
$xleft( {n{T_s}} ight) = cos left( {Omega n{T_s}} ight) = cos left( {omega n} ight)$
因此有
$omega = Omega {T_s} = Omega /{f_s}$
其中${f_s} = 1/{T_s}$为采样频率，$Omega $为模拟角频率，$omega $为数字角频率。结论是：数字频率是模拟频率关于采样频率的归一化。
根据Nyquist采样定理，当采样频率大于信号最高频率2倍时，可以由样点完全恢复原来的信号。我们不妨假设信号的最高频率正好是采样频率的一半，即[{f_{max }} = {f_s}/2]，对应的最高模拟角频率${Omega _{max }} = 2pi {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。将模拟频率指标代入式计算出对应的归一化数字频率指标
这样，我们就得到了数字域的四个指标：${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}$确定之后，方程也就完全确定了，所谓的滤波，就是要对给定的输入$xleft[ n ight]$，求系统的响应$yleft[ n ight]$。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) 系统的描述方法(LCCDE、单位脉冲响应、频率响应、系统函数、零极点图)
(3) 滤波器设计方法(代码和可视化方法)
(4) 滤波的实施方法(滤波器类型的选择——IIR、FIR；滤波器结构的选择，量化对系数的影响)