例说信号处理与滤波嚣设计
信号分析与处理概述
博客对图片和公式的支持不是很好,相应的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]} \]
如果仅局限于基础(经典)的数字信号处理算法研究,全部内容都是围绕这个线性常系数差分方程(Linear Constant-Coefficient Difference Equation,LCCDE)展开的。那么,我们需要研究LCCDE什么呢?很显然,我们可以将LCCDE看成一个系统,它描述了输入$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$。类似的概念推广到LCCDE描述的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$就是放大倍数(放大器的重要指标,从应用者角度而言,其实就是唯一关心的指标)。我们有理由相信,对于LCCDE描述的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]$的步骤。根据图 11,我们可以假想有这样的一个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]} \]
卷积和只包含了移位、数乘和相加运算,看来,数字信号处理的计算十分简单。如果再上升一个层次,卷积和能很直观解释系统如何改变输入信号的吗?从图 13的结果很难看出系统是如何将$x\left[ n \right]$变成$y\left[ n \right]$的,而且也很难说明这种变化意味着什么,因此,卷积和在解释系统行为特性时很不直观,这是卷积和运算本身决定的。
向量运算的启示
了解运算背后的物理或几何含义对理解问题的本质至关重要。例如向量的乘法和加法:设两个向量$\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}} \right) + i\left( {{y_1} + {y_2}} \right)$
但如果要求两个向量的乘法呢?代数运算的结果是
相比两个向量的加法,向量的乘法运算要复杂一些,更重要的是,根据运算结果,我们很难看清向量乘法的几何意义。如果将向量$\vec A$和$\vec B$用极坐标表示,情况又会如何呢?设
则两个向量的乘法
计算简单,几何意义明确——模相乘,相位相加。
图 14 向量加的直角坐标表示
图 15 向量乘法的极坐标表示
向量运算的例子告诉我们,选择不同的坐标系,不仅影响运算的复杂度,在解释运算的几何意义时也各有千秋。上例告诉我们,如果是向量的加法,直角坐标系比极坐标系方便;如果是向量的乘法,很显然,极坐标系比直角坐标系方便。
回到数字信号处理的话题,卷积和是让无数人困惑的公式,然而,它又是经典数字信号处理算法的根源。既然卷积和运算就是滤波,我们如何评判滤波的效果呢?从什么角度理解系统的滤波行为更好呢?
再想想知了的叫声,之所以烦人,是因为包含了大量的高频分量。看来,用频率这一参数来解释信号特性很符合人的直观感觉。我们有更好的方式来洞悉卷积和的物理意义吗?能不能站在频(率)域的角度来解释卷积和呢?当然可以,因为我们有卷积定理(性质):
通过傅里叶变换,将时域中的卷积和运算换成了频域中的相乘运算。那么频域中的相乘运算有什么好处呢?不要忘了,我们主要是为了直观解释滤波的行为和特性。
图 16是一段滤波前知了声$x\left[ n \right]$的时域波形,单从图形上看,似乎很难说出这一段信号的特点,运用傅里叶变换,得到$x\left[ n \right]$的幅度谱$\left| {X\left( {{e^{j\omega }}} \right)} \right|$如图 17所示,很明显,该图说明信号$x\left[ n \right]$包含了丰富的高频成分。现在,我们希望设计一个滤波器,可以滤除2500Hz以上的频率分量,目的是仅保留红框内的频率分量,因此要求滤波器的频率响应(幅频响应)是
滤波器的幅频特性如图 18所示,截止频率${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|\]
从图形上来看就更直观了:图 17与图 18相乘,得到了滤波后信号的幅度谱,滤波前、后信号频谱的变化一目了然,很显然,它是按照滤波器频率响应$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)$,可以很方便地解释滤波器的特性(主要是指频率选择性)。有什么定量指标来描述频率选择性吗?通常有四个指标——横坐标两个,纵坐标两个。以低通滤波器为例来说明。如图 21所示,横坐标两个指标为通带截止频率${\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)$与LCCDE的阶数$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。
乍一看,这些指标与图 21中的指标完全不同,实际需求的指标往往同本例,频率指标为模拟频率,衰减指标以dB为单位,而设计数字滤波器时,往往需要数字频率指标。如何进行指标转换呢?直接用公式计算${\omega _p} = 2\pi {F_p} = 4000\pi $,${\omega _s} = 2\pi {F_s} = 5000\pi $吗?图 21告诉我们,无论通带截止频率,还是阻带截止频率都不会超过$\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。将模拟频率指标代入式计算出对应的归一化数字频率指标
这样,我们就得到了数字域的四个指标:${\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) 系统的描述方法(LCCDE、单位脉冲响应、频率响应、系统函数、零极点图)
(3) 滤波器设计方法(代码和可视化方法)
(4) 滤波的实施方法(滤波器类型的选择——IIR、FIR;滤波器结构的选择,量化对系数的影响)