MATLAB数字信号处理课件

第5章使用MATLAB实现数字信号处理本章主要内容如下：5．1 数字信号处理基本内容及相应的MATLAB工具5．2 信号通过系统的时域分析5．3 信号通过系统的频域和Z域分析5．4 滤波器设计5．5 频谱分析第5章使用MATLAB实现数字信号处理本章主要内容如下：15.1数字信号处理基本内容

及相应的MATLAB工具数字信号处理的基本内容通常分为两部分：离散时间信号与系统分析主要涉及离散时间信号与系统的时域、频域表示，以及信号通过系统的时域、频域分析及其变换域分析。

MATLAB函数库中提供了filter,conv,convmtx,fft,ifft,freqz,impz,zplane等与之相应的函数。数字滤波器设计和谱分析数字滤波器设计包括了无限冲激响应(IIR)和有限冲激响应（FIR）滤波器设计，谱分析又可进一步分为线性谱分析和非线性谱分析。MATLAB为此提供了多种成熟算法的相应函数以及极为丰富的设计工具。5.1数字信号处理基本内容

及相应的MATLAB工具数字信25.2 时域分析卷积，滤波，单位冲激响应5.2 时域分析卷积，滤波，单位冲激响应35.2.1卷积MATLAB提供

conv函数实现标准的一维信号卷积:

例如，若系统h(n)为 >>h=[111]输入序列x(n)为 >>x=[111]则x(n)经过系统h(n)后的MATLAB实现为： >>conv(h,x)或

conv([111],[111])执行后即得到y(n)为 ans= 12321注意：使用conv函数时，h(n)和x(n)都必须是有限长的,，否则不能使用conv函数。5.2.1卷积MATLAB提供conv函数实现标准的一维4例5-1时域离散序列的卷积计算

与图形显示

例5-1(教材p63)

：已知离散信号x(n)和h(n)，求y(n)=x(n)*h(n)，并用图形表示。例5-1时域离散序列的卷积计算

与图形显示例5-1(教5例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);fork=m+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([03001.2]),title(‘x1(n)’) %绘图……（以下省略）例5-1的MATLAB程序Nh=20;Nx=10;m=5; 65.2.2滤波数字滤波器的系统函数H(z)用如下式表示:在MATLAB中，用向量b，a来表示滤波器的系数b(i)和

a(i)。

5.2.2滤波数字滤波器的系统函数H(z)用如下式表示:7滤波器分类当n

=

0，m≠0时，称为AR滤波器，即自回归（AutoRecurrence）滤波器，具无限冲激响应（IIR），也即其单位采样响应h(n)具无限长度；若m

=

0，a(1)≠0，称为MA滤波器，即滑动平均（MovingAverage）滤波器，其单位采样响应h(n)是有限长度，故称有限冲激响应（FIR）滤波器；如果n、m都大于零，称为ARMA滤波器，而其冲激响应也为IIR。滤波器分类当n

=

0，m≠0时，称为AR滤波器，即自回归（8filter函数

MATLAB提供了filter函数来对离散信号进行滤波，表达信号通过系统后的结果。与conv不同的是，filter函数可适用于无限冲激响应系统的情况，但信号仍须是有限长的。例如，一个单极点的低通滤波器系数如下：

>>b=1; %分子系数向量b(i) >>a=[1-0.9]; %分母系数向量a(i)

如果用filter函数实现对信号x滤波，只要调用： >>y=filter(b,a,x);

就可给出输入x经过滤波以后的输出y。filter函数MATLAB提供了filter函数来对离95.2.3单位冲激响应数字滤波器的单位冲激响应定义为输入为单位样本序列时数字滤波器的响应，即：h(n)=T[δ(n)] 其中：5.2.3单位冲激响应数字滤波器的单位冲激响应定义为输入10单位冲激响应的MATLAB实现MATLAB近似实现单位采样信号的方法为：

imp=[1;zeros(p,1)]; %zeros(p,1)产生p个零元素组成的列向量，p是正整数。使用imp后，滤波器的冲激响应可近似得到为：

h=filter(b,a,imp);impz函数可以直接求出数字滤波器的单位冲激响应，即：impz(b,a) 该命令将同时绘出滤波器的单位冲激响应，教材p66图5-2。单位冲激响应的MATLAB实现MATLAB近似实现单位采样信115．3频域和Z域分析

频率响应，零极点分析

5．3频域和Z域分析频率响应，零极点分析125.3.1频率响应MATLAB数字信号处理工具箱有很多函数提供对模拟和数字滤波器的频率响应分析。其中，freqz函数和freqs函数分别返回数字和模拟滤波器的频率响应。工具箱中通常使用的单位频率是Nyquist频率，即采样频率的1/2。注意：就数字滤波器函数来说，其频域指标中的所有频率都以Nyquist频率进行归一化。因此Nyquist频率也称归一化频率。5.3.1频率响应MATLAB数字信号处理工具箱有很多函13关于Nyquist频率的说明例如：系统采样频率为1000Hz，则若数字滤波器的截止频率等于300Hz，经Nyquist频率归一化后，其归一化频率就是300/500=0.6。若将归一化频率转换成数字信号处理教科书中所使用的数字频率(rad)，需乘以π；反之，若乘以采样频率fs的一半，则将归一化频率转换回了模拟域频率(Hz)。归一化频率f应满足0<f<1。对于频率响应而言，归一化频率和模拟频率都可使用。关于Nyquist频率的说明例如：系统采样频率为1000Hz14freqz命令

freqz

函数提供了基于FFT算法所得到的数字滤波器H(z)的频率响应最常用的形式：[h,w]=freqz(b,a,p)

其意义是，对于数字滤波器： freqz函数接受滤波器系数向量b和a，以及一个用于指定所需计算的频率响应样点数的整数p，返回一个与由p个频率样本点构成的实频率向量相对应的复频率响应向量h。

freqz命令

