跟小D每日学口语

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

信号分析与处理概述

博客对图片和公式的支持不是很好,相应的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;滤波器结构的选择,量化对系数的影响)

 

 

 

posted on 2016-07-11 18:37  湘厦人  阅读(4267)  评论(0编辑  收藏  举报

导航