MATLAB数字信号处理课件.ppt

1、第5章 使用MATLAB实现数字信号处理 本章主要内容如下：51 数字信号处理基本内容及相应的MATLAB工具 52 信号通过系统的时域分析 53 信号通过系统的频域和Z域分析 54 滤波器设计 55 频谱分析 5.1 数字信号处理基本内容及相应的MATLAB工具数字信号处理的基本内容通常分为两部分：离散时间信号与系统分析 主要涉及离散时间信号与系统的时域、频域表示，以 及信号通过系统的时域、频域分析及其变换域分析。MATLAB函数库中提供了filter,conv,convmtx,fft,ifft,freqz,impz,zplane等等与之相应的函数。数字滤波器设计和谱分析 数字滤波器设计包括

2、了无限冲激响应(IIR)和有限冲激响应（FIR）滤波器设计，谱分析又可进一步分为线性谱分析和非线性谱分析。MATLAB为此提供了多种成熟算法的相应函数以及极为丰富的设计工具。5.2 时域分析卷积，滤波，单位冲激响应5.2.1 卷积 MATLAB提供 conv函数实现标准的一维信号卷积:例如，若系统h(n)为h=1 1 1输入序列x(n)为x=1 1 1则x(n)经过系统h(n)后的MATLAB实现为：conv(h,x)或 conv(1 1 1,1 1 1)执行后即得到y(n)为ans=1 2 3 2 1注意：使用conv 函数时，h(n)和x(n)都必须是有限长的,，否则不能使用conv 函数

