离散傅里叶变换(DFT)

1、1一、一、DFT的定义的定义 设设x(n)是一个长度为是一个长度为M的有限长序列，的有限长序列， 则定义则定义x(n)的的N点离散傅里叶变换为点离散傅里叶变换为 X(k)的离散傅里叶逆变换为的离散傅里叶逆变换为 式中，式中， ，N称为称为DFT变换区间长度，变换区间长度，10( ) ( )( ), k=0, 1, , N-1 (3.1.1)NknNnX kDFT x nx n W 101( )( )( ), 0, 1, , -1 (3.1.2)NknNkx nIDFT X kX k WnNN 7.3 离散傅里叶变换（离散傅里叶变换（DFT）2jNNWe NM 2二、二、DFT和和Z变换的关系变

2、换的关系 设序列设序列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)的傅里叶变的傅里

3、叶变换换X(ej)在区间在区间0，2 上的上的N点等间隔采样点等间隔采样。3 实际上，实际上， 任何周期为任何周期为N的周期序列的周期序列 都可以看都可以看作长度为作长度为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.1

4、0)NR 有限长序列有限长序列x(n)的离散傅里叶变换的离散傅里叶变换X(k)，正好是，正好是x(n)的的周期延拓序列周期延拓序列x(n)N的离散傅里叶级数系数的离散傅里叶级数系数 的主值序的主值序列列，即，即X(k) 4图图 3.1.2 有限长序列及其周期延拓有限长序列及其周期延拓 57.3 离散傅里叶变换（离散傅里叶变换（DFT） 定义定义DFT：用类似于例用类似于例7.9中的方法，可把（中的方法，可把（7.3）式写成矩阵乘法运算。）式写成矩阵乘法运算。Xk=xn*Wnk其中，其中，xn为序列行向量，为序列行向量，Wnk是一是一 NN 阶方阵，阶方阵， 而而 称为旋转因子。称为旋转因子。W

5、nk 的的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 6用矩阵乘法计算用矩阵乘法计算N点点DFT的程序如下的程序如下: MATLAB程序程序q73a.m%用矩阵乘法计算用矩阵乘法计算N点点DFTclear;close allxn=input(请输入序列请输入序列x= );N = length(xn); %n=0:N-1; k=n; nk

6、=n*k;% 生成生成NN方阵方阵WN=exp(-j*2*pi/N); % 旋转因子旋转因子Wnk=WN.nk; % 产生旋转因子矩阵产生旋转因子矩阵Xk=xn*Wnk; % 计算计算N点点DFT 这种方法计算这种方法计算DFT概念清楚、编程简单，但概念清楚、编程简单，但占用内存大、运行占用内存大、运行速度低速度低，所以不实用。，所以不实用。MATLAB基础部分提供了基础部分提供了fft、ifft、fft2、ifft2等等快速计算傅里叶变换的函数，使等等快速计算傅里叶变换的函数，使DFT的运算速度量提高了若的运算速度量提高了若干数量级，在后面的例题中均直接调用这些函数。干数量级，在后面的例题中

7、均直接调用这些函数。7例例7.15 序列的离散傅立叶变换序列的离散傅立叶变换 求复正弦序列求复正弦序列 余弦序列余弦序列 正弦序列正弦序列的离散傅立叶变换，分别按的离散傅立叶变换，分别按N =16和和N =8进行计算。进行计算。绘出幅频特性曲线，绘出幅频特性曲线， 进行比较讨论。进行比较讨论。 解：直接产生序列解：直接产生序列x1n、x2n和和x3n，调用，调用fft函数求解函数求解 j81( )e( )nNx nRn 2( )cos( )8NxnnRn 3( )sin( )8Nx nnRn 8%第七章例第七章例7.15程序程序q715 % DFT计算计算clear;close allN=16

