离散傅里叶变换DFTPPT课件

1、1二、DFT和Z变换的关系 设序列x(n)的长度为N， 其Z变换和DFT分别为：1010( ) ( )( )( ) ( )( )0kN-1NnnNknNnX zZT x nx n zX kDFT x nx n W 比较上面二式可得关系式22( )( ),0-1(3.1.3)( )(),0-1(3.1.4)jkNz ejkNX kX zkNX kX ekN （3.1.3）式表明序列x(n)的N点DFT是x(n)的Z变换在单位圆上的N点等间隔采样； （3.1.4）式则说明X(k)为x(n)的傅里叶变换X(ej)在区间0，2 上的N点等间隔采样。第1页/共68页2 实际上， 任何周期为N的周期序列

2、都可以看作长度为N的有限长序列x(n)的周期延拓序列， 即( )()(3.1.5)( )( )( )(3.1.6)mNx nx nmNx nx nRn 为了以后叙述方便， 将(3.1.5)式用如下形式表示： ( )( )Nx nx n 三、DFT的隐含周期性 (),kk mNNNWWk m N 均为整数 ( )x n X(k)X(k)(k)(3.1.10)NR 有限长序列x(n)的离散傅里叶变换X(k)，正好是x(n)的周期延拓序列x(n)N的离散傅里叶级数系数 的主值序列，即X(k) 第2页/共68页3图 3.1.2 有限长序列及其周期延拓 第3页/共68页47.3 离散傅里叶变换（DFT）

3、 定义DFT：用类似于例7.9中的方法，可把（7.3）式写成矩阵乘法运算。Xk=xn*Wnk其中，xn为序列行向量，Wnk是一 NN 阶方阵， 而 称为旋转因子。Wnk 的MATLAB表示： Wnk = WN.(0:N-1* 0:N-1)10( )DFT ( )( ) , 0, , 1NknNnX kx nx n WkN nk0 00 10 (1)1 01 11 (1)(1) 0(1) 1(1) (1)NNNNNNNNNNNNNNNWWWWWWWWWW 2jNNWe 第4页/共68页5用矩阵乘法计算N点DFT的程序如下: MATLAB程序q73a.m%用矩阵乘法计算N点DFTclear;clo

4、se allxn=input(请输入序列x= );N = length(xn); %n=0:N-1; k=n; nk=n*k;% 生成NN方阵WN=exp(-j*2*pi/N); % 旋转因子Wnk=WN.nk; % 产生旋转因子矩阵Xk=xn*Wnk; % 计算N点DFT 这种方法计算DFT概念清楚、编程简单，但占用内存大、运行速度低，所以不实用。MATLAB基础部分提供了fft、ifft、fft2、ifft2等等快速计算傅里叶变换的函数，使DFT的运算速度量提高了若干数量级，在后面的例题中均直接调用这些函数。第5页/共68页6例7.15 序列的离散傅立叶变换 求复正弦序列 余弦序列 正弦序

5、列的离散傅立叶变换，分别按N =16和N =8进行计算。绘出幅频特性曲线， 进行比较讨论。 解：直接产生序列x1n、x2n和x3n，调用fft函数求解 j81( )e( )nNx nRn 2( )cos( )8NxnnRn 3( )sin( )8Nx nnRn 第6页/共68页7%第七章例7.15程序q715 % DFT计算clear;close allN=16;N1=8;%产生序列x1(n),计算DFTx1(n)n=0:N-1;x1n=exp(j*pi*n/8); %产生x1(n)X1k=fft(x1n,N); %计算N点DFTx1(n)Xk1=fft(x1n,N1); %计算N1点DFTx