freqz函数提供了基于FFT算法所得到的15freqz的命令形式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)用参数选项指定采样频率fsfreqz的命令形式h=freqz(b,a,w)16linspace的命令形式linspace函数的命令形式有两种：linspace（x1,x2） %产生[x1,x2]范围内经线性分割得到的100点的向量linspace（x1,x2，N） %产生[x1,x2]范围内经线性分割得到的N点的向量

例如：w=linspace（0，pi）%

w由（0，pi）上100个等分点组成 w=linspace（0，1000）

%w由（0，1000）上100个等分点组成linspace的命令形式linspace函数的命令形式有两17freqz函数无返回输出参数格式的调用freqz函数还能以无返回输出参数的格式调用： freqz（b，a，256）或freqz（b，a，256，2000）无返回输出参数调用freqz函数的好处是，MATLAB能够自动绘出频率在（0~π）范围内的幅频特性和相频特性图。注意，无返回输出参数调用freqz函数绘出的相频特性不能正确给定在处的值，这是因为在MATLAB中，属于下半个单位圆。

freqz函数无返回输出参数格式的调用freqz函数还能以无18频率响应的实例例：先构成一个截止频率为400Hz的9阶巴特沃思（Butterworth）低通数字滤波器，求出其系数b，a，再求出其256点频率响应。指定的采样频率fs=2000Hz。实现1：先调用butter函数，再调用freqz函数；实现2：无返回输出参数，调用freqz函数；实现3：为得到完整的（0~π）范围内的幅相特性图，可用带返回输出参数的格式调用freqz，在得到各频率点上的频率响应的数值后，使用MATLAB所提供的abs函数和angle函数分别提取幅频特性与相频特性。注意：为了得到连续的相频曲线，需要使用unwrap函数，该函数通过在发生跳变后的各处自动加上，从而除去了相位的跳变，这称为相位的“解卷绕”。

频率响应的实例例：先构成一个截止频率为400Hz的9阶巴特沃19实现1

[b,a]=butter(9,400/1000); %括号中第一个输入是阶数，第二个是归一化截止频率

[h,f]=freqz(b,a,256,2000);或

h=freqz（b，a，256，2000）；实现1 [b,a]=butter(9,400/1000)20实现2实现1中第二条命令的形式可改写为freqz（b，a，256）