3、。例5-1 时域离散序列的卷积计算与图形显示 例5-1(教材p63)：已知离散信号x(n)和h(n)，求y(n)=x(n)*h(n)，并用图形表示。)()(nRnx101)()(5nRnx102)()(nRnx101)(.)(nR90nh20n1)(.)(nR90nh20n2)()()(nhnxny111)()()(nhnxny222例5-1的MATLAB程序Nh=20;Nx=10;m=5;%设定Nx，Nh和位移值mn=0:Nh-1;h1=(0.9).n;%产生h1(n)h2=h1;nx=0:Nx-1;x1=ones(1,Nx);%产生x1(n)x2=zeros(1,Nx+m);for k=m

4、+1:m+Nx%产生x2(n)=x1(n-m)x2(k)=x1(k-m);end%产生x2(n)y1=conv(x1,h1);%计算y1(n)=x1(n)*h1(n)y2=conv(x2,h2);%计算y2(n)=x2(n)*h2(n)subplot(3,2,1)stem(nx,x1,.)axis(0 30 0 1.2),title(x1(n)%绘图（以下省略）5.2.2 滤波 数字滤波器的系统函数H(z)用如下式表示:在MATLAB中，用向量b，a来表示滤波器的系数b(i)和 a(i)。m1n1z1maz2a1az1nbz2b1bzH)()()()()()()(滤波器分类 当n=0，m0时，

5、称为AR滤波器，即自回归（Auto Recurrence）滤波器，具无限冲激响应（IIR），也即其单位采样响应h(n)具无限长度；若m=0，a(1)0，称为MA滤波器，即滑动平均（Moving Average）滤波器，其单位采样响应h(n)是有限长度，故称有限冲激响应（FIR）滤波器；如果n、m都大于零，称为ARMA滤波器，而其冲激响应也为IIR。filter函数 MATLAB提供了 filter函数来对离散信号进行滤波，表达信号通过系统后的结果。与conv不同的是，filter函数可适用于无限冲激响应系统的情况，但信号仍须是有限长的。例如，一个单极点的低通滤波器系数如下：b=1;%分子系数向

6、量b(i)a=1-0.9;%分母系数向量a(i)如果用filter函数实现对信号x滤波，只要调用：y=filter(b,a,x);就可给出输入x经过滤波以后的输出y。5.2.3 单位冲激响应 数字滤波器的单位冲激响应定义为输入为单位样本序列时数字滤波器的响应，即：h(n)=T(n)其中：1n01n1)n(单位冲激响应的MATLAB实现 MATLAB近似实现单位采样信号的方法为：imp=1;zeros(p,1);%zeros(p,1)产生 p个零元素组成的列向量，p是正整数。使用imp后，滤波器的冲激响应可近似得到为：h=filter(b,a,imp);impz 函数可以直接求出数字滤波器的单位

7、冲激响应，即：impz(b,a)该命令将同时绘出滤波器的单位冲激响应，教材p66图5-2。53 频域和Z域分析 频率响应，零极点分析 5.3.1 频率响应 MATLAB数字信号处理工具箱有很多函数提供对模拟和数字滤波器的频率响应分析。其中，freqz 函数和freqs 函数分别返回数字和模拟滤波器的频率响应。工具箱中通常使用的单位频率是Nyquist频率，即采样频率的1/2。注意：就数字滤波器函数来说，其频域指标中的所有频率都以Nyquist频率进行归一化。因此Nyquist频率也称归一化频率。关于Nyquist频率的说明 例如例如：系统采样频率为1000Hz，则若数字滤波器的截止频率等于30

8、0Hz，经Nyquist频率归一化后，其归一化频率就是300/500=0.6。若将归一化频率转换成数字信号处理教科书中所使用的数字频率(rad)，需乘以；反之，若乘以采样频率fs的一半，则将归一化频率转换回了模拟域频率(Hz)。归一化频率f 应满足0f1。对于频率响应而言，归一化频率和模拟频率都可使用。freqz命令 freqz 函数提供了基于FFT算法所得到的数字滤波器H(z)的频率响应 最常用的形式：h,w=freqz(b,a,p)其意义是，对于数字滤波器：freqz 函数接受滤波器系数向量b和a，以及一个用于指定所需计算的频率响应样点数的整数p，返回一个与由p个频率样本点构成的实频率向量

9、相对应的复频率响应向量h。m1n1z)1m(az)2(a)1(az)1n(bz)2(b)1(b)z(Hfreqz的命令形式 h=freqz(b,a,w)采用上面的形式时，需先对频率样本点向量w作出定义。通常的做法是使用linspace函数。h,w=freqz(b,a)w和p没有定义，默认w由（0）上均分的512点构成，频率单位为rad/sample。h,w=freqz(b,a,p，whole)使用带参数whole选项的命令形式 h,w=freqz(b,a,p,fs)用参数选项指定采样频率fs linspace的命令形式 linspace函数的命令形式有两种：linspacelinspace（x

10、1,x2）%产生x1,x2范围内经线性分割得到的100点的向量linspacelinspace（x1,x2，N）%产生x1,x2范围内经线性分割得到的N点的向量例如：例如：w=linspace（0，pi）%w由（0，pi）上100个等分点组成 w=linspace（0，1000）%w由（0，1000）上100个等分点组成freqz函数无返回输出参数格式的调用 freqz函数还能以无返回输出参数的格式调用：freqz（b，a，256）或 freqz（b，a，256，2000）无返回输出参数调用freqz函数的好处是，MATLAB能够自动绘出频率在（0）范围内的幅频特性和相频特性图。注意，无返回输

11、出参数调用freqz函数绘出的相频特性不能正确给定在 处的值，这是因为在MATLAB中，属于下半个单位圆。频率响应的实例 例：先构成一个截止频率为400Hz的9阶巴特沃思（Butterworth）低通数字滤波器，求出其系数b，a，再求出其256点频率响应。指定的采样频率fs=2000Hz。实现1：先调用butter函数，再调用freqz函数；实现2：无返回输出参数，调用freqz函数；实现3：为得到完整的（0）范围内的幅相特性图，可用带返回输出参数的格式调用freqz，在得到各频率点上的频率响应的数值后，使用MATLAB所提供的abs函数和angle函数分别提取幅频特性与相频特性。注意：为了得

12、到连续的相频曲线，需要使用unwrap函数，该函数通过在发生跳变后的各处自动加上，从而除去了相位的跳变，这称为相位的“解卷绕”。实现1b,a=butter(9,400/1000);%括号中第一个输入是阶数，第二个是归一化截止频率h,f=freqz(b,a,256,2000);或h=freqz（b，a，256，2000）；实现2 实现1中第二条命令的形式可改写为 freqz（b，a，256）或freqz（b，a，256，2000）可自动绘出频率在（0）范围内的幅频特性和相频特性图。注意无返回输出参数调用freqz函数绘出的相频特性不能正确给定在=处的值，因为=属于下半个单位圆。实现3b,a=bu

13、tter(9,400/1000);w=linspace(0,1000);h=freqz(b,a,w,2000);mag=abs(h);pha=angle(h);%提取滤波器频率响应的幅度mag和相位phasemilogy(w,mag);title(Magnitude);figure;plot(w,pha);title(Phase);出现了幅度为2的跳变卷绕解卷绕：unwrap函数 为了得到连续的相频曲线，需要使用unwrap函数，该函数通过在发生2跳变后的各处自动加上2，从而除去了相位的跳变，这称为相位的“解卷绕”。解卷绕freqs函数 freqs函数用于求取由b和a系数向量定义的模拟滤波器H

14、（s）的频率响应。使用方法类似freqz函数。与第二章（p32例2-1）采用数组相除方法求取频率响应相比，使用freqs 函数要方便很多。5.3.2 零极点分析 zplane 函数用于画出线性系统在Z平面上的零极点。有两种使用方法：1、在已知零极点时，例如某滤波器的零点为-1/2，一对共轭极点为 和 时，只要输入命令zer=-0.5;pol=0.9*exp(j*2*pi*-0.3 0.3);zplane(zer,pol)即可画出零极点。（见p70图5-6）)3.0(2je9.0)3.0(2je9.0用zplane函数绘制零极点图 另一种情况：已知系统的系统函数系数向量b 和 a，则可通过调用

15、zplane(b,a)绘出零极点。这种情形下，zplane 函数先求得系统函数的零点和极点，然后绘出零极点图。5.4 滤波器设计滤波器技术指标，IIR滤波器设计，FIR滤波器设计5.4.1 滤波器的技术指标 设计滤波器之前，要定义一组滤波器的技术指标。滤波器的技术指标通常采用幅度指标给出，以保证物理上的可实现性。幅度指标以容差方式给出，如p71图5-7所示。具体指标：通带截止频率Wp；阻带截止频率Ws；通带波纹p；阻带衰减s。其中通带波纹和阻带衰减也可以用分贝数给出)1log(20pp)log(20ss滤波器以容差方式给出的幅度指标|H(j)|10通带通带过渡带过渡带阻带阻带p s s p 1

16、 5.4.2 IIR 滤波器设计 MATLAB工具箱提供了几种模拟滤波器原型的产生函数,如巴特沃思(Butterworth),切比雪夫(Chebyshev)滤波器等。还提供了从模拟低通滤波器原型转换为高通、带通和带阻的转换函数，模拟滤波器转换为数字滤波器的双线性变换法和冲激响应不变法，模拟IIR滤波器的阶数选择函数以及数字滤波器直接设计函数等等，使用起来非常方便。本节以巴特沃思滤波器为例 介绍 MATLAB工具箱提供的滤波器设计函数 使用方法。1巴特沃思滤波器设计函数 b,a=butter(n,Wn,options)缺省情况下，返回一个用有理分式表示的低通数字滤波器系统函数的分子分母系数向量b

17、和a。设计的技术指标只需要指定一个归一化截止频率Wn（这里Wn的单位为：Nyquist 频率=1 Hz）。参数选项options：选项空缺返回低通滤波器；high则表示返回的是高通数字滤波器；至于带通和带阻滤波器设计，则输入的Wn必须是二个元素的向量；选项options 空缺返回带通滤波器，带阻滤波器需加参数options为stop。Butter函数举例 B,A=butter(N,Wn)返回系统函数的长度为N+1的分子系数B和分母系数A。式中Wn为归一化截止频率，0 Wn 1.0。B,A=butter(N,Wn,high)设计一个高通数字滤波器。如果Wn=W1 W2，且options空缺，则返

18、回一个2N阶的带通滤波器，其通带为：W1 W W2。B,A=butter(N,Wn,stop)，且Wn=W1 W2 设计2N阶的带阻滤波器。butter(N,Wn,s)，butter(N,Wn,high,s)或 butter(N,Wn,stop,s)调用butter()函数附加一个options为s时，表示设计的是一个模拟滤波器，而指定的截止频率Wn单位此时必须为rad/s。2.巴特沃思滤波器阶次选择函数 buttord函数用于给出满足给定频域指标要求的最小阶次滤波器的设计。其调用格式为：n,Wn=buttord(Wp,Ws,Rp,Rs)%用于数字滤波器 n,Wn=buttord(Wp,Ws,

19、Rp,Rs,s)%用于模拟滤波器 该函数返回符合技术指标的最小阶数N以及Butterworth滤波器3dB截止频率Wn。设计指标同前：通带起伏不超过Rp，阻带衰减不小于Rs，Wp和Ws分别为归一化通带和阻带截止频率。Buttord函数应用：例5-1 例：设计一个数字带通滤波器，通带为1000 到 2000 Hz，止带的起始位置离开通带两边500Hz，也即分别为500 Hz和2500Hz，设采样率为10KHz，通带起伏为1dB,止带最小衰减为60dB。解 MATLAB实现如下：n,Wn=buttord(1000 2000/5000,500 2500/5000,1,60)得到：n=12Wn=0.1

20、951 0.4080 因此b,a=butter(n,Wn);给出了满足要求的最小阶次巴特沃思带通滤波器。3经典的IIR滤波器设计方法 一种常用的IIR滤波器设计方法是利用已经很成熟的模拟滤波器设计结果。设计步骤如下：（1）按照一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。（2）根据转换后的技术指标设计模拟低通滤波器的H(s)。（3）再按一定规则将H(s)转换成H(z)。步骤(3)中最常用的方法:冲激响应不变法（也克称为脉冲响应不变法）和双线性变换法。冲激响应不变法的原理 冲激响应不变法是通过对模拟滤波器的单位冲激响应ha(nTs),进行采样，将其作为数字滤波器的单位冲激响

21、应h(n)。即：据此设计得到的数字滤波器的频率响应 与模拟滤波器的频率响应 之间的关系为：因此存在着频率混迭。从而限制了其适用范围。)nT(h)n(hsa)Tk2jTj(HT1)e(Hsksasj)e(Hj)j(Ha注意：S域频率轴与Z域单位圆之间不存在一一对应关系。冲激响应不变设计的impinvar()函数 impinvar()函数的调用格式是：Bz,Az=impinvar(Bs,As,fs)这一函数把具有Bs,As模拟滤波器的转换成了采样频率为fs的数字滤波器的Bz,Az。如果没有给定采样频率fs，函数默认为1Hz。双线性变换法原理 双线性变换法所定义的s平面与z平面的映射关系为：给定模拟

22、滤波器的系统函数，数字滤波器为：因此S域内的有理分式函数在Z域内仍然保持了有理分式函数形式，从而保证了数字滤波器的物理可实现性。11sz1z1T2s 11sz1z1T2sasHzH)(双线性变换法原理（续）双线性变换将S平面上的频率轴一对一映射到了Z平面的单位圆上，消除了频率混叠问题：但Z平面单位圆和S域频率轴之间存在严重的非线性关系；由此也对其适用范围产生了限制；并且，在滤波器设计中，在频域指标转换时需要进行Pre-warping。)2T(tan2s1双线性变换设计的bilinear函数 bilinear()函数的调用格式有三种：Bz Az=bilinear(Bs,As,Fs)%把具有Bs,

23、As模拟滤波器的转换成数字滤波器的Bz,Az。Fs是采样频率。Zd Pd Kd=bilinear(Z,P,K,Fs)%把模拟滤波器的零点、极点，转换成数字滤波器的零点、极点模型，Fs是采样频率。Ad Bd Cd Dd=bilinear(A,B,C,D,Fs)%把模拟滤波器的状态方程模型转换成数字滤波器的状态方程模型。Fs是采样频率。注意：这里的采样频率均为实际频率。注意：这里的采样频率均为实际频率。例5-3 IIR数字滤波器的设计示例 例5-3中直接给出了模拟低通滤波器技术指标（）。如果给定的是数字低通滤波器技术指标，则应先将数字滤波器的技术指标折合成相应的模拟低通滤波器指标，然后采用示例所用

24、的设计方法进行滤波器设计。|)j(H|a例5-3 设计要求与指标（1）设计一个巴特沃思模拟低通滤波器，指标如下：通带内 ，波纹小于10dB；阻带内 ，衰减不小于20dB。（2）用双线性变换法，将上述模拟低通滤波器变换为数字低通滤波器，采样间隔分别取T=1s，1.5s，比较设计结果。（3）用冲激响应不变法，将上述的模拟低通滤波器变换为数字低通滤波器，采样间隔分别取T=1s，1.5s，比较设计结果。s/rad2.0ps/rad3.0s例5-3的解 由于给定的设计指标是模拟滤波器的频域指标，所以，首先选用buttord()函数进行最小阶数的滤波器设计。N,Wn=buttord(Wp,Ws,Rp,Rs

25、,s)注意这里Wp，Ws单位为rad/s。此题中，Wp=0.2*pi,Ws=0.3*pi,Rp=10,Rs=20，所以可用N,Wn=buttord(0.2*pi,0.3*pi,10,20,s);得到实现该频域指标要求的巴特沃思低通滤波 器最低阶次N和3dB截止频率Wn。例5-3的解（续）再使用butter()函数，得到模拟低通滤波器的系统函数的系数向量Bs和As：Bs,As=butter(N,Wn,s);根据采样间隔分别为T=1s，1.5s，可得采样频率分别为：Fs1=1Hz，Fs2=2/3Hz，调用bilinear()函数和impinvar()函数，可以得到相应的两种数字滤波器系统函数的系数

26、向量Bz和Az。即：Bz,Az=bilinear(Bs，As，Fs)或 Bz,Az=impinvar(Bs，As，Fs)注意，由于低通巴特沃思滤波器并不是带限的，因此对其用冲激响应不变法设计数字滤波器时，将会发生混叠。（参见下张ppt）出现混叠无混叠5.4.3 FIR滤波器设计 FIR滤波器具有IIR滤波器所没有的线性相位特性，所以在工程中得到了广泛的应用。MATLAB中除了cremez函数以外，所有用于FIR滤波器设计的函数都是用于设计线性相位的FIR滤波器的。FIR滤波器的线性相位条件 FIR数字滤波器的系统函数为：线性相位条件：FIR数字滤波器的h(n)为实数，且满足以下任一条件 h(n

27、)=h(N-1-n)（关于n=(N-1)/2偶对称）h(n)-h(N-1-n)（关于n=(N-1)/2奇对称）nNnddznhzH10)()(根据h(n)的对称情况，分四种类型h(n)奇对称，N为奇数 h(n)奇对称，N为偶数 h(n)偶对称，N为奇数 h(n)偶对称，N为偶数 FIR线性相位滤波器的频率响应 具有线性相位的FIR滤波器的频率响应函数可以用下式表示：式中 或 ，。注意，式中N为脉冲响应h(n)的长度，是振幅（Amplitude）响应函数，不是幅度（Magnitude）响应函数 ，可以有正有负，与其相关的相位响应也是一个连续线性函数。对于前述四种h(n),其幅度函数有不同的特点。

28、)(jrje)(H)e(H0221N)(Hr|)jH(e|FIR滤波器设计：窗函数法 基本原理基本原理：将理想低通滤波器的无限长冲激响应序列截断（等效于加矩形窗）后乘以窗函数，用得到的有限长序列去逼近理想低通滤波器。设所要求的理想数字滤波器的频率响应为 ,是与其对应的单位脉冲响应，因此 由于 是矩形频率特性，故 一定是无限长的非因果序列。而所要设计的是有限长的FIR滤波器，因此用有限长的 去逼近无限长的 最直接的方法就是截断，即用有限长的窗函数 来截取 得到：这就是窗函数设计方法。)e(Hjd)(nhddeHnhnjdd)(21)()e(Hjd)(nhd)(nh)(nhd)(nw)()()(n

29、hnwnhd)(nhd窗函数法设计步骤（1）设所要求的FIR滤波器理想频率响应为：（2）由性能指标确定窗函数W(n)和窗口长度N，由过渡带宽近似于窗口函数主瓣宽求得窗口长度N；（3）由 求得实际滤波器的单位脉冲响应h(n)，即求取所设计的FIR滤波器h(n)的系数向量；（4）检验滤波器的性能。njndjde)n(h)e(H)n(W)n(h)n(hdFIR滤波器设计：fir1函数 MATLAB中使用窗函数法设计线性相位FIR滤波器的常用设计函数为fir1()，调用格式如下：b=fir1(n,Wn)b=fir1(n,Wn,ftype)b=fir1(n,Wn,window)b=fir1(n,Wn,f

30、type,window)fir1命令格式的说明 fir1命令中n为FIR滤波器的阶数，对于高通、带阻滤波器，n应该取偶数；Wn为滤波器的截止频率，0Wn1,对于带通、带阻滤波器，Wn=W1 W2，且满足W1W2,对于多带滤波器，Wn=W1 W2 W3，且满足0W1W21；b为FIR滤波器的系数向量，长度为n+1；ftype 是滤波器的类型，默认时为低通或带通滤波器，high 为高通滤波器，stop 为带阻滤波器；Window为窗函数，是长度为n+1的列向量。例5-4：低通滤波器 例5-4 理想数字低通滤波器的截止频率为0（rad），在频率小于0处，其幅度为1，在频率为0至处，幅度为0，故它的冲

31、激响应序列h(n)为：其冲激响应为无限长，因此不是因果可实现的。用窗截取其有限长度，再加以因果化，就可以得到一个可实现的线性相位FIR 滤波器。)n(csinde21de)(H21)n(h0o00njnj例5-4的解 例如取长度51，低通截止频率0=0.4，则得到:b=0.4*sinc(0.4*(-25:25);这里用了最简单的矩形窗。下面的命令将显示这个滤波器的频率响应（见p77图5-9）：H,w=freqz(b,1,512,2);%b为H(z)分子系数，1为分母系数，512个 样本，采样频率2Hzplot(w,abs(H);gridtitle(Truncated Sinc Lowpass

32、FIR Filter);图5-9 窗函数法低通滤波器的频率特性 使用其它窗函数 可以使用其他形式的窗函数；例如改用长度51的海明窗，即对H(z)的系数乘以海明窗函数，则命令为：b=0.4*sinc(0.4*(-25:25);b=b.*hamming(51);%对H(z)的系数乘以海明窗函数 H,w=freqz(b,1,512,2);plot(w,abs(H),grid plot(w,20*log10(abs(H),title(Hamming-Windowed Truncated Sinc LP FIR Filter);下一张图将给出其幅频特性，注意与矩形窗情况比较。图5-10 使用海明窗得到的

33、低通滤波器特性直接使用FIR滤波器设计函数fir1()b=fir1(51,0.4);%fir1()中窗函数（window）项空缺代表使用海明窗H,W=freqz(b,1,512);plot(W/pi,abs(H);gridtitle(Hamming_windowed FIR LPF)fir2()函数 函数 fir2()也是FIR 滤波器设计函数，用于设计具有任意形状频率响应的FIR 数字滤波器。例如：n=50;f=0 0.4 0.5 1;%给定通带、止带范围分别为0 0.4和0.5 1，1为Nyquist频率m=1 1 0 0;%给定通带，止带频率响应幅度要求b=fir2(n,f,m)将返回行

34、向量b 为包含 n+1 个系数的n 阶FIR 滤波器，其频率响应幅度符合给定的f 和m 要求。类似于fir1()函数，可以选择窗，命令形式为b=fir2(n,f,m,window)FIR滤波器设计：频率采样设计法 频率采样设计法是从频域出发，把给定的理想频率响应 等间隔采样N个点得到。再对 作IDFT，得到 然后将h（n）作为所设计的滤波器的单位冲激响应。)e(Hjd)k(H)k(H)e(HdN2jd)k(Hd)n(h1N0kknN2je)k(HN1)n(h1N,1,0k频率采样设计法的有关说明 对 进行频率采样得到的H(k)，可以表示成：为保证设计出的FIR滤波器具有线性相位，h(n)应在0

35、,N-1上具有对称性。由此进行的理论分析表明，频率采样值的幅度与相位也须满足一定的约束条件，分为4种类型。因此，频率采样法设计FIR滤波器可归结为从满足约束条件的H(k)中求出h(n)问题。)k(je)k(H)k(H)e(Hjd1N,1,0k频率采样设计法的MATLAB实现 从满足约束条件的H(k)中求出h(n)，可利用MATLAB中的ifft函数。注意：从时域角度看,求出的h(n)与欲设计的hd(n)是有差别的，两者关系如下：从频域角度看，由频率采样设计所得到的滤波器的频率特性是由H(k)内插形成的，因此对 的逼近程度与 特性的平滑性有关。rNd)n(R)rNn(h)n(h)e(Hjd)e(

36、Hjd存在时域混叠 例5-5：频率采样法设计举例 例5-5 所要设计的FIR低通数字滤波器技术指标为：3.0p，dB5Rp4.0s，dB40Rs例5-5的解20,12,11k)k21(212010,1,0kk2120)k(1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1)k(H在给定的图5-12理想低通滤波器特性上选择N=21点进行采样得到|H(k)|，如图5-13。附加相位约束条件)k(je)k(H)k(H对H(k)求反变换IDFT就可得h(n)。图5-14 例5-5频率采样法设计低通滤波器(N=21)最小阻带衰减仅为14dB 例5-5的改进设计 分析分析：

37、从图5-14幅度响应曲线可见，N=21时，最小的阻带衰减为14dB，这在绝大多数情况下是不能接受的。但仅仅增大采样点数，比如N=61时，最小的阻带衰减仍为18dB，仍然不能满足通常的止带指标要求。为此，除了增大N以外，还应采用在过渡带插入过渡采样点的办法来除去原始 的特性中的突变。|)e(H|jd图5-16 例5-5的改进设计(N=61，增加两个过渡采样点)增加两个过渡采样点后，在N=61时，最小的阻带衰减为45dB，已满足通常的止带衰减要求 5.5 频谱分析 频谱分析的主要作用是对信号的频率构成进行分析，了解其在频带内的分布情况；本课程内仅对确定性信号进行，不涉及随机信号的频谱分析；因此主要

38、工具为FFT，不涉及非线性谱分析或即参量谱分析的内容。离散傅里叶变换（DFT）在离散信号与系统分析中，DFT有着频谱分析仪的特性。对连续信号与系统通过时域采样，也可以使用DFT进行谱分析。因此，DFT在频谱分析中起了核心的作用。另一方面，从技术上而言，DFT运算存在着高效的快速算法FFT，因此信号的频谱分析中将主要使用MATLAB的fft函数以及ifft函数。使用DFT对连续信号进行谱分析 工程实际应用中，常会遇到连续时间信号，其频谱也是连续函数。为了利用DFT对信号进行频谱分析，须先对 进行时域采样，得到 ，再对 进行DFT，得到的 则是 的傅里叶变换 在频域区间 上的N点等间隔采样。对带限

39、信号来说，X(k)与 在 范围内 上的采样值除相差一个因子外完全相同。以上说明构成了用DFT对连续信号进行谱分析的理论基础。)t(xa)n(x)nT(x)n(xa)k(X)e(Xj2,0)j(Xa)2/,0(sNT2k)n(x用DFT进行连续信号频谱分析时注意点（1）若信号持续时间有限长，则其频谱无限宽。（2）若信号频谱有限，则其信号持续时间必定无限长。（3）实际应用中，所能观察到的信号长度总是有限的。（4）因此，用DFT对连续时间信号进行频谱分析时，采样后产生的频谱混叠失真是不可避的。（5）得到的谱分析结果总是近似的，近似的程度取决于信号带宽、采样率以及信号的截取长度。例5-6：连续信号的频

40、谱分析 已知模拟信号 对x(t)进行采样，得 ，分别以T=0.5s，0.1s进行采样，样本数分别取N=20，100，以使被分析信号等长，用DFT计算 的 并与 比较。0t00t5.00te)t(xtnTt)t(x)n(x)k(XN19k0NT2k)j(XT1例5-6的分析 原信号x(t)的傅立叶变换为 其幅度谱为 ，幅频特性如图5-17所示。由图可知，|X(j)|的-3dB点为=1，而=2时已衰减到-16dB。因此，如取采样频率s=4，造成的频率混叠效应会比较小，此时相当于fs=2,T=0.5s。在T=0.1s时，混迭将更小。j11dte)t(x)j(Xtj211)j(X例5-6解:x(t)采

41、样后所得信号的N点DFT1Nn0)n(21e1e)n(R)rNn(21)rNn(ue)n(xNTnTNrT)rNn(N1Nk021ee11e)n(x)k(XkN2jT1N0nknN2jNN按DFT定义求解题给连续信号x(t)经采样后x(n)的N点DFT实际情况下，信号表达式事先难以确知，因此解析式通常是无法得到的。两种采样率所得结果的比较例5-6的结果分析 实际应用情况下，通常是截取连续信号的一段样本进行DFT分析。显然，由于xa(t)本身非带限，且又仅截取了0,T0的一段，因此，采样后一定有频率混叠（aliasing）存在，再加上截取后的样本值x(n)并不等于前面式中所示的xN(n)，所以，

42、x(n)在数字频率0,内的DFT将不可能与相应模拟频率内的Xa(j)完全相同。造成误差的两个来源中，时域中x(n)与xN(n)的误差在实际应用中无法消除，但由于采样率不够而引入的频率混叠可以通过减小T而得到减轻。比较例5-6在两种采样率下的结果，可以明显看到提高采样率的作用，但其代价为计算量的增加。其他说明 实际工作中，被分析数据是对信号进行截取有限长度得到的，这相当于被分析的是原信号乘上了一个窗函数。这种时域上的截断，在频域上就是所得频谱是原信号的频谱和窗函数频谱的卷积，若信号频率不位于DFT的谱线位置，将会产生频谱的泄漏。在对实际信号进行DFT分析时，泄漏效应是难以避免的，只能通过选择好的

43、窗函数来抑制其影响。用DFT分析信号时，要求被分析的样本序列在所截取的范围内可视为是稳定的。利用DFT进行频谱分析时，DFT的计算结果是样本序列以窗口时间为周期，进行延拓后所得周期波形的傅里叶级数的系数，因此与信号的真实频谱仍有差别。fft函数与ifft函数 MATLAB提供fft()函数来计算DFT，fft()函数是用机器语言而不是以MATLAB指令写成的，因此它的运行速度很快。fft()函数的命令形式：Y=fft(x)-计算信号x的快速离散傅里叶变换Y。当x为矩阵（多通道信号）时，计算x中每一列信号的离散傅立叶变换。Y=fft(x,n)-计算信号x的n点FFT，当x的长度大于n时，截断x，否则补零。MATLAB还提供ifft()函数用于FFT的反变换计算，命令格式类似fft()函数。例5-7 当信号受到噪声污染后，很难从时域中看出它包含的频率成分，这时通常可以用FFT来分析信号频率成分。例5-7中给出了一个仿真产生的受到均匀分布随机噪声干扰的信号，例中用fft()函数进行频谱分析来找出其主要频率成分，并通过filter()函数滤去噪声的干扰。所得结果分别示于图5-20与图5-21中。5.6 练习 课后完成p91p92练习3-5。