8、;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点点DFTx1(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,

9、N); %计算计算N点点DFTx3(n)Xk3=fft(x3n,N1); %计算计算N1点点DFTx1(n)9%绘图绘图subplot(2,3,1);stem(n,abs(X1k),.);title(16点点 DFTx1(n);xlabel(k);ylabel(|X1(k)|)subplot(2,3,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

10、: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(|X2(k)|)subplot(2,3,6);stem(k,abs(Xk3),.);title(8点点 DFTx3(n);xlabel(k);ylabel(|X3(k)|) 10 在截取在截取16点时，得到的是完整的余弦波形；而截取点时，得到的是完整的余弦波形；而截取8点时，点时，得到的是半

11、截的余弦波形，当然有大量的谐波成分。得到的是半截的余弦波形，当然有大量的谐波成分。11例例7.16 验证验证N点点DFT的物理意义的物理意义（1）绘出幅频曲线和相频曲线。绘出幅频曲线和相频曲线。（2）计算并图示）计算并图示x(n)的的8点点DFT。（3）计算并图示）计算并图示x(n)的的16点点DFT。解解: 序列序列x(n)的的点点DFT的物理意义是的物理意义是 在在0，2 上进行上进行点等间隔采样。点等间隔采样。程序先密集采样，绘制出幅频曲线图。然后再分别做程序先密集采样，绘制出幅频曲线图。然后再分别做8点和点和16点点DFT来验证这个采样关系。来验证这个采样关系。 程序略。程序略。441

12、( )( ),() ( )1jjjex nR nX eFT x ne 求求得得，()jX e 12图图7.16 离散傅里叶变换与傅里叶变换的采样关系离散傅里叶变换与傅里叶变换的采样关系13频域采样定理频域采样定理 频域采样是指对有限长序列的傅氏变换在频域离散频域采样是指对有限长序列的傅氏变换在频域离散化得到化得到X(k)的过程。的过程。 本节讨论两个基本问题：本节讨论两个基本问题：a：频域采样（：频域采样（DFT）不失真的条件，即由）不失真的条件，即由 X (k)不失不失真地恢复真地恢复x(n)的条件的条件b：用：用X (k)表示表示X (z)和和 的插值公式的插值公式(内插公式内插公式)()

13、jX e 14图图3-3-1 时域恢复示意图时域恢复示意图 结论：若序列长度为结论：若序列长度为L，频域采样点数（，频域采样点数（DFT的长度）的长度）为为N，且，且LN，则频域采样后可不失真地恢复原序，则频域采样后可不失真地恢复原序列列 ；但；但若若LN，则频域采样后不能不失真地恢复，则频域采样后不能不失真地恢复原序列原序列 。( )x n( )x n15例例7.17 频域与时域采样对偶性频域与时域采样对偶性 （1）产生三角波序列）产生三角波序列 （2）对）对M = 40，计算，计算x(n)的的64点点DFT，并图示，并图示x(n)和和 X(k) = DFTx(n)，k = 0, 1, ,

14、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 ， ，3 31 1132( )xn16 % 第七章例第七章例7.17程序程序q717 % 时域与频域采样的对偶性验证时域与频域采样的对偶性验证 clear;close all M=40;N=64;n=0:M; %产生

15、产生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(n)的周期延拓序列的周期延拓序列 subplot(3,2,1);stem(n,xn,.) title(40点三角波序列点三角波序列x(n);xla

16、bel(n);ylabel(x(n)17 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,.) title(32点点IDFTX1(k);xlabel(n);ylabel(x1(n) subplot(3,1,3);s

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

18、证。1( )IDFT( )x nX k 20用用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如果如果21 由此可见，由此可见， 循环卷积既可在时域直接计算，循环卷积既可在时域直接计算， 也可以按也可以按照图照图3.4.1所示的计算框图，所示的计算框图， 在频域计算。在频域计算。 由于由于DFT有快有快速算法速算法FFT， 当当N很大时

19、，很大时， 在频域计算的速度快得多，在频域计算的速度快得多， 因因而常用而常用DFT(FFT)计算循环卷积。计算循环卷积。 图图 3.4.1 用用DFT计算循环卷积计算循环卷积 22 在实际应用中，需要计算两个序列的线性卷积，在实际应用中，需要计算两个序列的线性卷积， 与与计算循环卷积一样，计算循环卷积一样， 为了提高运算速度，为了提高运算速度， 也希望用也希望用DFT(FFT)计算线性卷积。计算线性卷积。 而而DFT只能直接用来计算循环只能直接用来计算循环卷积，卷积， 为此导出线卷积和循环卷积之间的关系以及循环卷为此导出线卷积和循环卷积之间的关系以及循环卷积与线性卷积相等的条件。积与线性卷积

20、相等的条件。 假设假设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) 23 其中，其中，LmaxN, M， 所以所以1010( )()()( )() ()( )NcLmqNLqmy nh mx nmqL Rnh m x nmqL Rn ( )(),Lqx nx nqL 对照式对照式(3.4.

21、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为周期进行周期延拓才无混叠现象。为周期进行周期延拓才无混叠现象。24图图 3.4.2 线性卷积与循环卷积线性卷积与循环卷积 0123451234h(n) x(n)nL 60123451234nL 8

22、67h(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 1025 利用利用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

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

24、卷积速度大大提高。框图如下: L点FETL点FETL点FETL点FFTL点FFTL 点IFFT28 % 第七章例第七章例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

25、L=pow2(nextpow2(M+N-1);%取取L为大于等于且最接近为大于等于且最接近（N+M-1）的）的2的正次幂的正次幂 tic, %快速卷积计时开始快速卷积计时开始29 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,.);

26、ylabel(h(n) subplot(2,1,2);ny=1:L; stem(ny,real(yn),.);ylabel(y(n) tic, yn=conv(xn,hn);%直接调用函数直接调用函数conv计算卷积与快速卷积计算卷积与快速卷积比较比较 toc3031 4. 用用DFT进行谱分析的误差问题进行谱分析的误差问题 DFT(实际中用实际中用FFT计算计算)可用来对连续信可用来对连续信号和数字信号进行谱分析。号和数字信号进行谱分析。 (1) 混叠现象。混叠现象。 (2) 栅栏效应。栅栏效应。 (3) 截断效应。截断效应。 根据傅里叶变换的频域卷积根据傅里叶变换的频域卷积定理有定理有 (

27、)1() ( )()()21()()2jjjNjjNY eFT y nX eReX eRed 32 幅度谱幅度谱RN()曲线如图曲线如图3.4.11所示所示(RN（）以以2为周期，为周期， 只画低频部分只画低频部分)。 图中，图中，|2/N的的部分称为主瓣，部分称为主瓣， 其余部分称为旁瓣。其余部分称为旁瓣。 例如，例如， x(n)=cos(0n), 0=/4其频谱为其频谱为 其中其中 1()2() ( )sin(/ 2)()( )( )sin(/ 2)jNjjjNNNX eFT x nNReFT RneRe () (2)(2)44jlX ell 33 图图 3.4.11 矩形窗函数的幅度谱矩

28、形窗函数的幅度谱 34图图 3.4.12 加矩形窗前后的频谱加矩形窗前后的频谱 cos()4n 泄露泄露谱间干扰谱间干扰35例例7.19 用用DFT求连续信号频谱求连续信号频谱 在计算机上用在计算机上用DFT对模拟信号进行谱分析时，只能以有限对模拟信号进行谱分析时，只能以有限大的采样频率大的采样频率fs对模拟信号采样有限点样本序列（等价于对模拟信号采样有限点样本序列（等价于截取模拟信号一段进行采样）作截取模拟信号一段进行采样）作DFT变换，得到模拟信号变换，得到模拟信号的近似频谱。其误差主要来自以下因素：的近似频谱。其误差主要来自以下因素： 截断效应（频谱泄露和谱间干扰）截断效应（频谱泄露和谱

29、间干扰） 频谱混叠失真频谱混叠失真 因素因素使谱分辨率（能分辨开的两根谱线间的最小间距）使谱分辨率（能分辨开的两根谱线间的最小间距）降低，并产生谱间干扰；降低，并产生谱间干扰； 因素因素使折叠频率（使折叠频率（ fs / 2）附近的频谱产生较大失真。）附近的频谱产生较大失真。 36例例7.19 用用DFT求连续信号频谱求连续信号频谱 加大截取长度加大截取长度Tp可可提高频率分辨率提高频率分辨率；选择合适的窗函数选择合适的窗函数可可降降低谱间干扰低谱间干扰；而；而频谱混叠失真频谱混叠失真要通过要通过提高采样频率提高采样频率fs和（或）和（或）预滤波预滤波（在采样之前滤除折叠频率以外的频率成分）来

30、改善。（在采样之前滤除折叠频率以外的频率成分）来改善。 编写程序编写程序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 37 % 第七章例第七章例7.19程序程序q719 % 用用DFT作谱分析作谱分析 clear;close all fs=400; T=1/fs;%采样频率为采样频率为400Hz Tp=

31、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实现实现 fk=0:4095/4096/T; subplot(3,2,2*m-1) pl

32、ot(fk,abs(Xk)/max(abs(Xk);ylabel(st(m,:) if m=1 title(矩形窗截取矩形窗截取);end end38 %加加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=0:4095/4096/T; su

33、bplot(3,2,2*m) plot(fk,abs(Xk)/max(abs(Xk);ylabel(st(m,:) if m=1 title(Hamming窗截取窗截取);end end39图图7.1940例例7.19 用用DFT求连续信号频谱求连续信号频谱 如图如图7.19所示。图中所示。图中X1(jf)，X4(jf)和和X8(jf)分别表示分别表示Tp = 0.04s，0.04*4s和和0.04*8s时的谱分析结果。由图可见，由时的谱分析结果。由图可见，由于截断使原频谱中的于截断使原频谱中的单频谱线展宽单频谱线展宽（又称之为（又称之为泄漏泄漏），）， Tp 越大泄漏越小，频率分辨率越高越大

34、泄漏越小，频率分辨率越高。 Tp = 0.04s时，时，25Hz与与50Hz两根谱线已分辨不清了。所以两根谱线已分辨不清了。所以实际谱分析的实际谱分析的截取时间截取时间Tp是由频率分辨率决定的是由频率分辨率决定的。 另外，在本应为零的频段上出现了一些参差不齐的小谱包另外，在本应为零的频段上出现了一些参差不齐的小谱包（称为（称为谱间干扰谱间干扰）。）。谱间干扰的大小取决于加窗的类型谱间干扰的大小取决于加窗的类型。用矩形窗比用用矩形窗比用Hamming窗的频率分辨率高（泄漏小），窗的频率分辨率高（泄漏小），但谱间干扰刚好相反。但谱间干扰刚好相反。 417.5 FIR数字滤波器设计数字滤波器设计 滤

35、波器的特性指标滤波器的特性指标 用绝对值用绝对值1,1,2 2表示表示; ; 用分贝用分贝Rp,RsRp,Rs表示表示1-11p1012s101120log0,120log01RR 42（1）窗函数法设计）窗函数法设计FIR滤波器滤波器 先根据先根据 c和和N求出相应的理想滤波器单位脉冲响求出相应的理想滤波器单位脉冲响应应hd(n)。 第二步要选择合适的窗函数第二步要选择合适的窗函数w(n)来截取来截取hd(n) 的适的适当长度（即阶数），以保证实现要求的阻带衰减；最当长度（即阶数），以保证实现要求的阻带衰减；最后得到后得到FIR滤波器单位脉冲响应滤波器单位脉冲响应h(n)=hd(n).*w(

36、n),即其即其系数。系数。 )()(sin(de )e (21)(cjjdccnnHnhnd43（2）等波纹最佳一致逼近法）等波纹最佳一致逼近法（2）等波纹最佳一致逼近法）等波纹最佳一致逼近法: 信号处理工具箱采用信号处理工具箱采用remez算法实现线性相位算法实现线性相位FIR数数字滤波器的等波纹最佳一致逼近设计。其优点是，设计指字滤波器的等波纹最佳一致逼近设计。其优点是，设计指标相同时，使滤波器阶数最低；或阶数相同时，使通带最标相同时，使滤波器阶数最低；或阶数相同时，使通带最平坦，阻带最小衰减最大；通带和阻带均为等波纹形式，平坦，阻带最小衰减最大；通带和阻带均为等波纹形式，最适合设计片段常

37、数特性的滤波器。其调用格式如下：最适合设计片段常数特性的滤波器。其调用格式如下： b = remez(N, f, m, w, ftype)其中其中N由由remezord函数求出函数求出: N, fo, mo, w = remezord(f, m, dev, Fs)输入变元输入变元dev为各逼近频段允许的波纹振幅。为各逼近频段允许的波纹振幅。remez函数可直接调用函数可直接调用remezord返回的参数如下：返回的参数如下： b=remez(N, fo, mo, w) 44例例7.25 窗函数法设计数字滤波器窗函数法设计数字滤波器 分别用矩形窗和分别用矩形窗和Hamming窗设计线性相位窗设计

38、线性相位FIR低通低通滤波器。要求通带截止频率滤波器。要求通带截止频率 c = /4，单位脉冲响应，单位脉冲响应h(n)的长度的长度N = 21。绘出。绘出h(n)及其幅频响应特性曲线。及其幅频响应特性曲线。 先求出相应的理想滤波器（本例应为理想低通）单先求出相应的理想滤波器（本例应为理想低通）单位脉冲响应位脉冲响应hd(n)，再根据阻带最小衰减选择合适的，再根据阻带最小衰减选择合适的窗函数窗函数w(n)，最后得到，最后得到FIR滤波器单位脉冲响应滤波器单位脉冲响应h(n)=hd(n).*w(n)。 45例例7.25 窗函数法设计数字滤波器窗函数法设计数字滤波器 本题中，本题中， c = /4

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

40、(n)n图图7.25-1图图7.25-247工具箱设计函数工具箱设计函数fir1和和fir2 MATLAB提供了提供了基于窗函数法的基于窗函数法的FIR滤波器设滤波器设计函数计函数fir1和和fir2,其功能及用法如下。其功能及用法如下。 fir1功能：标准频率响应形状。功能：标准频率响应形状。格式：格式：b=fir1(N, wc, ftype, window)。当当wc=wc1，wc2时，得到的是带通滤波器。时，得到的是带通滤波器。 当当ftype=high时，设计高通时，设计高通FIR滤波器；滤波器； 当当ftype=stop时，设计带阻时，设计带阻FIR滤波器。滤波器。 fir2功能：任

41、意频率响应形状。功能：任意频率响应形状。格式：格式：b = fir2(N, f, m, window) 48例例7.26 窗函数法设计带通滤波器窗函数法设计带通滤波器 使用使用fir1函数函数b = fir1(N, wc, window) 编程编程参数参数 c为行向量为行向量 c = lp/ ， hp/ 根据阻带最小衰减根据阻带最小衰减Rs = 60dB选择窗函数类型和阶次。选择窗函数类型和阶次。可以查上面列出的可以查上面列出的“窗函数设计滤波器时的阶数选窗函数设计滤波器时的阶数选择表择表”。选。选blackman窗，其滤波器阻带最小衰减可窗，其滤波器阻带最小衰减可达到达到74dB，其窗口长度

42、，其窗口长度M由过渡带宽度由过渡带宽度B = 0.15 决决定，定，Blackman窗设计的滤波器过渡带宽度为窗设计的滤波器过渡带宽度为12 /M，故故M取取80。因。因M = N+1，所以滤波器阶数，所以滤波器阶数N = 79。 49例例7.27 用用remez函数设计低通滤波器函数设计低通滤波器 解：先由题意计算设计参数解：先由题意计算设计参数 f = 1/4,5/16，m = 1,0; dev的计算稍复杂一些，由于的计算稍复杂一些，由于 所以所以 有了这几个参数就可以调用有了这几个参数就可以调用remezord和和remez函函数了数了.)2(devlog(20 , ) 1 (dev1)

43、 1 (dev120lospAgRpps/20/20(/20)dev(1)(101)/(101) , dev(2)10RRA50例例7.27 用用remez函数设计低通滤波器函数设计低通滤波器 横线为横线为-3dB，两，两条竖线分别位于条竖线分别位于频率频率 /4和和5 /16。显然，通带指标显然，通带指标稍有富裕，过渡稍有富裕，过渡带宽度和阻带最带宽度和阻带最小衰减刚好满足小衰减刚好满足指标要求。指标要求。 程序输出的幅频特性程序输出的幅频特性51例例7.28 remez函数设计高通滤波器函数设计高通滤波器 观察等波纹逼近法中加权系数观察等波纹逼近法中加权系数w( )及滤波器阶数及滤波器阶数

44、N的作用和的作用和影响。期望逼近的滤波器通带为影响。期望逼近的滤波器通带为3 /4， ，阻带为，阻带为0，23 /32。 解解:在滤波器设计中，技术指标越高，实现滤波器的阶数也在滤波器设计中，技术指标越高，实现滤波器的阶数也就越高。另外，对固定的阶数，通带与阻带指标可以互换，就越高。另外，对固定的阶数，通带与阻带指标可以互换，过渡带宽度与通带波纹和阻带衰减指标可以互换。过渡带宽度与通带波纹和阻带衰减指标可以互换。 取取f = 0, 3/4, 23/32, 1，m = 0, 0, 1, 1。其余参数分三种情况。其余参数分三种情况进行设计，进行设计，N = 30，w = 1, 1；N = 30，w

45、 = 1, 5; N = 60，w = 1, 1。 52例例7.28 remez函数设计高通滤波器函数设计高通滤波器 程序运行结果如图程序运行结果如图 由图可见，由图可见，w较大的频段逼近精度较高；较大的频段逼近精度较高；w较小的频段较小的频段逼近精度较低。逼近精度较低。N较大时逼近精度较高较大时逼近精度较高,N较小时逼近精较小时逼近精度较低度较低 。537.6 IIR数字滤波器设计数字滤波器设计 IIR数字滤波器设计的主要方法是数字滤波器设计的主要方法是先设计低通模拟滤波器，先设计低通模拟滤波器，进行频率变换，将其转换为相应的（高通、带通等）模拟滤进行频率变换，将其转换为相应的（高通、带通等

46、）模拟滤波器，再转换为高通、带通或带阻数字滤波器波器，再转换为高通、带通或带阻数字滤波器。对设计的全。对设计的全过程的各个步骤，过程的各个步骤，MATLAB都提供了相应的工具箱函数，都提供了相应的工具箱函数，使使IIR数字滤波器设计变得非常简单。本节主要结合例题介数字滤波器设计变得非常简单。本节主要结合例题介绍这些绍这些IIR滤波器设计的工具箱函数。滤波器设计的工具箱函数。 IIR数字滤波器的设计步骤由以下的流程图来表示。下面以数字滤波器的设计步骤由以下的流程图来表示。下面以巴特沃斯滤波器设计函数为典型，介绍此流程图中函数的功巴特沃斯滤波器设计函数为典型，介绍此流程图中函数的功能和用法。能和用

47、法。54IIR数字滤波器设计流程图数字滤波器设计流程图模拟低通滤波器原型设计模拟低通滤波器原型设计Buttap,cheb1ap,cheb2apbesselap，ellipap函数函数频率变换频率变换(变为高通，带变为高通，带通，带阻等通，带阻等)lp2lp，lp2hp，lp2bp，lp2bs模拟数字变换模拟数字变换bilinearimpinvar合为一步的设计函数合为一步的设计函数 butter，cheb1，cheb2，ellip， besself求最小阶数求最小阶数Nbuttord, cheb1ordcheb2ord,ellipord设计指标设计指标wp, ws, Rp, RsN，wc滤波器

48、系数滤波器系数B，A滤波器系数滤波器系数B，A55巴特沃斯滤波器设计流程巴特沃斯滤波器设计流程 （1）求最小阶数）求最小阶数N的函数的函数buttord N, wc = buttord (wp, ws, Rp, Rs, s) 根据滤波器指标根据滤波器指标wp，ws，Rp，Rs，求出巴特沃斯模拟滤波器的阶数，求出巴特沃斯模拟滤波器的阶数N及及频率参数频率参数wc，此处，此处wp，ws及及wc均以弧度均以弧度/秒为单位。秒为单位。 （2）得到）得到N后后,调用设计函数调用设计函数buttap z，p，k = buttap(N) 得到得到z, p, k后，很容易求出滤波后，很容易求出滤波器系数器系数

49、B，A。 （3）调用模拟频率变换函数）调用模拟频率变换函数lp2lp Bt, At = lp2lp(B, A, wo) （4）调用模拟数字变换函数）调用模拟数字变换函数 Bd, Ad = bilinear (B, A, Fs) 56集成的数字滤波器设计函数集成的数字滤波器设计函数 把（把（2）、（）、（3）、（）、（4）合为一步的数字滤波器设计函数）合为一步的数字滤波器设计函数butter(N, wc, ftype) B, A = butter (N, wc) 设计低通或带通数字滤波器系数设计低通或带通数字滤波器系数B，A（当为带通滤波器时，第（当为带通滤波器时，第(1)类函数由类函数由wp

50、= wp1, wp2会自会自动生成动生成wc = w1, w2）。）。 B, A = butter (N, wc, high) 设计高通数字滤波器系数设计高通数字滤波器系数B，A。 B, A = butter (N, wc, stop) 设计带阻数字滤波器系数设计带阻数字滤波器系数B，A。 butter(N, wc, ftype)还有零极增益和状态空间形式，读者可还有零极增益和状态空间形式，读者可用用help命令查阅。命令查阅。57例例7.29 巴特沃斯模拟滤波器设计巴特沃斯模拟滤波器设计 设计一个设计一个低通巴特沃斯模拟滤波器低通巴特沃斯模拟滤波器，指标如下。，指标如下。通带频率：通带频率：

51、fp = 3400Hz，最大衰减：，最大衰减：Rp = 3dB阻带频率：阻带频率：fs = 4000Hz，最小衰减：，最小衰减：As = 40dB 解解:它的系统函数完全由阶数它的系统函数完全由阶数N和和3dB截止频率截止频率 c决定。而决定。而N和和 c是由滤波器设计指标决定的。是由滤波器设计指标决定的。 取取 c = c1，通带指标刚好，阻带指标富裕；，通带指标刚好，阻带指标富裕； 取取 c = c2，则阻带指标刚好，通带指标富裕。，则阻带指标刚好，通带指标富裕。 MATLAB工具箱函数工具箱函数buttord，butter就是根据以上公式就是根据以上公式编写的。因此就无需再记忆这些公式了

52、。编写的。因此就无需再记忆这些公式了。 1)(10 , 1)10(21s0.1sc2210.1pc1pNANR-58模拟转换为数字：脉冲响应不变法模拟转换为数字：脉冲响应不变法 模拟滤波器离散化的基本方法有模拟滤波器离散化的基本方法有脉冲响应不变法脉冲响应不变法和和双线性双线性变换法变换法。 脉冲响应不变法及脉冲响应不变法及impinvar函数函数单极点的单极点的N阶模拟滤波器阶模拟滤波器Ha(s)，用部分分式展开为，用部分分式展开为脉冲响应不变法的数字化结果为脉冲响应不变法的数字化结果为 工具箱函数工具箱函数impinvar可实现以上计算，格式为可实现以上计算，格式为 Bz, Az = im

53、pinvar(B, A, Fs) NkkkssAsH1a)()(NkTskzAzHk11)e1 ()(59模拟转换为数字：双线性变换法模拟转换为数字：双线性变换法 双线性变换法函数双线性变换法函数bilinear 双线性变换法的原理是用双线性变换法的原理是用 代换代换Ha(s)中中的的s值，得到值，得到H(z)。bilinear函数用来实现这个转换。函数用来实现这个转换。其使用格式为其使用格式为 Bz, Az = bilinear(B, A, Fs) 脉冲响应不变法的缺点是存在频率混叠失真。双线脉冲响应不变法的缺点是存在频率混叠失真。双线性变换法可完全消除频率混叠失真，缺点是存在非性变换法可完

54、全消除频率混叠失真，缺点是存在非线性频率失真。线性频率失真。 11112zzTs60例例7.30 模拟低通转换为数字低通模拟低通转换为数字低通 已知一模拟滤波器的系统函数为已知一模拟滤波器的系统函数为 分别用脉冲响应不变法和双线性变换法将分别用脉冲响应不变法和双线性变换法将Ha(s)转转换成数字滤波器系统函数换成数字滤波器系统函数H(z)，并图示，并图示Ha(s)和和H(z)的幅频响应曲线。的幅频响应曲线。 程序中的核心语句是以下两条：程序中的核心语句是以下两条： d,c=impinvar(b,a,Fs) % 用用impinvar 函数离散化函数离散化 f,e=bilinear(b,a,Fs)

55、% 用用bilinear 函数离散化函数离散化10001000)(assH61例例7.30 模拟低通转换为数字低通模拟低通转换为数字低通 图形结果如图图形结果如图7.30所示。由图所示。由图7.30（b）可见，对）可见，对脉冲响应不变法，采样频率脉冲响应不变法，采样频率Fs越高（越高（T越小），越小），混叠越小；由图混叠越小；由图7.30（c）可见，对双线性变换法，）可见，对双线性变换法，无频率混叠，但存在非线性失真。无频率混叠，但存在非线性失真。62例例7.31 切比雪夫切比雪夫数字滤波器设计数字滤波器设计 解：切比雪夫解：切比雪夫型滤波器通带内为等波纹，阻带内单调下降；型滤波器通带内为等波

56、纹，阻带内单调下降；切比雪夫切比雪夫型型滤波器通带内为单调下降，阻带内等波纹。滤波器通带内为单调下降，阻带内等波纹。调用调用cheb2ord函数和函数和cheby2函数使切比雪夫函数使切比雪夫型设计变得型设计变得非常简单。非常简单。 先先用用N, wc = cheb2ord(wp, ws, Rp, Rs)求出求出N和和 wc，提供函数提供函数cheby2的输入变元，再的输入变元，再由由B, A = cheby2(N, Rp, wc)设计切比雪夫设计切比雪夫型数字滤波器型数字滤波器。B和和A分别为分别为H(z)的分的分子和分母多项式系数。子和分母多项式系数。 对对切比雪夫切比雪夫型型滤波器，同样

57、有相应的工具箱函数滤波器，同样有相应的工具箱函数cheb1ord和和cheby1。 63例例7.32 椭圆带通数字滤波器设计椭圆带通数字滤波器设计 设计一个设计一个椭圆带通数字滤波器椭圆带通数字滤波器，要求如下，要求如下 Rp = 1 dB， Rs = 60 dB 解：调用解：调用ellipord函数和函数和ellip函数，代入参数函数，代入参数wp = 0.25，0.45；ws = 0.15，0.55；Rp = 1；Rs = 60;即可得到本题的程序即可得到本题的程序 。lplsupus50.55，64例例7.33 高通和带通数字滤波器设计高通和带通数字滤波器设计 用双线性变换法从用双线性变换