或freqz（b，a，256，2000）可自动绘出频率在（0~π）范围内的幅频特性和相频特性图。注意无返回输出参数调用freqz函数绘出的相频特性不能正确给定在ω=π处的值，因为ω=π属于下半个单位圆。实现2实现1中第二条命令的形式可改写为21实现3[b,a]=butter(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π的跳变——卷绕实现3[b,a]=butter(9,400/1000);22解卷绕：unwrap函数为了得到连续的相频曲线，需要使用unwrap函数，该函数通过在发生2π跳变后的各处自动加上2π，从而除去了相位的跳变，这称为相位的“解卷绕”。解卷绕解卷绕：unwrap函数为了得到连续的相频曲线，需要使用u23freqs函数freqs函数用于求取由b和a系数向量定义的模拟滤波器 H（s）的频率响应。使用方法类似freqz函数。与第二章（p32例2-1）采用数组相除方法求取频率响应相比，使用freqs函数要方便很多。freqs函数freqs函数用于求取由b和a系数向量定义的模245.3.2零极点分析zplane函数用于画出线性系统在Z平面上的零极点。有两种使用方法：

1、在已知零极点时，例如某滤波器的零点为-1/2，一对共轭极点为和时，只要输入命令zer=-0.5;pol=0.9*exp(j*2*pi*[-0.30.3]');zplane(zer,pol)

即可画出零极点。（见p70图5-6）5.3.2零极点分析zplane函数用于画出线性系统25用zplane函数绘制零极点图另一种情况：已知系统的系统函数系数向量b和a，则可通过调用zplane(b,a)绘出零极点。这种情形下，zplane函数先求得系统函数的零点和极点，然后绘出零极点图。用zplane函数绘制零极点图另一种情况：已知系统的系统函数265.4 滤波器设计滤波器技术指标，IIR滤波器设计，FIR滤波器设计5.4 滤波器设计滤波器技术指标，275.4.1滤波器的技术指标设计滤波器之前，要定义一组滤波器的技术指标。滤波器的技术指标通常采用幅度指标给出，以保证物理上的可实现性。幅度指标以容差方式给出，如p71图5-7所示。具体指标：通带截止频率Wp；阻带截止频率Ws；通带波纹δp；阻带衰减δs。其中通带波纹和阻带衰减也可以用分贝数给出5.4.1滤波器的技术指标设计滤波器之前，要定义一组滤波28滤波器以容差方式给出的幅度指标

|H(jw)|10通带过渡带阻带pwswsdpd-1w滤波器以容差方式给出的幅度指标|H(jw)|10通带过渡295.4.2IIR滤波器设计MATLAB工具箱提供了几种模拟滤波器原型的产生函数,如巴特沃思(Butterworth),切比雪夫(Chebyshev)滤波器等。还提供了从模拟低通滤波器原型转换为高通、带通和带阻的转换函数，模拟滤波器转换为数字滤波器的双线性变换法和冲激响应不变法，模拟IIR滤波器的阶数选择函数以及数字滤波器直接设计函数等等，使用起来非常方便。本节以巴特沃思滤波器为例

介绍MATLAB工具箱提供的滤波器设计函数使用方法。5.4.2IIR滤波器设计MATLAB工具箱提供了几301．巴特沃思滤波器设计函数

[b,a]

=

butter(n,Wn,options)缺省情况下，返回一个用有理分式表示的低通数字滤波器系统函数的分子分母系数向量b和a。

设计的技术指标只需要指定一个归一化截止频率Wn（这里Wn的单位为：Nyquist

频率

=

1Hz）。参数选项options：选项空缺返回低通滤波器；‘high’则表示返回的是高通数字滤波器；至于带通和带阻滤波器设计，则输入的Wn必须是二个元素的向量；选项options空缺返回带通滤波器，带阻滤波器需加参数options为'stop'。

1．巴特沃思滤波器设计函数[b,a]

=

butter(n31Butter函数举例[B,A]=butter(N,Wn)返回系统函数的长度为N+1的分子系数B和分母系数A。式中Wn为归一化截止频率，0<Wn<1.0。[B,A]=butter(N,Wn,'high')设计一个高通数字滤波器。如果Wn=[W1W2]，且options空缺，则返回一个2N阶的带通滤波器，其通带为：W1<W<W2。[B,A]=butter(N,Wn,'stop')，且Wn=[W1W2]设计2N阶的带阻滤波器。butter(N,Wn,'s')，butter(N,Wn,'high','s')或

butter(N,Wn,'stop','s')调用butter()函数附加一个options为‘s’时，表示设计的是一个模拟滤波器，而指定的截止频率Wn单位此时必须为rad/s。Butter函数举例[B,A]=butter(N,Wn322.巴特沃思滤波器阶次选择函数buttord函数用于给出满足给定频域指标要求的最小阶次滤波器的设计。其调用格式为：

[n,Wn]

=

buttord(Wp,Ws,Rp,Rs) %用于数字滤波器[n,Wn]

=

buttord(Wp,Ws,Rp,Rs,'s')%用于模拟滤波器该函数返回符合技术指标的最小阶数N以及Butterworth滤波器3dB截止频率Wn。设计指标同前：通带起伏不超过Rp，阻带衰减不小于Rs，Wp和Ws分别为归一化通带和阻带截止频率。

2.巴特沃思滤波器阶次选择函数buttord函数用于给出满足33Buttord函数应用：例5-1

例：设计一个数字带通滤波器，通带为1000到2000Hz，止带的起始位置离开通带两边500Hz，也即分别为500Hz和2500Hz，设采样率为10KHz，通带起伏为1dB,止带最小衰减为60dB。解MATLAB实现如下：

[n,Wn]=buttord([10002000]/5000,[5002500]/5000,1,60)得到： n= 12 Wn= 0.19510.4080因此 [b,a]=butter(n,Wn);给出了满足要求的最小阶次巴特沃思带通滤波器。Buttord函数应用：例5-1例：设计一个数字带通滤波器343．经典的IIR滤波器设计方法一种常用的IIR滤波器设计方法是利用已经很成熟的模拟滤波器设计结果。设计步骤如下：（1）按照一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。 （2）根据转换后的技术指标设计模拟低通滤波器的H(s)。 （3）再按一定规则将H(s)转换成H(z)。步骤(3)中最常用的方法:冲激响应不变法（也克称为脉冲响应不变法）和双线性变换法。

3．经典的IIR滤波器设计方法一种常用的IIR滤波器设计方法35冲激响应不变法的原理

冲激响应不变法是通过对模拟滤波器的单位冲激响应ha(nTs),进行采样，将其作为数字滤波器的单位冲激响应h(n)。即：据此设计得到的数字滤波器的频率响应与模拟滤波器的频率响应之间的关系为：因此存在着频率混迭。从而限制了其适用范围。注意：S域频率轴与Z域单位圆之间不存在一一对应关系。冲激响应不变法的原理冲激响应不变法是通过对模拟滤波器的单位36冲激响应不变设计的impinvar()函数impinvar()函数的调用格式是：

[Bz,Az]=impinvar(Bs,As,fs)这一函数把具有[Bs,As]模拟滤波器的转换成了采样频率为fs的数字滤波器的[Bz,Az]。如果没有给定采样频率fs，函数默认为1Hz。

冲激响应不变设计的impinvar()函数impinvar(37双线性变换法原理双线性变换法所定义的s平面与z平面的映射关系为：给定模拟滤波器的系统函数，数字滤波器为：因此S域内的有理分式函数在Z域内仍然保持了有理分式函数形式，从而保证了数字滤波器的物理可实现性。双线性变换法原理双线性变换法所定义的s平面与z平面的映射关系38双线性变换法原理（续）双线性变换将S平面上的频率轴一对一映射到了Z平面的单位圆上，消除了频率混叠问题：但Z平面单位圆和S域频率轴之间存在严重的非线性关系；由此也对其适用范围产生了限制；并且，在滤波器设计中，在频域指标转换时需要进行Pre-warping。

双线性变换法原理（续）双线性变换将S平面上的频率轴一对一映射39双线性变换设计的bilinear函数bilinear()函数的调用格式有三种：

[BzAz]=bilinear(Bs,As,Fs)%把具有[Bs,As]模拟滤波器的转换成数字滤波器的[Bz,Az]。Fs是采样频率。

[ZdPdKd]=bilinear(Z,P,K,Fs)%把模拟滤波器的零点、极点，转换成数字滤波器的零点、极点模型，Fs是采样频率。[AdBdCdDd]=bilinear(A,B,C,D,Fs)%把模拟滤波器的状态方程模型转换成数字滤波器的状态方程模型。Fs是采样频率。注意：这里的采样频率均为实际频率。双线性变换设计的bilinear函数bilinear()函数40例5-3IIR数字滤波器的设计示例

例5-3中直接给出了模拟低通滤波器技术指标（）。如果给定的是数字低通滤波器技术指标，则应先将数字滤波器的技术指标折合成相应的模拟低通滤波器指标，然后采用示例所用的设计方法进行滤波器设计。

例5-3IIR数字滤波器的设计示例例5-3中直接给出了41例5-3设计要求与指标（1）设计一个巴特沃思模拟低通滤波器，指标如下： 通带内 ，波纹小于10dB； 阻带内，衰减不小于20dB。（2）用双线性变换法，将上述模拟低通滤波器变换为数字低通滤波器，采样间隔分别取T=1s，1.5s，比较设计结果。（3）用冲激响应不变法，将上述的模拟低通滤波器变换为数字低通滤波器，采样间隔分别取T=1s，1.5s，比较设计结果。

例5-3设计要求与指标（1）设计一个巴特沃思模拟低通滤波42例5-3的解由于给定的设计指标是模拟滤波器的频域指标，所以，首先选用buttord()函数进行最小阶数的滤波器设计。 [N,Wn]

=

buttord(Wp,Ws,Rp,Rs,'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的解由于给定的设计指标是模拟滤波器的频域指标，所以，43例5-3的解（续）再使用butter()函数，得到模拟低通滤波器的系统函数的系数向量Bs和As： [Bs,As]=butter(N,Wn,'s');根据采样间隔分别为T=1s，1.5s，可得采样频率分别为：Fs1=1Hz，Fs2=2/3Hz，调用bilinear()函数和impinvar()函数，可以得到相应的两种数字滤波器系统函数的系数向量Bz和Az。即： [Bz,Az]=bilinear(Bs，As，Fs) 或 [Bz,Az]=impinvar(Bs，As，Fs)注意，由于低通巴特沃思滤波器并不是带限的，因此对其用冲激响应不变法设计数字滤波器时，将会发生混叠。（参见下张ppt）

例5-3的解（续）再使用butter()函数，得到模拟低通44出现混叠无混叠出现混叠无混叠455.4.3FIR滤波器设计FIR滤波器具有IIR滤波器所没有的线性相位特性，所以在工程中得到了广泛的应用。

MATLAB中除了cremez函数以外，所有用于FIR滤波器设计的函数都是用于设计线性相位的FIR滤波器的。5.4.3FIR滤波器设计FIR滤波器具有IIR滤波器所46FIR滤波器的线性相位条件

FIR数字滤波器的系统函数为：线性相位条件：FIR数字滤波器的h(n)为实数，且满足以下任一条件

h(n)=h(N-1-n)（关于n=(N-1)/2偶对称）h(n)＝-h(N-1-n)（关于n=(N-1)/2奇对称）FIR滤波器的线性相位条件FIR数字滤波器的系统函数为：47根据h(n)的对称情况，分四种类型h(n)奇对称，N为奇数

h(n)奇对称，N为偶数

h(n)偶对称，N为奇数

h(n)偶对称，N为偶数

根据h(n)的对称情况，分四种类型h(n)奇对称，N为奇数48FIR线性相位滤波器的频率响应具有线性相位的FIR滤波器的频率响应函数可以用下式表示：式中或，。注意，式中N为脉冲响应h(n)的长度，是振幅（Amplitude）响应函数，不是幅度（Magnitude）响应函数，可以有正有负，与其相关的相位响应也是一个连续线性函数。对于前述四种h(n),其幅度函数有不同的特点。FIR线性相位滤波器的频率响应具有线性相位的FIR滤波器的频49FIR滤波器设计：窗函数法基本原理：将理想低通滤波器的无限长冲激响应序列截断（等效于加矩形窗）后乘以窗函数，用得到的有限长序列去逼近理想低通滤波器。设所要求的理想数字滤波器的频率响应为, 是与其对应的单位脉冲响应，因此

由于是矩形频率特性，故一定是无限长的非因果序列。而所要设计的是有限长的FIR滤波器，因此用有限长的去逼近无限长的最直接的方法就是截断，即用有限长的窗函数来截取得到：

这就是窗函数设计方法。

FIR滤波器设计：窗函数法基本原理：将理想低通滤波器的无限长50窗函数法设计步骤（1）设所要求的FIR滤波器理想频率响应为：（2）由性能指标确定窗函数W(n)和窗口长度N，由过渡带宽近似于窗口函数主瓣宽求得窗口长度N；（3）由求得实际滤波器的单位脉冲响应h(n)，即求取所设计的FIR滤波器h(n)的系数向量；（4）检验滤波器的性能。窗函数法设计步骤（1）设所要求的FIR滤波器理想频率响应为：51FIR滤波器设计：fir1函数MATLAB中使用窗函数法设计线性相位FIR滤波器的常用设计函数为fir1()，调用格式如下：b=fir1(n,Wn)b=fir1(n,Wn,'ftype')b=fir1(n,Wn,window)b=fir1(n,Wn,'ftype',window)FIR滤波器设计：fir1函数MATLAB中使用窗函数法设计52fir1命令格式的说明fir1命令中n为FIR滤波器的阶数，对于高通、带阻滤波器，n应该取偶数；Wn为滤波器的截止频率，0<Wn<1,对于带通、带阻滤波器，Wn=[W1W2]，且满足W1<W2,对于多带滤波器，Wn=[W1W2W3…]，且满足0<W1<W2<…<1；b为FIR滤波器的系数向量，长度为n+1；‘ftype’是滤波器的类型，默认时为低通或带通滤波器，‘high’为高通滤波器，‘stop’为带阻滤波器；Window为窗函数，是长度为n+1的列向量。fir1命令格式的说明fir1命令中n为FIR滤波器的阶数，53例5-4：低通滤波器例5-4理想数字低通滤波器的截止频率为ω0

（rad），在频率小于ω0处，其幅度为1，在频率为ω0至π处，幅度为0，故它的冲激响应序列h(n)为：其冲激响应为无限长，因此不是因果可实现的。用窗截取其有限长度，再加以因果化，就可以得到一个可实现的线性相位FIR滤波器。例5-4：低通滤波器例5-4理想数字低通滤波器的截止频率54例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个样本，采样频率2Hz

plot(w,abs(H));grid

title('TruncatedSincLowpassFIRFilter');例5-4的解例如取长度51，低通截止频率ω0=0.4π，则55图5-9窗函数法低通滤波器的频率特性

图5-9窗函数法低通滤波器的频率特性56使用其它窗函数可以使用其他形式的窗函数；例如改用长度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)),gridplot(w,20*log10(abs(H))),title('Hamming-WindowedTruncatedSincLPFIRFilter');下一张图将给出其幅频特性，注意与矩形窗情况比较。使用其它窗函数可以使用其他形式的窗函数；57图5-10使用海明窗得到的低通滤波器特性图5-10使用海明窗得到的低通滤波器特性58直接使用FIR滤波器设计函数fir1()

b=fir1(51,0.4);%fir1()中窗函数（window）项空缺代表使用海明窗[H,W]=freqz(b,1,512);plot(W/pi,abs(H));gridtitle(‘Hamming_windowedFIRLPF’)直接使用FIR滤波器设计函数fir1()b=fir1(559fir2()函数函数fir2()也是FIR滤波器设计函数，用于设计具有任意形状频率响应的FIR数字滤波器。例如： n=50;f=[00.40.51];%给定通带、止带范围分别为[00.4]和[0.51]，1为Nyquist频率m=[1100];%给定通带，止带频率响应幅度要求b=fir2(n,f,m) 将返回行向量b为包含n+1个系数的n阶FIR滤波器，其频率响应幅度符合给定的f和m要求。类似于fir1()函数，可以选择窗，命令形式为

b=fir2(n,f,m,window)fir2()函数函数fir2()也是FIR滤波器设计函60FIR滤波器设计：频率采样设计法频率采样设计法是从频域出发，把给定的理想频率响应等间隔采样N个点得到。再对作IDFT，得到然后将h（n）作为所设计的滤波器的单位冲激响应。

FIR滤波器设计：频率采样设计法频率采样设计法是从频域出发，61频率采样设计法的有关说明对进行频率采样得到的H(k)，可以表示成：为保证设计出的FIR滤波器具有线性相位，h(n)应在[0,N-1]上具有对称性。由此进行的理论分析表明，频率采样值的幅度与相位也须满足一定的约束条件，分为4种类型。因此，频率采样法设计FIR滤波器可归结为从满足约束条件的H(k)中求出h(n)问题。

频率采样设计法的有关说明对进行频率采62频率采样设计法的MATLAB实现从满足约束条件的H(k)中求出h(n)，可利用MATLAB中的ifft函数。注意：从时域角度看,求出的h(n)与欲设计的hd(n)是有差别的，两者关系如下：从频域角度看，由频率采样设计所得到的滤波器的频率特性是由H(k)内插形成的，因此对的逼近程度与特性的平滑性有关。

存在时域混叠

频率采样设计法的MATLAB实现从满足约束条件的H(k)中求63例5-5：频率采样法设计举例例5-5所要设计的FIR低通数字滤波器技术指标为：，，例5-5：频率采样法设计举例例5-5所要设计的FIR低通64例5-5的解在给定的图5-12理想低通滤波器特性上选择N=21点进行采样得到|H(k)|，如图5-13。附加相位约束条件对H(k)求反变换IDFT就可得h(n)。例5-5的解在给定的图5-12理想低通滤波器特性上选择N=265图5-14例5-5频率采样法设计低通滤波器(N=21)

最小阻带衰减仅为14dB图5-14例5-5频率采样法设计低通滤波器(N=21)66例5-5的改进设计分析：从图5-14幅度响应曲线可见，N=21时，最小的阻带衰减为14dB，这在绝大多数情况下是不能接受的。但仅仅增大采样点数，比如N=61时，最小的阻带衰减仍为18dB，仍然不能满足通常的止带指标要求。为此，除了增大N以外，还应采用在过渡带插入过渡采样点的办法来除去原始的特性中的突变。

例5-5的改进设计分析：从图5-14幅度响应曲线可见，N=267图5-16例5-5的改进设计

(N=61，增加两个过渡采样点)增加两个过渡采样点后，在N=61时，最小的阻带衰减为45dB，已满足通常的止带衰减要求图5-16例5-5的改进设计

(N=61，增加两个过渡采685.5频谱分析频谱分析的主要作用是对信号的频率构成进行分析，了解其在频带内的分布情况；本课程内仅对确定性信号进行，不涉及随机信号的频谱分析；因此主要工具为FFT，不涉及非线性谱分析或即参量谱分析的内容。

5.5频谱分析频谱分析的主要作用是对信号的频率构成进行69离散傅里叶变换（DFT）在离散信号与系统分析中，DFT有着频谱分析仪的特性。对连续信号与系统通过时域采样，也可以使用DFT进行谱分析。因此，DFT在频谱分析中起了核心的作用。另一方面，从技术上而言，DFT运算存在着高效的快速算法FFT，因此信号的频谱分析中将主要使用MATLAB的fft函数以及ifft函数。

离散傅里叶变换（DFT）在离散信号与系统分析中，DFT有着频70使用DFT对连续信号进行谱分析工程实际应用中，常会遇到连续时间信号，其频谱也是连续函数。为了利用DFT对信号进行频谱分析，须先对进行时域采样，得到，再对进行DFT，得到的则是的傅里叶变换在频域区间上的N点等间隔采样。对带限信号来说，X(k)与在范围内上的采样值除相差一个因子外完全相同。以上说明构成了用DFT对连续信号进行谱分析的理论基础。使用DFT对连续信号进行谱分析工程实际应用中，常会遇到连续时71用DFT进行连续信号频谱分析时注意点

（1）若信号持续时间有限长，则其频谱无限宽。（2）若信号频谱有限，则其信号持续时间必定无限长。（3）实际应用中，所能观察到的信号长度总是有限的。（4）因此，用DFT对连续时间信号进行频谱分析时，采样后产生的频谱混叠失真是不可避的。（5）得到的谱分析结果总是近似的，近似的程度取决于信号带宽、采样率以及信号的截取长度。

用DFT进行连续信号频谱分析时注意点（1）若信号持续时间有72例5-6：连续信号的频谱分析

已知模拟信号对x(t)进行采样，得，分别以T=0.5s，0.1s进行采样，样本数分别取N=20，100，以使被分析信号等长，用DFT计算的并与比较。例5-6：连续信号的频谱分析已知模拟信号73例5-6的分析原信号x(t)的傅立叶变换为其幅度谱为，幅频特性如图5-17所示。由图可知，|X(jΩ)|的-3dB点为Ω=1，而Ω=2π时已衰减到-16dB。因此，如取采样频率Ωs=4π，造成的频率混叠效应会比较小，此时相当于fs=2,T=0.5s。在T=0.1s时，混迭将更小。例5-6的分析原信号x(t)的傅立叶变换为74例5-6解:x(t)采样后所得信号的N点DFT按DFT定义求解题给连续信号x(t)经采样后x(n)的N点DFT实际情况下，信号表达式事先难以确知，因此解析式通常是无法得到的。

例5-6解:x(t)采样后所得信号的N点DFT按DFT定义75两种采样率所得结果的比较两种采样率所得结果的比较76例5-6的结果分析实际应用情况下，通常是截取连续信号的一段样本进行DFT分析。显然，由于xa(t)本身非带限，且又仅截取了[0,T0]的一段，因此，采样后一定有频率混叠（aliasing）存在，再加上截取后的样本值x(n)并不等于前面式中所示的xN(n)，所以，x(n)在数字频率[0,π]内的DFT将不可能与相应模拟频率内的Xa(jΩ)完全相同。造成误差的两个来源中，时域中x(n)与xN(n)的误差在实际应用中无法消除，但由于采样率不够而引入的频率混叠可以通过减小T而得到减轻。比较例5-6在两种采样率下的结果，可以明显看到提高采样率的作用，但其代价为计算量的增加。例5-6的结果分析实际应用情况下，通常是截取连续信号的一段样77其他说明实际工作中，被分析数据是对信号进行截取有限长度得到的，这相当于被分析的是原信号乘上了一个窗函数。这种时域上的截断，在频域上就是所得频谱是原信号的频谱和窗函数频谱的卷积，若信号频率不位于DFT的谱线位置，将会产生频谱的泄漏。在对实际信号进行DFT分析时，泄漏效应是难以避免的，只能通过选择好的窗函数来抑制其影响。用DFT分析信号时，要求被分析的样本序列在所截取的范围内可视为是稳定的。利用DFT进行频谱分析时，DFT的计算结果是样本序列以窗口时间为周期，进行延拓后所得周期波形的傅里叶级数的系数，因此与信号的真实频谱仍有差别。

其他说明实际工作中，被分析数据是对信号进行截取有限长度得到的78fft函数与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()函数。fft函数与ifft函数MATLAB提供fft()函数来计79例5-7

当信号受到噪声污染后，很难从时域中看出它包含的频率成分，这时通常可以用FFT来分析信号频率成分。例5-7中给出了一个仿真产生的受到均匀分布随机噪声干扰的信号，例中用fft()函数进行频谱分析来找出其主要频率成分，并通过filter()函数滤去噪声的干扰。所得结果分别示于图5-20与图5-21中。

例5-7当信号受到噪声污染后，很难从时域中看出它包含的频率805.6 练习课后完成p91~p92练习3-5。5.6 练习课后完成p91~p92练习3-5。81第5章使用MATLAB实现数字信号处理本章主要内容如下：5．1 数字信号处理基本内容及相应的MATLAB工具5．2 信号通过系统的时域分析5．3 信号通过系统的频域和Z域分析5．4 滤波器设计5．5 频谱分析第5章使用MATLAB实现数字信号处理本章主要内容如下：825.1数字信号处理基本内容

及相应的MATLAB工具数字信号处理的基本内容通常分为两部分：离散时间信号与系统分析主要涉及离散时间信号与系统的时域、频域表示，以及信号通过系统的时域、频域分析及其变换域分析。

MATLAB函数库中提供了filter,conv,convmtx,fft,ifft,freqz,impz,zplane等与之相应的函数。数字滤波器设计和谱分析数字滤波器设计包括了无限冲激响应(IIR)和有限冲激响应（FIR）滤波器设计，谱分析又可进一步分为线性谱分析和非线性谱分析。MATLAB为此提供了多种成熟算法的相应函数以及极为丰富的设计工具。5.1数字信号处理基本内容

及相应的MATLAB工具数字信835.2 时域分析卷积，滤波，单位冲激响应5.2 时域分析卷积，滤波，单位冲激响应845.2.1卷积MATLAB提供

conv函数实现标准的一维信号卷积:

例如，若系统h(n)为 >>h=[111]输入序列x(n)为 >>x=[111]则x(n)经过系统h(n)后的MATLAB实现为： >>conv(h,x)或

conv([111],[111])执行后即得到y(n)为 ans= 12321注意：使用conv函数时，h(n)和x(n)都必须是有限长的,，否则不能使用conv函数。5.2.1卷积MATLAB提供conv函数实现标准的一维85例5-1时域离散序列的卷积计算

与图形显示

例5-1(教材p63)

：已知离散信号x(n)和h(n)，求y(n)=x(n)*h(n)，并用图形表示。例5-1时域离散序列的卷积计算

与图形显示例5-1(教86例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);fork=m+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([03001.2]),title(‘x1(n)’) %绘图……（以下省略）例5-1的MATLAB程序Nh=20;Nx=10;m=5; 875.2.2滤波数字滤波器的系统函数H(z)用如下式表示:在MATLAB中，用向量b，a来表示滤波器的系数b(i)和

a(i)。

5.2.2滤波数字滤波器的系统函数H(z)用如下式表示:88滤波器分类当n

=

0，m≠0时，称为AR滤波器，即自回归（AutoRecurrence）滤波器，具无限冲激响应（IIR），也即其单位采样响应h(n)具无限长度；若m

=

0，a(1)≠0，称为MA滤波器，即滑动平均（MovingAverage）滤波器，其单位采样响应h(n)是有限长度，故称有限冲激响应（FIR）滤波器；如果n、m都大于零，称为ARMA滤波器，而其冲激响应也为IIR。滤波器分类当n

=

0，m≠0时，称为AR滤波器，即自回归（89filter函数

MATLAB提供了filter函数来对离散信号进行滤波，表达信号通过系统后的结果。与conv不同的是，filter函数可适用于无限冲激响应系统的情况，但信号仍须是有限长的。例如，一个单极点的低通滤波器系数如下：

>>b=1; %分子系数向量b(i) >>a=[1-0.9]; %分母系数向量a(i)

如果用filter函数实现对信号x滤波，只要调用： >>y=filter(b,a,x);

就可给出输入x经过滤波以后的输出y。filter函数MATLAB提供了filter函数来对离905.2.3单位冲激响应数字滤波器的单位冲激响应定义为输入为单位样本序列时数字滤波器的响应，即：h(n)=T[δ(n)] 其中：5.2.3单位冲激响应数字滤波器的单位冲激响应定义为输入91单位冲激响应的MATLAB实现MATLAB近似实现单位采样信号的方法为：

imp=[1;zeros(p,1)]; %zeros(p,1)产生p个零元素组成的列向量，p是正整数。使用imp后，滤波器的冲激响应可近似得到为：

h=filter(b,a,imp);impz函数可以直接求出数字滤波器的单位冲激响应，即：impz(b,a) 该命令将同时绘出滤波器的单位冲激响应，教材p66图5-2。单位冲激响应的MATLAB实现MATLAB近似实现单位采样信925．3频域和Z域分析

频率响应，零极点分析

5．3频域和Z域分析频率响应，零极点分析935.3.1频率响应MATLAB数字信号处理工具箱有很多函数提供对模拟和数字滤波器的频率响应分析。其中，freqz函数和freqs函数分别返回数字和模拟滤波器的频率响应。工具箱中通常使用的单位频率是Nyquist频率，即采样频率的1/2。注意：就数字滤波器函数来说，其频域指标中的所有频率都以Nyquist频率进行归一化。因此Nyquist频率也称归一化频率。5.3.1频率响应MATLAB数字信号处理工具箱有很多函94关于Nyquist频率的说明例如：系统采样频率为1000Hz，则若数字滤波器的截止频率等于300Hz，经Nyquist频率归一化后，其归一化频率就是300/500=0.6。若将归一化频率转换成数字信号处理教科书中所使用的数字频率(rad)，需乘以π；反之，若乘以采样频率fs的一半，则将归一化频率转换回了模拟域频率(Hz)。归一化频率f应满足0<f<1。对于频率响应而言，归一化频率和模拟频率都可使用。关于Nyquist频率的说明例如：系统采样频率为1000Hz95freqz命令

freqz

函数提供了基于FFT算法所得到的数字滤波器H(z)的频率响应最常用的形式：[h,w]=freqz(b,a,p)

其意义是，对于数字滤波器： freqz函数接受滤波器系数向量b和a，以及一个用于指定所需计算的频率响应样点数的整数p，返回一个与由p个频率样本点构成的实频率向量相对应的复频率响应向量h。

freqz命令

freqz函数提供了基于FFT算法所得到的96freqz的命令形式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)用参数选项指定采样频率fsfreqz的命令形式h=freqz(b,a,w)97linspace的命令形式linspace函数的命令形式有两种：linspace（x1,x2） %产生[x1,x2]范围内经线性分割得到的100点的向量linspace（x1,x2，N） %产生[x1,x2]范围内经线性分割得到的N点的向量

例如：w=linspace（0，pi）%

w由（0，pi）上100个等分点组成 w=linspace（0，1000）

%w由（0，1000）上100个等分点组成linspace的命令形式linspace函数的命令形式有两98freqz函数无返回输出参数格式的调用freqz函数还能以无返回输出参数的格式调用： freqz（b，a，256）或freqz（b，a，256，2000）无返回输出参数调用freqz函数的好处是，MATLAB能够自动绘出频率在（0~π）范围内的幅频特性和相频特性图。注意，无返回输出参数调用freqz函数绘出的相频特性不能正确给定在处的值，这是因为在MATLAB中，属于下半个单位圆。

freqz函数无返回输出参数格式的调用freqz函数还能以无99频率响应的实例例：先构成一个截止频率为400Hz的9阶巴特沃思（Butterworth）低通数字滤波器，求出其系数b，a，再求出其256点频率响应。指定的采样频率fs=2000Hz。实现1：先调用butter函数，再调用freqz函数；实现2：无返回输出参数，调用freqz函数；实现3：为得到完整的（0~π）范围内的幅相特性图，可用带返回输出参数的格式调用freqz，在得到各频率点上的频率响应的数值后，使用MATLAB所提供的abs函数和angle函数分别提取幅频特性与相频特性。注意：为了得到连续的相频曲线，需要使用unwrap函数，该函数通过在发生跳变后的各处自动加上，从而除去了相位的跳变，这称为相位的“解卷绕”。

频率响应的实例例：先构成一个截止频率为400Hz的9阶巴特沃100实现1

[b,a]=butter(9,400/1000); %括号中第一个输入是阶数，第二个是归一化截止频率

[h,f]=freqz(b,a,256,2000);或

h=freqz（b，a，256，2000）；实现1 [b,a]=butter(9,400/1000)101实现2实现1中第二条命令的形式可改写为freqz（b，a，256）

或freqz（b，a，256，2000）可自动绘出频率在（0~π）范围内的幅频特性和相频特性图。注意无返回输出参数调用freqz函数绘出的相频特性不能正确给定在ω=π处的值，因为ω=π属于下半个单位圆。实现2实现1中第二条命令的形式可改写为102实现3[b,a]=butter(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π的跳变——卷绕实现3[b,a]=butter(9,400/1000);103解卷绕：unwrap函数为了得到连续的相频曲线，需要使用unwrap函数，该函数通过在发生2π跳变后的各处自动加上2π，从而除去了相位的跳变，这称为相位的“解卷绕”。解卷绕解卷绕：unwrap函数为了得到连续的相频曲线，需要使用u104freqs函数freqs函数用于求取由b和a系数向量定义的模拟滤波器 H（s）的频率响应。使用方法类似freqz函数。与第二章（p32例2-1）采用数组相除方法求取频率响应相比，使用freqs函数要方便很多。freqs函数freqs函数用于求取由b和a系数向量定义的模1055.3.2零极点分析zplane函数用于画出线性系统在Z平面上的零极点。有两种使用方法：

1、在已知零极点时，例如某滤波器的零点为-1/2，一对共轭极点为和时，只要输入命令zer=-0.5;pol=0.9*exp(j*2*pi*[-0.30.3]');zplane(zer,pol)

即可画出零极点。（见p70图5-6）5.3.2零极点分析zplane函数用于画出线性系统106用zplane函数绘制零极点图另一种情况：已知系统的系统函数系数向量b和a，则可通过调用zplane(b,a)绘出零极点。这种情形下，zplane函数先求得系统函数的零点和极点，然后绘出零极点图。用zplane函数绘制零极点图另一种情况：已知系统的系统函数1075.4 滤波器设计滤波器技术指标，IIR滤波器设计，FIR滤波器设计5.4 滤波器设计滤波器技术指标，1085.4.1滤波器的技术指标设计滤波器之前，要定义一组滤波器的技术指标。滤波器的技术指标通常采用幅度指标给出，以保证物理上的可实现性。幅度指标以容差方式给出，如p71图5-7所示。具体指标：通带截止频率Wp；阻带截止频率Ws；通带波纹δp；阻带衰减δs。其中通带波纹和阻带衰减也可以用分贝数给出5.4.1滤波器的技术指标设计滤波器之前，要定义一组滤波109滤波器以容差方式给出的幅度指标

|H(jw)|10通带过渡带阻带pwswsdpd-1w滤波器以容差方式给出的幅度指标|H(jw)|10通带过渡1105.4.2IIR滤波器设计MATLAB工具箱提供了几种模拟滤波器原型的产生函数,如巴特沃思(Butterworth),切比雪夫(Chebyshev)滤波器等。还提供了从模拟低通滤波器原型转换为高通、带通和带阻的转换函数，模拟滤波器转换为数字滤波器的双线性变换法和冲激响应不变法，模拟IIR滤波器的阶数选择函数以及数字滤波器直接设计函数等等，使用起来非常方便。本节以巴特沃思滤波器为例

介绍MATLAB工具箱提供的滤波器设计函数使用方法。5.4.2IIR滤波器设计MATLAB工具箱提供了几1111．巴特沃思滤波器设计函数

[b,a]

=

butter(n,Wn,options)缺省情况下，返回一个用有理分式表示的低通数字滤波器系统函数的分子分母系数向量b和a。

设计的技术指标只需要指定一个归一化截止频率Wn（这里Wn的单位为：Nyquist

频率

=

1Hz）。参数选项options：选项空缺返回低通滤波器；‘high’则表示返回的是高通数字滤波器；至于带通和带阻滤波器设计，则输入的Wn必须是二个元素的向量；选项options空缺返回带通滤波器，带阻滤波器需加参数options为'stop'。

1．巴特沃思滤波器设计函数[b,a]

=

butter(n112Butter函数举例[B,A]=butter(N,Wn)返回系统函数的长度为N+1的分子系数B和分母系数A。式中Wn为归一化截止频率，0<Wn<1.0。[B,A]=butter(N,Wn,'high')设计一个高通数字滤波器。如果Wn=[W1W2]，且options空缺，则返回一个2N阶的带通滤波器，其通带为：W1<W<W2。[B,A]=butter(N,Wn,'stop')，且Wn=[W1W2]设计2N阶的带阻滤波器。butter(N,Wn,'s')，butter(N,Wn,'high','s')或

butter(N,Wn,'stop','s')调用butter()函数附加一个options为‘s’时，表示设计的是一个模拟滤波器，而指定的截止频率Wn单位此时必须为rad/s。Butter函数举例[B,A]=butter(N,Wn1132.巴特沃思滤波器阶次选择函数buttord函数用于给出满足给定频域指标要求的最小阶次滤波器的设计。其调用格式为：

[n,Wn]

=

buttord(Wp,Ws,Rp,Rs) %用于数字滤波器[n,Wn]

=

buttord(Wp,Ws,Rp,Rs,'s')%用于模拟滤波器该函数返回符合技术指标的最小阶数N以及Butterworth滤波器3dB截止频率Wn。设计指标同前：通带起伏不超过Rp，阻带衰减不小于Rs，Wp和Ws分别为归一化通带和阻带截止频率。

2.巴特沃思滤波器阶次选择函数buttord函数用于给出满足114Buttord函数应用：例5-1

例：设计一个数字带通滤波器，通带为1000到2000Hz，止带的起始位置离开通带两边500Hz，也即分别为500Hz和2500Hz，设采样率为10KHz，通带起伏为1dB,止带最小衰减为60dB。解MATLAB实现如下：

[n,Wn]=buttord([10002000]/5000,[5002500]/5000,1,60)得到： n= 12 Wn= 0.19510.4080因此 [b,a]=butter(n,Wn);给出了满足要求的最小阶次巴特沃思带通滤波器。Buttord函数应用：例5-1例：设计一个数字带通滤波器1153．经典的IIR滤波器设计方法一种常用的IIR滤波器设计方法是利用已经很成熟的模拟滤波器设计结果。设计步骤如下：（1）按照一定规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。 （2）根据转换后的技术指标设计模拟低通滤波器的H(s)。 （3）再按一定规则将H(s)转换成H(z)。步骤(3)中最常用的方法:冲激响应不变法（也克称为脉冲响应不变法）和双线性变换法。

3．经典的IIR滤波器设计方法一种常用的IIR滤波器设计方法116冲激响应不变法的原理

冲激响应不变法是通过对模拟滤波器的单位冲激响应ha(nTs),进行采样，将其作为数字滤波器的单位冲激响应h(n)。即：据此设计得到的数字滤波器的频率响应与模拟滤波器的频率响应之间的关系为：因此存在着频率混迭。从而限制了其适用范围。注意：S域频率轴与Z域单位圆之间不存在一一对应关系。冲激响应不变法的原理冲激响应不变法是通过对模拟滤波器的单位117冲激响应不变设计的impinvar()函数impinvar()函数的调用格式是：

[Bz,Az]=impinvar(Bs,As,fs)这一函数把具有[Bs,As]模拟滤波器的转换成了采样频率为fs的数字滤波器的[Bz,Az]。如果没有给定采样频率fs，函数默认为1Hz。

冲激响应不变设计的impinvar()函数impinvar(118双线性变换法原理双线性变换法所定义的s平面与z平面的映射关系为：给定模拟滤波器的系统函数，数字滤波器为：因此S域内的有理分式函数在Z域内仍然保持了有理分式函数形式，从而保证了数字滤波器的物理可实现性。