MATLAB基础与应用教程(人民邮电出版社-蔡旭辉)第七章b.ppt

1、1,一、DFT的定义 设x(n)是一个长度为M的有限长序列， 则定义x(n)的N点离散傅里叶变换为 X(k)的离散傅里叶逆变换为 式中， ，N称为DFT变换区间长度,7.3 离散傅里叶变换（DFT,2,二、DFT和Z变换的关系 设序列x(n)的长度为N， 其Z变换和DFT分别为,比较上面二式可得关系式,3.1.3）式表明序列x(n)的N点DFT是x(n)的Z变换在单位圆上的N点等间隔采样； （3.1.4）式则说明X(k)为x(n)的傅里叶变换X(ej)在区间0，2 上的N点等间隔采样,3,实际上， 任何周期为N的周期序列 都可以看作长度为N的有限长序列x(n)的周期延拓序列， 即,为了以后叙述

2、方便， 将(3.1.5)式用如下形式表示,三、DFT的隐含周期性,均为整数,有限长序列x(n)的离散傅里叶变换X(k)，正好是x(n)的周期延拓序列x(n)N的离散傅里叶级数系数 的主值序列，即,4,图 3.1.2 有限长序列及其周期延拓,5,7.3 离散傅里叶变换（DFT,定义DFT： 用类似于例7.9中的方法，可把（7.3）式写成矩阵乘法运算。 Xk=xn*Wnk 其中，xn为序列行向量，Wnk是一 NN 阶方阵， 而 称为旋转因子。 Wnk 的MATLAB表示： Wnk = WN.(0:N-1* 0:N-1,6,用矩阵乘法计算N点DFT的程序如下: MATLAB程序q73a.m %用矩阵

3、乘法计算N点DFT clear;close all xn=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的运算速度量提高了若干数量级，在后面的例题中均直接调用这些函数,7,例7.15 序列的离散傅立叶

4、变换,求复正弦序列 余弦序列 正弦序列 的离散傅立叶变换，分别按N =16和N =8进行计算。绘出幅频特性曲线， 进行比较讨论。 解：直接产生序列x1n、x2n和x3n，调用fft函数求解,8,第七章例7.15程序q715 % DFT计算 clear;close all N=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点DFTx1(n) %产生序列x2(n),计算DFTx2(n) x2n=cos(pi

5、*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,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点

6、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(|X2(k)|) subplot(2,3,6

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

8、变换的采样关系,13,频域采样定理,频域采样是指对有限长序列的傅氏变换在频域离散化得到X(k)的过程。 本节讨论两个基本问题： a：频域采样（DFT）不失真的条件，即由 X (k)不失真地恢复x(n)的条件 b：用X (k)表示X (z)和 的插值公式(内插公式,14,图3-3-1 时域恢复示意图,结论：若序列长度为L，频域采样点数（DFT的长度）为N，且LN，则频域采样后可不失真地恢复原序列 ；但若LN，则频域采样后不能不失真地恢复原序列,15,例7.17 频域与时域采样对偶性,1）产生三角波序列 （2）对M = 40，计算x(n)的64点DFT，并图示x(n)和 X(k) = DFTx(n

9、)，k = 0, 1, , 63。 （3）对（2）中所得X(k)在 0，2 上进行32点抽样得 （4）求 的32点IDFT，即 （5）绘出 的波形图，评述它与x(n)的关系,16, 第七章例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点IFFT

10、X1(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);xlabel(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);yl

11、abel(|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);stem(nc,xc,.) title(x1(n)的周期延拓序列);xlabel(n);ylabel(x(mod(n,32) set(gcf,color,w,18,19,由于频域在0，2上的采样点数N（N = 32）小于x(n)的长度M（M = 40），所以，产生时域混叠现象，不能由X1(k)恢复原序列x(n)。只有满足NM时，可由频域采样X1(k)得到原序列x(n)。

12、这就是频域采样定理。对NM的情况，请读者自己编程验证,20,用DFT计算线性卷积,0kL-1,则由时域循环卷积定理有 Y(k)=DFTy(n)=X1(k)X2(k), 0kL-1,如果,21,由此可见， 循环卷积既可在时域直接计算， 也可以按照图3.4.1所示的计算框图， 在频域计算。 由于DFT有快速算法FFT， 当N很大时， 在频域计算的速度快得多， 因而常用DFT(FFT)计算循环卷积,图 3.4.1 用DFT计算循环卷积,22,在实际应用中，需要计算两个序列的线性卷积， 与计算循环卷积一样， 为了提高运算速度， 也希望用DFT(FFT)计算线性卷积。 而DFT只能直接用来计算循环卷积，

13、 为此导出线卷积和循环卷积之间的关系以及循环卷积与线性卷积相等的条件。 假设h(n)和x(n)都是有限长序列， 长度分别是N和M。 它们的线性卷积和循环卷积分别表示如下,3.4.1,3.4.2,23,其中，LmaxN, M， 所以,对照式(3.4.1)可以看出， 上式中,3.4.3,即,3.4.3)式说明，yc(n)等于yl(n)以L为周期的周期延拓序列的主值序列。 yl(n)长度为N+M-1，只有当循环卷积长度L N+M-1时，yl(n)以L为周期进行周期延拓才无混叠现象,24,图 3.4.2 线性卷积与循环卷积,25,利用DFT进行线性卷积的步骤如下： 1.将序列x(n)和h(n)补零延长

14、，使其长度 2.做x(n)和h(n)的长为L点的DFT得到X(k)和H(k)，求它们的积Y(k)=X(k)H(k) 3.求Y(k)的IDFT并取前 N1 点获得线性卷积的结果为,LN1=N+M-1,y(n)=IDFTY(k)，0nN1,26,图 3.4.3 用DFT计算线性卷积框图,27,例7.18 快速卷积,快速卷积就是根据DFT的循环卷积性质，将时域卷积转换为频域相乘，最后再进行IDFT得到时域卷积序列y(n)。其中时域和频域之间的变换均用FFT实现，所以使卷积速度大大提高。框图如下,28, 第七章例7.18程序q718 % 快速卷积计算 clear;close all xn=input(

15、请输入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, %快速卷积计时开始,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

16、点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计算卷积与快速卷积比较 toc,30,31,4. 用DFT进行谱分析的误差问题 DFT(实际中用FFT计算)可用来对连续信号和数字信号进行谱分析。 (1) 混叠现象。 (2) 栅栏效应。 (3) 截断

17、效应。 根据傅里叶变换的频域卷积定理有,32,幅度谱RN()曲线如图3.4.11所示(RN（）以2为周期， 只画低频部分)。 图中，|2/N的部分称为主瓣， 其余部分称为旁瓣。 例如， x(n)=cos(0n), 0=/4其频谱为,其中,33,图 3.4.11 矩形窗函数的幅度谱,34,图 3.4.12 加矩形窗前后的频谱,泄露 谱间干扰,35,例7.19 用DFT求连续信号频谱,在计算机上用DFT对模拟信号进行谱分析时，只能以有限大的采样频率fs对模拟信号采样有限点样本序列（等价于截取模拟信号一段进行采样）作DFT变换，得到模拟信号的近似频谱。其误差主要来自以下因素： 截断效应（频谱泄露和谱

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

19、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实现 fk=0:4095/4096/T; subplot(3,

20、2,2*m-1) plot(fk,abs(Xk)/max(abs(Xk);ylabel(st(m,:) if m=1 title(矩形窗截取);end end,38,加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; subplot(3,2,2*m) plot(

21、fk,abs(Xk)/max(abs(Xk);ylabel(st(m,:) if m=1 title(Hamming窗截取);end end,39,图7.19,40,例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两根谱线已分辨不清了。所以实际谱分析的截取时间Tp是由频率分辨率决定的。 另外，在本应为零的频段上出现了一些参差不齐的

22、小谱包（称为谱间干扰）。谱间干扰的大小取决于加窗的类型。用矩形窗比用Hamming窗的频率分辨率高（泄漏小），但谱间干扰刚好相反,41,7.5 FIR数字滤波器设计,滤波器的特性指标 用绝对值1,2表示; 用分贝Rp,Rs表示,42,1）窗函数法设计FIR滤波器,先根据c和N求出相应的理想滤波器单位脉冲响应hd(n)。 第二步要选择合适的窗函数w(n)来截取hd(n) 的适当长度（即阶数），以保证实现要求的阻带衰减；最后得到FIR滤波器单位脉冲响应h(n)=hd(n).*w(n),即其系数,43,2）等波纹最佳一致逼近法,2）等波纹最佳一致逼近法: 信号处理工具箱采用remez算法实现线性相位

23、FIR数字滤波器的等波纹最佳一致逼近设计。其优点是，设计指标相同时，使滤波器阶数最低；或阶数相同时，使通带最平坦，阻带最小衰减最大；通带和阻带均为等波纹形式，最适合设计片段常数特性的滤波器。其调用格式如下： 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窗设计

24、线性相位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，N = 21，所以线性相位理想低通滤波器的单位脉冲响应为: 为了满足线性相位FIR滤波器条件h(n) = h(N-1-n)，要求 = (N-1)/2 = 10。 信号处理工具箱中有窗生成函数boxcar，hamming，hanni

25、ng和blackman等,46,例7.25 窗函数法设计数字滤波器,对两种窗函数的设计结果分别如右图7.25-1和图7.25-2所示,图7.25-1,图7.25-2,47,工具箱设计函数fir1和fir2,MATLAB提供了基于窗函数法的FIR滤波器设计函数fir1和fir2,其功能及用法如下。 fir1功能：标准频率响应形状。 格式：b=fir1(N, wc, ftype, window)。 当wc=wc1，wc2时，得到的是带通滤波器。 当ftype=high时，设计高通FIR滤波器； 当ftype=stop时，设计带阻FIR滤波器。 fir2功能：任意频率响应形状。 格式：b = fir

26、2(N, f, m, window,48,例7.26 窗函数法设计带通滤波器,使用fir1函数b = fir1(N, wc, window) 编程 参数c为行向量c = lp/，hp/ 根据阻带最小衰减Rs = 60dB选择窗函数类型和阶次。可以查上面列出的“窗函数设计滤波器时的阶数选择表”。选blackman窗，其滤波器阻带最小衰减可达到74dB，其窗口长度M由过渡带宽度B = 0.15 决定，Blackman窗设计的滤波器过渡带宽度为12/M，故M取80。因M = N+1，所以滤波器阶数N = 79,49,例7.27 用remez函数设计低通滤波器,解：先由题意计算设计参数 f = 1/4

27、,5/16，m = 1,0; dev的计算稍复杂一些，由于 所以 有了这几个参数就可以调用remezord和remez函数了,50,例7.27 用remez函数设计低通滤波器,横线为-3dB，两条竖线分别位于频率/4和5/16。显然，通带指标稍有富裕，过渡带宽度和阻带最小衰减刚好满足指标要求,程序输出的幅频特性,51,例7.28 remez函数设计高通滤波器,观察等波纹逼近法中加权系数w()及滤波器阶数N的作用和影响。期望逼近的滤波器通带为3/4，阻带为0，23/32。 解:在滤波器设计中，技术指标越高，实现滤波器的阶数也就越高。另外，对固定的阶数，通带与阻带指标可以互换，过渡带宽度与通带波纹

28、和阻带衰减指标可以互换。 取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,52,例7.28 remez函数设计高通滤波器,程序运行结果如图 由图可见，w较大的频段逼近精度较高；w较小的频段逼近精度较低。N较大时逼近精度较高,N较小时逼近精度较低,53,7.6 IIR数字滤波器设计,IIR数字滤波器设计的主要方法是先设计低通模拟滤波器，进行频率变换，将其转换为相应的（高通、带通等）模拟滤波器，再转换为高通、带通或带阻数字滤波器。对设计的全过程的各

29、个步骤，MATLAB都提供了相应的工具箱函数，使IIR数字滤波器设计变得非常简单。本节主要结合例题介绍这些IIR滤波器设计的工具箱函数。 IIR数字滤波器的设计步骤由以下的流程图来表示。下面以巴特沃斯滤波器设计函数为典型，介绍此流程图中函数的功能和用法,54,IIR数字滤波器设计流程图,模拟低通滤波器原型设计 Buttap,cheb1ap,cheb2ap besselap，ellipap函数,频率变换(变为高通，带通，带阻等) lp2lp，lp2hp，lp2bp，lp2bs,模拟数字变换 bilinear impinvar,合为一步的设计函数 butter，cheb1，cheb2，ellip，

30、 besself,求最小阶数N buttord, cheb1ord cheb2ord,ellipord,设计指标 wp, ws, Rp, Rs,N，wc,滤波器系数 B，A,滤波器系数 B，A,55,巴特沃斯滤波器设计流程,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后，很容易求出滤波器系数B，A。 （3）调用

31、模拟频率变换函数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 = wp1, wp2会自动生成wc = w1, w2）。 B, A = butter (N, wc, high) 设计高通数字滤波器系数B，A。 B, A = butter (N, wc, s

32、top) 设计带阻数字滤波器系数B，A。 butter(N, wc, ftype)还有零极增益和状态空间形式，读者可用help命令查阅,57,例7.29 巴特沃斯模拟滤波器设计,设计一个低通巴特沃斯模拟滤波器，指标如下。 通带频率：fp = 3400Hz，最大衰减：Rp = 3dB 阻带频率：fs = 4000Hz，最小衰减：As = 40dB 解:它的系统函数完全由阶数N和3dB截止频率c决定。而N和c是由滤波器设计指标决定的。 取c = c1，通带指标刚好，阻带指标富裕； 取c = c2，则阻带指标刚好，通带指标富裕。 MATLAB工具箱函数buttord，butter就是根据以上公式编写

33、的。因此就无需再记忆这些公式了,58,模拟转换为数字：脉冲响应不变法,模拟滤波器离散化的基本方法有脉冲响应不变法和双线性变换法。 脉冲响应不变法及impinvar函数 单极点的N阶模拟滤波器Ha(s)，用部分分式展开为 脉冲响应不变法的数字化结果为 工具箱函数impinvar可实现以上计算，格式为 Bz, Az = impinvar(B, A, Fs,59,模拟转换为数字：双线性变换法,双线性变换法函数bilinear 双线性变换法的原理是用 代换Ha(s)中的s值，得到H(z)。bilinear函数用来实现这个转换。其使用格式为 Bz, Az = bilinear(B, A, Fs) 脉冲响

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

35、（T越小），混叠越小；由图7.30（c）可见，对双线性变换法，无频率混叠，但存在非线性失真,62,例7.31 切比雪夫数字滤波器设计,解：切比雪夫型滤波器通带内为等波纹，阻带内单调下降；切比雪夫型滤波器通带内为单调下降，阻带内等波纹。调用cheb2ord函数和cheby2函数使切比雪夫型设计变得非常简单。 先用N, wc = cheb2ord(wp, ws, Rp, Rs)求出N和 wc，提供函数cheby2的输入变元，再由B, A = cheby2(N, Rp, wc)设计切比雪夫型数字滤波器。B和A分别为H(z)的分子和分母多项式系数。 对切比雪夫型滤波器，同样有相应的工具箱函数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; 即可得到本题的程序,64,例7.33 高通