6、1(n)%产生序列x2(n),计算DFTx2(n)x2n=cos(pi*n/8);X2k=fft(x2n,N); %计算N点DFTx2(n)Xk2=fft(x2n,N1); %计算N1点DFTx1(n)%产生序列x3(n),计算DFTx3(n)x3n=sin(pi*n/8);X3k=fft(x3n,N); %计算N点DFTx3(n)Xk3=fft(x3n,N1); %计算N1点DFTx1(n)第7页/共68页8%绘图subplot(2,3,1);stem(n,abs(X1k),.);title(16点 DFTx1(n);xlabel(k);ylabel(|X1(k)|)subplot(2,3,

7、2);stem(n,abs(X2k),.);title(16点 DFTx2(n);xlabel(k);ylabel(|X2(k)|)subplot(2,3,3);stem(n,abs(X3k),.);title(16点 DFTx3(n);xlabel(k);ylabel(|X3(k)|)k=0:N1-1;subplot(2,3,4);stem(k,abs(Xk1),.);title(8点 DFTx1(n);xlabel(k);ylabel(|X1(k)|)subplot(2,3,5);stem(k,abs(Xk2),.);title(8点 DFTx2(n);xlabel(k);ylabel(|

8、X2(k)|)subplot(2,3,6);stem(k,abs(Xk3),.);title(8点 DFTx3(n);xlabel(k);ylabel(|X3(k)|) 第8页/共68页9 在截取16点时，得到的是完整的余弦波形；而截取8点时，得到的是半截的余弦波形，当然有大量的谐波成分。第9页/共68页10例7.16 验证N点DFT的物理意义（1）绘出幅频曲线和相频曲线。（2）计算并图示x(n)的8点DFT。（3）计算并图示x(n)的16点DFT。解: 序列x(n)的点DFT的物理意义是 在0，2 上进行点等间隔采样。程序先密集采样，绘制出幅频曲线图。然后再分别做8点和16点DFT来验证这个

9、采样关系。 程序略。441( )( ),() ( )1jjjex nR nX eFT x ne 求求得得，()jX e 第10页/共68页11图7.16 离散傅里叶变换与傅里叶变换的采样关系第11页/共68页12频域采样定理 频域采样是指对有限长序列的傅氏变换在频域离散化得到X(k)的过程。 本节讨论两个基本问题：a：频域采样（DFT）不失真的条件，即由 X (k)不失真地恢复x(n)的条件b：用X (k)表示X (z)和 的插值公式(内插公式)()jX e 第12页/共68页13图3-3-1 时域恢复示意图 结论：若序列长度为L，频域采样点数（DFT的长度）为N，且LN，则频域采样后可不失真

10、地恢复原序列 ；但若LN，则频域采样后不能不失真地恢复原序列 。( )x n( )x n第13页/共68页14例7.17 频域与时域采样对偶性 （1）产生三角波序列 （2）对M = 40，计算x(n)的64点DFT，并图示x(n)和 X(k) = DFTx(n)，k = 0, 1, , 63。 （3）对（2）中所得X(k)在 0，2 上进行32点抽样得 （4）求 的32点IDFT，即 （5）绘出 的波形图，评述它与x(n)的关系。0 / 2( ) / 2 nnMx nMnMnM 1( )(2 ) 0 131XkXkk ， ， ，1( )X k11( )IDFT( )n = 0 1x nXk ，

11、 ，3 31 1132( )xn第14页/共68页15 % 第七章例7.17程序q717 % 时域与频域采样的对偶性验证 clear;close all M=40;N=64;n=0:M; %产生M长三角波序列x(n) xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=xa,xb; Xk=fft(xn,64);%64点FFTx(n) X1k=Xk(1:2:N);%隔点抽取Xk得到X1(K) x1n=ifft(X1k,N/2);%32点IFFTX1(k)得到x1(n) nc=0:3*N/2;% 取97点为观察区 xc=x1n(mod(nc,N/2)+1);%x1(

12、n)的周期延拓序列 subplot(3,2,1);stem(n,xn,.) title(40点三角波序列x(n);xlabel(n);ylabel(x(n)第15页/共68页16 k=0:N-1; subplot(3,2,3);stem(k,abs(Xk),.) title(64点DFTx(n);xlabel(k);ylabel(|X(k)|) k=0:N/2-1; subplot(3,2,4);stem(k,abs(X1k),.) title(X(k)的隔点抽取);xlabel(k);ylabel(|X1(k)|) n1=0:N/2-1; subplot(3,2,2);stem(n1,x1n

13、,.) title(32点IDFTX1(k);xlabel(n);ylabel(x1(n) subplot(3,1,3);stem(nc,xc,.) title(x1(n)的周期延拓序列);xlabel(n);ylabel(x(mod(n,32) set(gcf,color,w)第16页/共68页17第17页/共68页18 由于频域在0，2 上的采样点数N（N = 32）小于x(n)的长度M（M = 40），所以，产生时域混叠现象，不能由X1(k)恢复原序列x(n)。只有满足NM时，可由频域采样X1(k)得到原序列x(n)。 这就是频域采样定理。对NM的情况，请读者自己编程验证。1( )IDF

14、T( )x nX k 第18页/共68页19用DFT计算线性卷积112120( )( )( )()()( )LLLmy nx nx nx m xnmRn 1122( )( )( )( )X kDFT x nXkDFT x n 0kL-1则由时域循环卷积定理有 Y(k)=DFTy(n)=X1(k)X2(k), 0kL-1如果第19页/共68页20 由此可见， 循环卷积既可在时域直接计算， 也可以按照图3.4.1所示的计算框图， 在频域计算。 由于DFT有快速算法FFT， 当N很大时， 在频域计算的速度快得多， 因而常用DFT(FFT)计算循环卷积。 图 3.4.1 用DFT计算循环卷积 第20页

15、/共68页21 在实际应用中，需要计算两个序列的线性卷积， 与计算循环卷积一样， 为了提高运算速度， 也希望用DFT(FFT)计算线性卷积。 而DFT只能直接用来计算循环卷积， 为此导出线卷积和循环卷积之间的关系以及循环卷积与线性卷积相等的条件。 假设h(n)和x(n)都是有限长序列， 长度分别是N和M。 它们的线性卷积和循环卷积分别表示如下： 1010( )( )( )() ()( )( )( )() ()( )NlmLcLLmy nh nx nh m x nmy nh nx nh m x nmRn (3.4.1) (3.4.2) 第21页/共68页22 其中，LmaxN, M， 所以101

16、0( )()()( )() ()( )NcLmqNLqmy nh mx nmqL Rnh m x nmqL Rn ( )(),Lqx nx nqL 对照式(3.4.1)可以看出， 上式中 ( )()( )clLqy ny nqL Rn (3.4.3) 10() ()()Nlmh m x nqLMy nqL 即(3.4.3)式说明，yc(n)等于yl(n)以L为周期的周期延拓序列的主值序列。 yl(n)长度为N+M-1，只有当循环卷积长度L N+M-1时，yl(n)以L为周期进行周期延拓才无混叠现象。第22页/共68页23图 3.4.2 线性卷积与循环卷积 0123451234h(n) x(n)

17、nL 60123451234nL 867h(n) x(n)0123451234nL 1067h(n) x(n)（ d ）（ e ）( f )0123451234nN M1 867h(n) x(n)*nM 5012341x(n)nN 401231h(n)（ a ）（ b ）（ c ）89* * 189 10第23页/共68页24 利用DFT进行线性卷积的步骤如下： 1.将序列x(n)和h(n)补零延长，使其长度 2.做x(n)和h(n)的长为L点的DFT得到X(k)和H(k)，求它们的积Y(k)=X(k)H(k) 3.求Y(k)的IDFT并取前 N N1 1 点获得线性卷积的结果为LN1=N+M

18、-1 y(n)=IDFTY(k)，0nN1 第24页/共68页25图 3.4.3 用DFT计算线性卷积框图 补L N个零点L点DFT补L M个零点L点DFTL点IDFTy(n)h(n)x(n)第25页/共68页26例7.18 快速卷积 快速卷积就是根据快速卷积就是根据DFT的循环卷积性质，将的循环卷积性质，将时域卷积转换为频域相乘，最后再进行时域卷积转换为频域相乘，最后再进行IDFT得到得到时域卷积序列时域卷积序列y(n)。其中时域和频域之间的变换。其中时域和频域之间的变换均用均用FFT实现，所以使卷积速度大大提高。框图实现，所以使卷积速度大大提高。框图如下如下: L点FETL点FETL点FE

19、TL点FFTL点FFTL 点IFFT第26页/共68页27 % 第七章例7.18程序q718 % 快速卷积计算 clear;close all xn=input(请输入x(n)序列：xn= 书上用 sin(0.4*1:15); hn=input(请输入h(n)长度：hn= 书上用 0.9.(1:20); M=length(xn); N=length(hn); nx=1:M; nh=1:N; %循环卷积等于线性卷积的条件：循环卷积区间长度L=M+N-1 L=pow2(nextpow2(M+N-1);%取L为大于等于且最接近（N+M-1）的2的正次幂 tic, %快速卷积计时开始第27页/共68页

20、28 Xk=fft(xn,L);%L点FFTx(n) Hk=fft(hn,L);%L点FFTh(n) Yk=Xk.*Hk; %频域相乘得Y(k) yn=ifft(Yk,L);%L点IFFT得到卷积结果y(n) toc %快速卷积计时结束 subplot(2,2,1),stem(nx,xn,.); ylabel(x(n) subplot(2,2,2),stem(nh,hn,.); ylabel(h(n) subplot(2,1,2);ny=1:L; stem(ny,real(yn),.);ylabel(y(n) tic, yn=conv(xn,hn);%直接调用函数conv计算卷积与快速卷积比较

21、 toc第28页/共68页29第29页/共68页30 4. 用DFT进行谱分析的误差问题 DFT(实际中用FFT计算)可用来对连续信号和数字信号进行谱分析。 (1) 混叠现象。 (2) 栅栏效应。 (3) 截断效应。 根据傅里叶变换的频域卷积定理有 ()1() ( )()()21()()2jjjNjjNY eFT y nX eReX eRed 第30页/共68页31 幅度谱RN()曲线如图3.4.11所示(RN（）以2为周期， 只画低频部分)。 图中，|2/N的部分称为主瓣， 其余部分称为旁瓣。 例如， x(n)=cos(0n), 0=/4其频谱为 其中 1()2() ( )sin(/ 2)(

22、)( )( )sin(/ 2)jNjjjNNNX eFT x nNReFT RneRe () (2)(2)44jlX ell 第31页/共68页32 图 3.4.11 矩形窗函数的幅度谱 第32页/共68页33图 3.4.12 加矩形窗前后的频谱 cos()4n 泄露谱间干扰第33页/共68页34例7.19 用DFT求连续信号频谱 在计算机上用DFT对模拟信号进行谱分析时，只能以有限大的采样频率fs对模拟信号采样有限点样本序列（等价于截取模拟信号一段进行采样）作DFT变换，得到模拟信号的近似频谱。其误差主要来自以下因素： 截断效应（频谱泄露和谱间干扰） 频谱混叠失真 因素使谱分辨率（能分辨开的

23、两根谱线间的最小间距）降低，并产生谱间干扰； 因素使折叠频率（ fs / 2）附近的频谱产生较大失真。 第34页/共68页35例7.19 用DFT求连续信号频谱 加大截取长度Tp可提高频率分辨率；选择合适的窗函数可降低谱间干扰；而频谱混叠失真要通过提高采样频率fs和（或）预滤波（在采样之前滤除折叠频率以外的频率成分）来改善。 编写程序q719.m验证截断效应及加窗的改善作用，先选取以下参数： 采样频率fs = 400Hz，T = 1/ fs 采样信号序列 对x(n)作4096点DFT作为xa(t)的近似频谱Xa(jf)。a( )() ( ), n=0,1,N-1x nx nT w n 第35页

24、/共68页36 % 第七章例7.19程序q719 % 用DFT作谱分析 clear;close all fs=400; T=1/fs;%采样频率为400Hz Tp=0.04; N=Tp*fs;%采样点数N N1=N, 4*N, 8*N;% 设定三种截取长度供调用 st=|X1(jf)|;|X4(jf)|;|X8(jf)|;% 设定三种标注语句供调用 %矩形窗截断 for m=1:3 n=1:N1(m); xn=cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T);%产生采样序列x(n) Xk=fft(xn,4096); %4096点DFT，用FFT实现

25、fk=0:4095/4096/T; subplot(3,2,2*m-1) plot(fk,abs(Xk)/max(abs(Xk);ylabel(st(m,:) if m=1 title(矩形窗截取);end end第36页/共68页37 %加hamming窗改善谱间干扰 for m=1:3 n=1:N1(m); wn=hamming(N1(m); %调用工具箱函数hamming产生N长hamming窗序列wn xn=(cos(200*pi*n*T)+sin(100*pi*n*T)+cos(50*pi*n*T).*wn; Xk=fft(xn,4096); %4096点DFT，用FFT实现 fk=

26、0:4095/4096/T; subplot(3,2,2*m) plot(fk,abs(Xk)/max(abs(Xk);ylabel(st(m,:) if m=1 title(Hamming窗截取);end end第37页/共68页38图7.19第38页/共68页39例7.19 用DFT求连续信号频谱 如图7.19所示。图中X1(jf)，X4(jf)和X8(jf)分别表示Tp = 0.04s，0.04*4s和0.04*8s时的谱分析结果。由图可见，由于截断使原频谱中的单频谱线展宽（又称之为泄漏）， Tp 越大泄漏越小，频率分辨率越高。 Tp = 0.04s时，25Hz与50Hz两根谱线已分辨不

27、清了。所以实际谱分析的截取时间Tp是由频率分辨率决定的。 另外，在本应为零的频段上出现了一些参差不齐的小谱包（称为谱间干扰）。谱间干扰的大小取决于加窗的类型。用矩形窗比用Hamming窗的频率分辨率高（泄漏小），但谱间干扰刚好相反。 第39页/共68页407.5 FIR数字滤波器设计 滤波器的特性指标 用绝对值1,21,2表示; ; 用分贝Rp,RsRp,Rs表示1-11p1012s101120log0,120log01RR 第40页/共68页41（1）窗函数法设计FIR滤波器 先根据 c和N求出相应的理想滤波器单位脉冲响应hd(n)。 第二步要选择合适的窗函数w(n)来截取hd(n) 的适当

28、长度（即阶数），以保证实现要求的阻带衰减；最后得到FIR滤波器单位脉冲响应h(n)=hd(n).*w(n),即其系数。 )()(sin(de )e (21)(cjjdccnnHnhnd第41页/共68页42（2）等波纹最佳一致逼近法（2）等波纹最佳一致逼近法: 信号处理工具箱采用remez算法实现线性相位FIR数字滤波器的等波纹最佳一致逼近设计。其优点是，设计指标相同时，使滤波器阶数最低；或阶数相同时，使通带最平坦，阻带最小衰减最大；通带和阻带均为等波纹形式，最适合设计片段常数特性的滤波器。其调用格式如下： b = remez(N, f, m, w, ftype)其中N由remezord函数求

29、出: N, fo, mo, w = remezord(f, m, dev, Fs)输入变元dev为各逼近频段允许的波纹振幅。remez函数可直接调用remezord返回的参数如下： b=remez(N, fo, mo, w) 第42页/共68页43例7.25 窗函数法设计数字滤波器 分别用矩形窗和Hamming窗设计线性相位FIR低通滤波器。要求通带截止频率 c = /4，单位脉冲响应h(n)的长度N = 21。绘出h(n)及其幅频响应特性曲线。 先求出相应的理想滤波器（本例应为理想低通）单位脉冲响应hd(n)，再根据阻带最小衰减选择合适的窗函数w(n)，最后得到FIR滤波器单位脉冲响应h(n

30、)=hd(n).*w(n)。 第43页/共68页44例7.25 窗函数法设计数字滤波器 本题中， c = /4，N = 21，所以线性相位理想低通滤波器的单位脉冲响应为: 为了满足线性相位FIR滤波器条件h(n) = h(N-1-n)，要求 = (N-1)/2 = 10。 信号处理工具箱中有窗生成函数boxcar，hamming，hanning和blackman等。 )(4)(sin)(dnnnh第44页/共68页45例7.25 窗函数法设计数字滤波器 对两种窗函数的设计结果分别如右图7.25-1和图7.25-2所示。 h(n)h(n)nh(n)h(n)n图7.25-1图7.25-2第45页/

31、共68页46工具箱设计函数fir1和fir2 MATLAB提供了基于窗函数法的FIR滤波器设计函数fir1和fir2,其功能及用法如下。 fir1功能：标准频率响应形状。格式：b=fir1(N, wc, ftype, window)。当wc=wc1，wc2时，得到的是带通滤波器。 当ftype=high时，设计高通FIR滤波器； 当ftype=stop时，设计带阻FIR滤波器。 fir2功能：任意频率响应形状。格式：b = fir2(N, f, m, window) 第46页/共68页47例7.26 窗函数法设计带通滤波器 使用fir1函数b = fir1(N, wc, window) 编程参

32、数 c为行向量 c = lp/ ， hp/ 根据阻带最小衰减Rs = 60dB选择窗函数类型和阶次。可以查上面列出的“窗函数设计滤波器时的阶数选择表”。选blackman窗，其滤波器阻带最小衰减可达到74dB，其窗口长度M由过渡带宽度B = 0.15 决定，Blackman窗设计的滤波器过渡带宽度为12 /M，故M取80。因M = N+1，所以滤波器阶数N = 79。 第47页/共68页48例7.27 用remez函数设计低通滤波器 解：先由题意计算设计参数 f = 1/4,5/16，m = 1,0; dev的计算稍复杂一些，由于 所以 有了这几个参数就可以调用remezord和remez函数

33、了.)2(devlog(20 , ) 1 (dev1) 1 (dev120lospAgRpps/20/20(/20)dev(1)(101)/(101) , dev(2)10RRA第48页/共68页49例7.27 用remez函数设计低通滤波器 横线为-3dB，两条竖线分别位于频率 /4和5 /16。显然，通带指标稍有富裕，过渡带宽度和阻带最小衰减刚好满足指标要求。 程序输出的幅频特性第49页/共68页50例7.28 remez函数设计高通滤波器 观察等波纹逼近法中加权系数w( )及滤波器阶数N的作用和影响。期望逼近的滤波器通带为3 /4， ，阻带为0，23 /32。 解:在滤波器设计中，技术指

34、标越高，实现滤波器的阶数也就越高。另外，对固定的阶数，通带与阻带指标可以互换，过渡带宽度与通带波纹和阻带衰减指标可以互换。 取f = 0, 3/4, 23/32, 1，m = 0, 0, 1, 1。其余参数分三种情况进行设计，N = 30，w = 1, 1；N = 30，w = 1, 5; N = 60，w = 1, 1。 第50页/共68页51例7.28 remez函数设计高通滤波器 程序运行结果如图 由图可见，w较大的频段逼近精度较高；w较小的频段逼近精度较低。N较大时逼近精度较高,N较小时逼近精度较低 。第51页/共68页527.6 IIR数字滤波器设计 IIR数字滤波器设计的主要方法是

35、先设计低通模拟滤波器，进行频率变换，将其转换为相应的（高通、带通等）模拟滤波器，再转换为高通、带通或带阻数字滤波器。对设计的全过程的各个步骤，MATLAB都提供了相应的工具箱函数，使IIR数字滤波器设计变得非常简单。本节主要结合例题介绍这些IIR滤波器设计的工具箱函数。 IIR数字滤波器的设计步骤由以下的流程图来表示。下面以巴特沃斯滤波器设计函数为典型，介绍此流程图中函数的功能和用法。第52页/共68页53IIR数字滤波器设计流程图模拟低通滤波器原型设计Buttap,cheb1ap,cheb2apbesselap，ellipap函数频率变换(变为高通，带通，带阻等)lp2lp，lp2hp，lp

36、2bp，lp2bs模拟数字变换bilinearimpinvar合为一步的设计函数 butter，cheb1，cheb2，ellip， besself求最小阶数Nbuttord, cheb1ordcheb2ord,ellipord设计指标wp, ws, Rp, RsN，wc滤波器系数B，A滤波器系数B，A第53页/共68页54巴特沃斯滤波器设计流程 （1）求最小阶数N的函数buttordN, wc = buttord (wp, ws, Rp, Rs, s) 根据滤波器指标wp，ws，Rp，Rs，求出巴特沃斯模拟滤波器的阶数N及频率参数wc，此处wp，ws及wc均以弧度/秒为单位。 （2）得到N后

37、,调用设计函数buttapz，p，k = buttap(N) 得到z, p, k后，很容易求出滤波器系数B，A。 （3）调用模拟频率变换函数lp2lpBt, At = lp2lp(B, A, wo) （4）调用模拟数字变换函数 Bd, Ad = bilinear (B, A, Fs) 第54页/共68页55集成的数字滤波器设计函数 把（2）、（3）、（4）合为一步的数字滤波器设计函数butter(N, wc, ftype) B, A = butter (N, wc) 设计低通或带通数字滤波器系数B，A（当为带通滤波器时，第(1)类函数由wp = wp1, wp2会自动生成wc = w1, w2

38、）。 B, A = butter (N, wc, high) 设计高通数字滤波器系数B，A。 B, A = butter (N, wc, stop) 设计带阻数字滤波器系数B，A。 butter(N, wc, ftype)还有零极增益和状态空间形式，读者可用help命令查阅。第55页/共68页56例7.29 巴特沃斯模拟滤波器设计 设计一个低通巴特沃斯模拟滤波器，指标如下。通带频率：fp = 3400Hz，最大衰减：Rp = 3dB阻带频率：fs = 4000Hz，最小衰减：As = 40dB 解:它的系统函数完全由阶数N和3dB截止频率 c决定。而N和 c是由滤波器设计指标决定的。 取 c

39、= c1，通带指标刚好，阻带指标富裕； 取 c = c2，则阻带指标刚好，通带指标富裕。 MATLAB工具箱函数buttord，butter就是根据以上公式编写的。因此就无需再记忆这些公式了。 1)(10 , 1)10(21s0.1sc2210.1pc1pNANR-第56页/共68页57模拟转换为数字：脉冲响应不变法 模拟滤波器离散化的基本方法有脉冲响应不变法和双线性变换法。 脉冲响应不变法及impinvar函数单极点的N阶模拟滤波器Ha(s)，用部分分式展开为脉冲响应不变法的数字化结果为 工具箱函数impinvar可实现以上计算，格式为 Bz, Az = impinvar(B, A, Fs)

40、 NkkkssAsH1a)()(NkTskzAzHk11)e1 ()(第57页/共68页58模拟转换为数字：双线性变换法 双线性变换法函数bilinear 双线性变换法的原理是用 代换Ha(s)中的s值，得到H(z)。bilinear函数用来实现这个转换。其使用格式为 Bz, Az = bilinear(B, A, Fs) 脉冲响应不变法的缺点是存在频率混叠失真。双线性变换法可完全消除频率混叠失真，缺点是存在非线性频率失真。 11112zzTs第58页/共68页59例7.30 模拟低通转换为数字低通 已知一模拟滤波器的系统函数为 分别用脉冲响应不变法和双线性变换法将Ha(s)转换成数字滤波器系

41、统函数H(z)，并图示Ha(s)和H(z)的幅频响应曲线。 程序中的核心语句是以下两条： d,c=impinvar(b,a,Fs)% 用impinvar 函数离散化 f,e=bilinear(b,a,Fs)% 用bilinear 函数离散化10001000)(assH第59页/共68页60例7.30 模拟低通转换为数字低通 图形结果如图7.30所示。由图7.30（b）可见，对脉冲响应不变法，采样频率Fs越高（T越小），混叠越小；由图7.30（c）可见，对双线性变换法，无频率混叠，但存在非线性失真。第60页/共68页61例7.31 切比雪夫数字滤波器设计 解：切比雪夫型滤波器通带内为等波纹，阻带

42、内单调下降；切比雪夫型滤波器通带内为单调下降，阻带内等波纹。调用cheb2ord函数和cheby2函数使切比雪夫型设计变得非常简单。 先用N, wc = cheb2ord(wp, ws, Rp, Rs)求出N和 wc，提供函数cheby2的输入变元，再由B, A = cheby2(N, Rp, wc)设计切比雪夫型数字滤波器。B和A分别为H(z)的分子和分母多项式系数。 对切比雪夫型滤波器，同样有相应的工具箱函数cheb1ord和cheby1。 第61页/共68页62例7.32 椭圆带通数字滤波器设计 设计一个椭圆带通数字滤波器，要求如下 Rp = 1 dB， Rs = 60 dB 解：调用ellipor