1、1第第 7章章 MATLAB在数字信号处理中的在数字信号处理中的应用应用2书山有路勤为径，书山有路勤为径，学海无涯苦作舟。学海无涯苦作舟。好好学习，天天向上好好学习，天天向上3引言n数字信号处理中包含离散时间信号与系统的一部分知识。n本章主要研究内容分为两部分：离散时间信号与系统时域分析离散时间信号与系统频域分析n在数字信号处理中还有一个非常大的知识点：数字滤波器设计，留待下学期开设本课程时，可以学习。4一、离散时间信号与系统时域分析5包含内容n信号分析典型信号表示impseq(),stepseq()7.1 序列相加，相乘7.2 序列合成和截取7.3 序列移位和周期延拓7.6离散序列卷积n系统
2、分析7.4系统响应7.5系统线性判定67.1 离散信号的产生及时域处理n时域离散信号用x(n)表示（离散序列）。x x和n同时使用才能完整地表示一个序列：时间变量n（表示采样位置）只能取整数；向量x（表示采样点的幅度）与n一一对应。n由于n序列是按整数递增的,可简单地用其初值ns决定,因为它的终值nf取决于ns 和x的长度length(x),故可写成：n=ns:nf 或n=ns:nslength(x)17n典型信号表示典型信号表示impseq(),stepseq()8例7.1 序列的相加和相乘n给出两个序列x1(n)和x2(n)。x1=0,1,2,3,4,3,2,1,0;n1=-2:6;x2=
3、2,2,0,0,0,-2,-2;n2=2:8;要求它们的和ya及乘积yp。解：编程的思路是把序列长度延拓到覆盖n1和n2的范围，这样才能把两序列的时间变量对应起来，然后进行对应元素的运算。9例程nx1=0,1,2,3,4,3,2,1,0;ns1=-2;%给定x1及ns1nx2=2,2,0,0,0,-2,-2;ns2=2;%给定x2及ns2nnf1=ns1+length(x1)-1;nf2=ns2+length(x2)-1;nny=min(ns1,ns2):max(nf1,nf2);%y(n)的时间变量nxa1=zeros(1,length(ny);xa2=xa1;%延拓序列初始化nxa1(fi
4、nd(ny=ns1)&(ny=ns2)&(ny=0;%产生单位阶跃序列（u(n-n1)）nx1=(n-n1)=0-(n-n1-N)=0%用阶跃序列之差产生矩形序列nx2=(n=n1)&(n=0)&(n=0;%产生输入信号x3(n)ny3=filter(B,A,x3);%对x3(n)的响应nx4=(n=0)&(n=0)&(nx=sym(an*cos(pi*n);Z=ztrans(x);simplify(Z)ans=z/(z+a)()cos()(nunanxn35【例例2 2】试用试用iztransiztrans函数求下列函数的函数求下列函数的z z反变换。反变换。Z=sym(8*z-19)/(z
5、2-5*z+6);x=iztrans(Z);simplify(x)ans=-19/6*charfcn0(n)+5*3(n-1)+3*2(n-1)-19/6*charfcn0(n)+5*3(n-1)+3*2(n-1)charfcn0(n)charfcn0(n)是函数在是函数在MATLABMATLAB符号工具箱中的表示，符号工具箱中的表示，反变换后的函数形式为：反变换后的函数形式为：3|65198)(2zzzzzX)()2335()(619)(11nunnxnn36例7.14 用符号运算工具箱解z变换n解：无限长度时间序列的z变换和逆z变换都属于符号运算的范围。MATLAB的symbolic（符号
6、运算）工具箱已提供了这种函数。如果读者已在计算机上安装了这个工具箱，可以键入以下程序。n MATLAB程序q714.m n其特点是程序的开始要指定符号自变量nsyms z n a N w0%规定z,n,a为符号变量37程序分析：nsyms z n aN w0%规定z,n,a为符号变量ny1=an;%给出y的表示式nY1=ztrans(y1)%求y的z变换nX1=-3*z-1/(2-5*z-1+2*z-2);%给出X的表示式nx1=iztrans(X1)%求X的逆Z变换nsimplify(x);%可显示变换后表达式的形式38n如果信号的如果信号的z z域表示式是有理函数，进行域表示式是有理函数，
7、进行z z反变换的反变换的另一个方法是对进行部分分式展开，然后求各简单另一个方法是对进行部分分式展开，然后求各简单分式的分式的z z反变换反变换.如果如果X(z)X(z)的有理分式表示为：的有理分式表示为：)()(1)(221122110zAzBzazazazbzbzbbzXnnmmrkkikrMkkknNMnnzzCzzAzBzX111101 1)(nMATLABMATLAB实现实现-方法方法2 2（residuezresiduez函数实函数实现反变换）现反变换）39 nMATLABMATLAB信号处理工具箱提供了一个对进行部分分式展开信号处理工具箱提供了一个对进行部分分式展开的函数的函数r
8、esiduezresiduez，其语句格式为：，其语句格式为：nR,P,K=residuez(B,A)R,P,K=residuez(B,A)n其中，其中，B B，A A分别表示分别表示X(z)X(z)的分子与分母多项式的系数向的分子与分母多项式的系数向量，分子与分母多项式按照量，分子与分母多项式按照z-1z-1升幂排列。升幂排列。nR R为部分分式的系数向量；为部分分式的系数向量；nP P为极点向量；为极点向量；nK K为多项式的系数。为多项式的系数。n若若X(z)X(z)为有理真分式，则为有理真分式，则K K为零。为零。40【例【例1 1】试用试用MATLABMATLAB命令进行部分分式展开
9、，并求出其命令进行部分分式展开，并求出其z z反变换。反变换。n解：解：MATLAB源程序为源程序为nB=18;nA=18,3,-4,-1;nR,P,K=residuez(B,A)nR=n 0.3600n 0.2400n 0.40005.0|431818)(321zzzzzX41 nP=n 0.5000n -0.3333n -0.3333nK=n 从运行结果可知，32pp 表示系统有一个二重极点。所以，X(z)的部分分式展开为2111)3330.314.03333.0124.05.0136.0)(zzzzX（)()3333.0)(1(4.0)3333.0(24.0)5.0(36.0)(nunn
10、xnnn42教材中求z的逆变换的方法n对于z变换分式n可以用部分分式法或长除法求其反变换。n用函数residuez可以求出它的极点留数分解其中r,p,k=residuez(B,A)n其反变换为：)1()()2()1()1()()2()1()()()()1(1)1(1NNMMzNAzNAzAAzMBzMBzBBzAzBzY1111)2()1()(1)()2(1)2()1(1)1()()(zkkzNpNrzprzprzAzB()(1)(1)()(2)(2)()(1)()nny nrpu nrpu nkn43注意几点：nr,p,k=residuez(B,A),A首项不能为零，即分母z最高次项不能为0
11、nk只有当分式为假分式时会有值，即对应冲激序列部分（FIR）；真分式，值为零n还有一些细节需要处理：44例7.7 有限序列的z和逆z变换两序列x1=1,2,3，n1=-1:1 及x2=2,4,3,5，n2=-2:1，求出x1与x2及其卷积x的z变换。n解：其z变换可写成n两个多项式乘积可用conv函数来求得。n数组要自己判别。n的起点ns=ns1+ns2=3，终点nf=nf1+nf2=2。n=ns:nf。由x和n即可得出X(z)。12112()23,()2435XzzzXzzzz)()()(21zXzXzX45程序分析nx1=1,2,3;ns1=-1;%设定x1和ns1nnf1=ns1+len
12、gth(x1)-1;%nf1可以算出nx2=2,4,3,5;ns2=-2;%设定x2和ns2nnf2=ns1+length(x1)-1;%nf2可以算出nx=conv(x1,x2)%求出xnn=(ns1+ns2):(nf1+nf2)%求出n4647例7.8 求z多项式分式的逆变换 n设系统函数为，输入例7.7中的x2信号，用z变换计算输出y(n)解：由例7.7可知，故Y(z)=X(z)W(z)=其中nsy=分母分子中z的最高幂次之差。调用 r,p,k=residuez(B,A)，可由B，A求出r，p，k，进而求逆z变换，得 2115.02.223)(zzzzW1225342)(zzzzX()(
13、)nsyzB zA z()(1)(1)()(1)()n nsyy nrpu nnsyknnsy48例7.8 z多项式分式逆变换(续)由程序算出nsy=-1，留数、极点分别为r=-57.7581 和 204.7581p=0.7791 和 0.3209k=-150 -30代入得111)2()1()2(1)2()1(1)1()(zkkzprzprzzY11()57.7 0.779(1)204.7 0.321(1)150(n+1)-30 (n)nny nu nu n 49Z反变换法求解系统输出(编程实现)nY(z)=X(z)W(z)nx=2,4,3,5;nsx=-2;%输入序列及初始时间nnfx=ns
14、x+length(x)-1;%计算序列终止时间nBw=-3;nsbw=-1;%系统函数的分子系数,及z的最高次数nAw=2,-2.2,0.5;nsaw=0;%系统函数的分母系数,及z的最高次数nB=conv(-3,x);%输入与分子z变换的多项式乘积nA=Aw;%分母不变ny(n)=IZTY(Z)nr,p,k=residuez(B,A)%求留数r,极点p及直接项knnf=input(终点时间nf=);%要求键入终点时间nn=min(nsx,nsy):nf;%生成总时间数组n%求无限序列yi和直接序列ydnyi=(r(1)*p(1).(n-nsy)+r(2)*p(2).(n-nsy).*step
15、seq(nsy,n(1),nf);nyd=k(1)*impseq(nsy,n(1),nf)+k(2)*impseq(-1-nsy,n(1),nf);ny=yi+yd;50另外一种方法(filter函数)：nfilter(B,A,x)51例7.9 离散时间傅里叶变换n取周期的正弦信号，作8点采样，求它的连续频谱。然后对该信号进行N个周期延拓，再求它的连续频谱。把N无限增大，比较分析其结果。n解：先求离散傅立叶变换的MATLAB子程序n最后得到X=x*exp(-j*w*n)。n有了子程序，本例就没有什么难度了。NKNNKKnjnjnjnjnjnjnjnjnjNKeeeeeeeeenxnxXXX21
16、2222111211)(,),()(,),(),(12152DTFT程序：nfunction X=dtft(x,w)n%计算离散时间傅立叶变换n%X=dtft(x,n,w),%X=在w频率点上的DTFT数组,%x=沿n的有限长度序列,%n=样本位置向量%w=频率点位置向量nn=1:length(x);ewn=exp(-n*w*i);X=x*ewn;53正弦信号x0的DTFT变换nx0=sin(2*pi*1:8/8)*5;%x0是8点行向量ndt=2*pi/8;nw=linspace(-2*pi,2*pi,1000)/dt;%w是1000点行向量nX0=dtft(x0,w)*dt;%求得频率响应
17、X054正弦信号延拓N个周期后，DTFTnN=4;%延拓周期数nx1=reshape(x0*ones(1,N),1,N*length(x0);%延拓后的时域信号x1nX1=dtft(x1,w)*dt;%求x1的频率响应X1n%延拓100次后n%x1=reshape(x0*ones(1,100),1,100*length(x0);%延拓后的时域信号x1n%X1=dtft(x1,w)/100*dt;%求x1的频率响应X155重复无穷次ndisp(重复无穷次的八点信号的离散时间傅立叶变换-傅立叶级数)npause,X2=fft(x0*dt);%离散傅立叶变换nw1=2*pi*0:length(x0)
18、-1/length(x0);%离散频点向量nsubplot(3,1,3),stem(-w1,w1,abs(X2),abs(X2),grid,naxis(min(w),max(w),0,max(abs(X2),grid56例7.9 离散时间傅里叶变换2n 程序运行结果执行程序q709并按提示键入N=4，所得图形如图7.10所示。N取得愈大，其峰值愈大，宽度愈窄。当N取得很大时，会出现内存不足的问题，这是用矩阵乘法做傅里叶变换的缺点。另外，因为那时峰值点处的宽度很窄，也会出现所选频点对不上峰值点的问题。所以对于N无限增大的情况，必须用fft函数来求。这时用连续频谱也没有意义了。这里用同样的横坐标把
19、几种频谱进行对比，使读者更好地理解其关系。57程序显示结果：58例7.10 时域采样频率与频谱混叠n分别以采样频率fs=1000Hz，400Hz和200Hz对xa(t)进行等间隔采样，计算并图示三种采样频率下的采样信号及其幅频特性n解：程序分别设定4种采样频率fs=10kHz，1kHz，400Hz和200Hz，对xa(t)进行采样，得到采样序列xa(t)，xa1(n)，xa2(n)，xa3(n)，画出其幅度频谱。采样时间区间均为0.1秒。为了便于比较，画出了幅度归一化的幅频曲线，如图7.11所示。59例7.10 采样频率与频谱混叠(续)由于由以上关系式可见，采样信号的频谱函数是原模拟信号频谱函
20、数的周期延拓，延拓周期为2/T。如果以频率f为自变量（=2f），则以采样频率fs=1/T为延拓周期。对频带限于fc的模拟信号xa(t)，只有当fs2fc时，采样后 才不会发生频谱混叠失真。这就是著名的采样定理 knTTTkXTnxnxXn2j1e)()(FT)e(ajj)e(j TX程序要点：nt=0:1/fs:0.1;%fs代不同的取样频率 nxa=exp(-a*t).*sin(b*t);nk=0:511;f=fs*k/512;n%由wk=2k/512=2fT求得模拟频率fnXa=dtft(xa,2*pi*k/512);%近似模拟信号频谱60结果161结果26263例7.12 梳状滤波器零极
21、点和幅特性n梳状滤波器系统函数有如下两种类型。n FIR型：n IIR型：n freqz 数字滤波器频率特性计算和绘制函数n zplane H(z)的零-极点图绘制。n解：调用函数freqz和zplane 很容易写出程序q712.m。NzzH1)(1NNNzazzH11)(264程序分析nfreqzfreqz(b,a,w)nzplane65频率响应nb=1,0,0,0,0,0,0,0,-1;na0=1;na1=1,0,0,0,0,0,0,0,-(0.8)8;na2=1,0,0,0,0,0,0,0,-(0.9)8;na3=1,0,0,0,0,0,0,0,-(0.98)8;nH,w=freqz(b
22、,a0);nH1,w1=freqz(b,a1);nH2,w2=freqz(b,a2);nH3,w3=freqz(b,a3);66零极点及绘图nsubplot(2,2,1);zplane(b,a0);title(FIR梳状滤波器零点图);nsubplot(2,2,2);zplane(b,a1);title(IIR梳状滤波器零、极点图,a=0.8);nsubplot(2,2,3);plot(w/pi,abs(H);title(FIR梳状滤波器幅频响应曲线);nylabel(幅度);xlabel(/);nsubplot(2,2,4);plot(w1/pi,abs(H1);title(IIR梳状滤波器
23、幅频响应曲线,a=0.8);nylabel(幅度);xlabel(/);nfigure(2);nsubplot(2,2,1);zplane(b,a2);title(IIR梳状滤波器零、极点图,a=0.9);nsubplot(2,2,2);zplane(b,a3);title(IIR梳状滤波器零、极点图,a=0.98);nsubplot(2,2,3);plot(w2/pi,abs(H2);title(IIR梳状滤波器幅频响应曲线,a=0.9);nylabel(幅度);xlabel(/);nsubplot(2,2,4);plot(w3/pi,abs(H3);title(IIR梳状滤波器幅频响应曲线
24、,a=0.98);nylabel(幅度);xlabel(/)67结果显示68结果显示69例7.13 低通滤波及时域卷积定理 n输入信号 x(n)=cos(0.04n)+cos(0.08n)+cos(0.4n)+0.3(n)，0n63 通过低通滤波器，计算滤波器对x(n)的响应输出y(n)，并图示x(n)和y(n)，观察滤波效果。n解：如前所述，只要求出H(z)=B(z)/A(z)的分子和分母多项式系数向量B和A，则可调用滤波器直接型实现函数filter对输入信号x(n)进行滤波。y=filter(B,A,x)70程序要点n输入信号构造%产生输入信号产生输入信号x(n)x(n)nn=0:255;
25、N=4096;n=0:255;N=4096;nx=cos(0.04x=cos(0.04*pipi*n)+cos(0.08n)+cos(0.08*pipi*n)+cos(0.n)+cos(0.4 4*pipi*n);n);nw=randn(size(x);%w=randn(size(x);%产生正态零均值噪声产生正态零均值噪声nx=x+0.3x=x+0.3*w;w;71程序要点n%求H(z)分子分母多项式系数向量B和Anb=1,2,1;%(1+z-1)2 的展开系数nB=0.0003738*conv(conv(b,b),b);%嵌套调用卷积函数convna1=1,-1.2686,0.7051;n
26、a2=1,-1.0106,0.3583;na3=1,-0.9044,0.2155;nA=conv(conv(a1,a2),a3);72程序要点n滤波及图形n%对x(n)滤波ny=filter(B,A,x);n%绘图nX=fft(x,N);nY=fft(y,N);nsubplot(3,2,1);stem(x,.)naxis(0,max(n)/4,min(x),max(x);line(0,max(n),0,0)ntitle(输入信号x(n);xlabel(n);ylabel(x(n)nsubplot(3,2,5);stem(y,.)naxis(0,max(n)/4,min(y),max(y);li
27、ne(0,max(n),0,0)ntitle(输出信号 y(n);xlabel(n);ylabel(y(n)nk=0:N-1;f=2*k/N;nsubplot(3,2,2);plot(f,abs(X)ntitle(输入信号x(n)的幅频曲线);xlabel(/);ylabel(|FTx(n)|)naxis(0,0.5,0,max(abs(X);nsubplot(3,2,6);plot(f,abs(Y)ntitle(输出信号y(n)的幅频曲线);xlabel(/);ylabel(|FTy(n)|)naxis(0,0.5,0,max(abs(Y);73结果显示：74程序要点nh,f=freqz(B
28、,A,N,whole);nfigure(2)nsubplot(3,2,1);plot(f/pi,abs(h)ntitle(滤波器幅频响应曲线);xlabel(/);ylabel(H幅度)naxis(0,0.5,0,max(abs(h);nYm=h.*X;nsubplot(3,2,2);plot(f/pi,abs(Ym)ntitle(|FTx(n)FTh(n)|);xlabel(/);ylabel(Ym幅度)naxis(0,0.5,0,max(abs(Ym);75显示结果76此处讲一下-离散系统时频域综合分析771 1、系统函数的零极点分析、系统函数的零极点分析离散时间系统的系统函数定义为系统零
29、状态响应的离散时间系统的系统函数定义为系统零状态响应的z z变换与变换与激励的激励的z z变换之比变换之比:)()()(zXzYzH如果系统函数如果系统函数)(zH的有理函数表示式为的有理函数表示式为11211121)(nnnnmmmmazazazabzbzbzbzH离散系统时频域综合分析离散系统时频域综合分析78在在MATLABMATLAB中系统函数的零极点就可通过函数中系统函数的零极点就可通过函数rootsroots得到，也可得到，也可借助借助DSPDSP工具箱中的函数工具箱中的函数tf2zptf2zp得到，得到，tf2zptf2zp的语句格式为：的语句格式为：R,P,K=tf2zp(B,
30、A)R,P,K=tf2zp(B,A)其中，其中，B B与与A A分别表示分子与分母多项式的系数向量。分别表示分子与分母多项式的系数向量。它的作用是将它的作用是将H(z)H(z)的有理分式表示式转换为零极点增益形式：的有理分式表示式转换为零极点增益形式：)()()()()(2121nmpzpzpzzzzzzzkzHMATLABMATLAB实现实现离散系统时频域综合分析离散系统时频域综合分析79离散系统时频域综合分析离散系统时频域综合分析【例【例1 1】已知一离散因果已知一离散因果LTILTI系统的系统函数为：系统的系统函数为：16.032.0)(2zzzzH试用试用MATLABMATLAB命令求
31、该系统的零极点。命令求该系统的零极点。802121216.0132.016.032.0)(zzzzzzzzHB=1,0.32;B=1,0.32;A=1,1,0.16;A=1,1,0.16;R,P,K=tf2zp(B,A)R,P,K=tf2zp(B,A)R=R=-0.3200 -0.3200P=P=-0.8000 -0.8000 -0.2000 -0.2000K=K=1 1因此，零点为因此，零点为32.0z，极点为：，极点为：8.01p2.02p。求解：81离散系统时频域综合分析若要获得系统函数的零极点分布图，可直接应用若要获得系统函数的零极点分布图，可直接应用zplanezplane函数，函数
32、，其语句格式为：其语句格式为：zplane(B,A)zplane(B,A)其中，其中，B B与与A A分别表示的分子和分母多项式的系数向量。分别表示的分子和分母多项式的系数向量。它的作用是在它的作用是在Z Z平面上画出单位圆、零点与极点。平面上画出单位圆、零点与极点。82离散系统时频域综合分析离散系统时频域综合分析【例【例2 2】已知一离散因果已知一离散因果LTILTI系统的系统函数为：系统的系统函数为：试用试用MATLABMATLAB命令绘出该系统的零极点分布图。命令绘出该系统的零极点分布图。68.052.136.0)(22zzzzH83求解：。2122268.052.1136.0168.0
33、52.136.0)(zzzzzzzHMATLABMATLAB源程序为：源程序为：B=1,0,-0.36;B=1,0,-0.36;A=1,-1.52,0.68;A=1,-1.52,0.68;zplane(B,A),grid onzplane(B,A),grid onlegend(legend(零点零点,极点极点)title(title(零极点分布图零极点分布图)程序运行结果如图所示。程序运行结果如图所示。84离散系统时频域综合分析离散系统时频域综合分析)(zH系统函数的零极点分布与其时域特性的关系系统函数的零极点分布与其时域特性的关系：在离散系统中，在离散系统中，z z变换建立了时域函数变换建立
34、了时域函数 与与z z域函数域函数之间的对应关系。因此，之间的对应关系。因此，z z变换的函数变换的函数 从形式可以从形式可以反映反映 的部分内在性质。的部分内在性质。我们通过讨论我们通过讨论H(z)H(z)的一阶极点情况，来说明系统函数的零极点的一阶极点情况，来说明系统函数的零极点分布与系统时域特性的关系。分布与系统时域特性的关系。)(zH)(nh)(nh85离散系统时频域综合分析离散系统时频域综合分析nMATLABMATLAB求解单位抽样响应求解单位抽样响应 可利用函数可利用函数filterfilter，filterfilter函数的常用语句格式为：函数的常用语句格式为：ny=filter
35、(b,a,x)y=filter(b,a,x)n表示由向量表示由向量b b和和a a组成的系统对输入组成的系统对输入x x进行滤波，系统的进行滤波，系统的输出为输出为y;y;)(nh2 2、系统时域响应分析、系统时域响应分析86离散系统时频域综合分析离散系统时频域综合分析nMATLABMATLAB另一种求单位抽样响应另一种求单位抽样响应 的方法是利用控制的方法是利用控制系统工具箱提供的函数系统工具箱提供的函数impzimpz来实现。来实现。impzimpz函数的常用函数的常用语句格式为语句格式为nimpz(b,a,N)impz(b,a,N)n其中，参数其中，参数N N通常为正整数，代表计算单位抽
36、样响应的通常为正整数，代表计算单位抽样响应的样值个数。样值个数。)(nh87离散系统时频域综合分析离散系统时频域综合分析【例【例1 1】试用试用MATLABMATLAB命令画出系统函数的零极点分布图、命令画出系统函数的零极点分布图、以及对应的时域单位抽样响应以及对应的时域单位抽样响应 的波形。的波形。)(nh18.0118.0)(ZzzzHb1=1,0;a1=1,-0.8;subplot(121)zplane(b1,a1)title(极点在单位圆内的正实数)subplot(122)impz(b1,a1,30);grid on;88结果显示：结果显示：89离散系统时频域综合分析离散系统时频域综合
37、分析3 3、离散时间、离散时间LTILTI系统的频率特性分析系统的频率特性分析 离散时间系统的频率响应定义为：离散时间系统的频率响应定义为：)(|)(|)(jjjeeHeH|)(|jeH)(其中，其中，称为离散时间系统的幅频特性；称为离散时间系统的相频特性。90离散系统时频域综合分析离散系统时频域综合分析MATLAB提供了求离散时间系统频响特性的函数freqz，调用freqz的格式主要有两种。一种形式为H,w=freqz(B,A,N)其中，B与A分别表示的分子和分母多项式的系数向量；N为正整数，默认值为512；返回值w包含 范围内的N个频率等分点；返回值H则是离散时间系统频率响应。091离散系
38、统时频域综合分析离散系统时频域综合分析另一种形式为：另一种形式为：H,w=freqz(B,A,N,whole)与第一种方式不同之处在于角频率的范围扩展到与第一种方式不同之处在于角频率的范围扩展到2092离散系统时频域综合分析离散系统时频域综合分析【例【例1 1】试用试用MATLABMATLAB命令绘制以下系统的频率响应曲线。命令绘制以下系统的频率响应曲线。8109.056.19028.096.0)(22zzzzzH解：利用函数解：利用函数freqz计算出计算出)(jeH然后利用函数然后利用函数abs和和angle分别求出幅频特性与相频特性，分别求出幅频特性与相频特性，最后利用最后利用plot命
39、令绘出曲线。命令绘出曲线。93离散系统时频域综合分析离散系统时频域综合分析nMATLAB源程序为：源程序为：nb=1-0.96 0.9028;na=1-1.56 0.8109;nH,w=freqz(b,a,400,whole);nHm=abs(H);nHp=angle(H);nsubplot(211)nplot(w,Hm),grid onnxlabel(omega(rad/s),ylabel(Magnitude)ntitle(离散系统幅频特性曲线离散系统幅频特性曲线)nsubplot(212)nplot(w,Hp),grid onnxlabel(omega(rad/s),ylabel(Phas
40、e)ntitle(离散系统相频特性曲线离散系统相频特性曲线)94结果显示：结果显示：95最后一个内容DFT变换967.3 离散傅里叶变换（DFT）n定义DFT：用类似于例7.9中的方法，可把（7.3）式写成矩阵乘法运算。其中，xn为序列行向量，Wnk是一NN阶方阵，而 称为旋转因子。10()DFT()(),0,1NknNnX kx nx n WkN()knnkXxW)1()1(1)1(0)1()1(11101)1(01000NNNNNNNNNNNNNNNWWWWWWWWWnkWNNW2je977.3 离散傅里叶变换（DFT）n用矩阵乘法计算N点DFT的程序如下。nn MATLAB程序q73a.
41、mn%用矩阵乘法计算N点DFTnclear;close allnxn=input(请输入序列x=);nN=length(xn);%nn=0:N-1;k=n;nk=n*k;%生成NN方阵nWN=exp(-j*2*pi/N);nWnk=WN.nk;%产生旋转因子矩阵nXk=xn*Wnk;%计算N点DFT98例7.15 序列的离散傅立叶变换n求复正弦序列 n余弦序列 n正弦序列的离散傅立叶变换，分别按N=16和N=8进行计算。绘出幅频特性曲线，进行比较讨论。)(e)(8j1nRnxNn)(cos)(2nRnnxN)(8sin)(3nRnnxN 99程序分析nclear;close allnN=16;
42、nN1=8;n%n%产生序列x1(n),计算DFTx1(n)nn=0:N-1;nx1n=exp(j*pi*n/8);%产生x1(n)nX1k=fft(x1n,N);%计算N点DFTx1(n)nXk1=fft(x1n,N1);%计算N1点DFTx1(n)n%n%产生序列x2(n),计算DFTx2(n)nx2n=cos(pi*n/8);nX2k=fft(x2n,N);%计算N点DFTx2(n)nXk2=fft(x2n,N1);%计算N1点DFTx1(n)n%n%产生序列x3(n),计算DFTx3(n)nx3n=sin(pi*n/8);nX3k=fft(x3n,N);%计算N点DFTx3(n)nXk
43、3=fft(x3n,N1);%计算N1点DFTx1(n)n%100例7.15 序列的离散傅立叶变换n在截取16点时，得到的是完整的余弦波形；而截取8点时，得到的是半截的余弦波形，当然有大量的谐波成分。101本题分析nDFT求解n结果会分析102例7.16 验证N点DFT的物理意义（1），绘出幅频曲线和相频曲线。（2）计算并图示x(n)的8点DFT。（3）计算并图示x(n)的16点DFT。解:序列x(n)的点DFT的物理意义是 在0，2上进行点等间隔采样。程序先密集采样，绘制出幅频曲线图。然后再分别做8点和16点DFT来验证这个采样关系。j4j4j1 e()(),(e)FT()1 ex nR n
44、Xx n求得)e(jX结果1：103结果2104结果3105106对本章对本章MATLAB在数字信号在数字信号处理中的应用的内容讲解，处理中的应用的内容讲解，告一段落。告一段落。107例7.10 时域采样频率与频谱混叠n分别以采样频率fs=1000Hz，400Hz和200Hz对xa(t)进行等间隔采样，计算并图示三种采样频率下的采样信号及其幅频特性n解：程序分别设定4种采样频率fs=10kHz，1kHz，400Hz和200Hz，对xa(t)进行采样，得到采样序列xa(t)，xa1(n)，xa2(n)，xa3(n)，画出其幅度频谱。采样时间区间均为0.1秒。为了便于比较，画出了幅度归一化的幅频曲
45、线，如图7.11所示。108例7.10 采样频率与频谱混叠(续)由于由以上关系式可见，采样信号的频谱函数是原模拟信号频谱函数的周期延拓，延拓周期为2/T。如果以频率f为自变量（=2f），则以采样频率fs=1/T为延拓周期。对频带限于fc的模拟信号xa(t)，只有当fs2fc时，采样后 才不会发生频谱混叠失真。这就是著名的采样定理 knTTTkXTnxnxXn2j1e)()(FT)e(ajj)e(j TX109例7.11 由离散序列恢复模拟信号n用时域内插公式其中模拟用理想低通滤波器恢复的过程，观察恢复波形，计算出最大恢复误差。解：这个公式与卷积公式相像，可以用向量和矩阵乘法来解决。nnTtgn
46、xtx)()()(as()sin()sinc()g tttF tTT110例7.11 由离散序列恢复模拟信号nxa=x*g(TNM)=x*G n其中 G=sinc(Fs*TNM)nM表示在两个采样点之间增加的间隔数，使输出更密，更接近模拟信号。1)T(Nt(M)1)T(Nt(2)1)T(Nt(1)2Tt(M)2Tt(2)2Tt(1)Tt(M)Tt(2)Tt(1)t(M)t(2)t(1)TNM111例7.16 验证N点DFT的物理意义（1），绘出幅频曲线和相频曲线。（2）计算并图示x(n)的8点DFT。（3）计算并图示x(n)的16点DFT。解:序列x(n)的点DFT的物理意义是 在0，2上进行
47、点等间隔采样。程序先密集采样，绘制出幅频曲线图。然后再分别做8点和16点DFT来验证这个采样关系。j4j4j1 e()(),(e)FT()1 ex nR nXx n求得)e(jX112例7.17 频域与时域采样对偶性 n（1）产生三角波序列n（2）对M=40，计算x(n)的64点DFT，并图示x(n)和X(k)=DFTx(n)，k=0,1,63。n（3）对（2）中所得X(k)在 0，2 上进行32点抽样得n（4）求的32点IDFT，即n n（5）绘出 的波形图，评述它与x(n)的关系。0 /2()/2 nnMx nMnMnM1()(2)k=0 1X kXk，31)(1kX11()IDFT()n
48、=0 1x nX k，31321)(nx113例7.17 频域与时域采样对偶性n由于频域在0，2上的采样点数N（N=32）小于x(n)的长度M（M=40），所以，产生时域混叠现象，不能由X1(k)恢复原序列x(n)。只有满足NM时，可由频域采样X1(k)得到原序列x(n)。n这就是频域采样定理。对NM的情况，请读者自己编程验证。)(IDFT)(1kXnx114例7.18 快速卷积快速卷积就是根据DFT的循环卷积性质，将时域卷积转换为频域相乘，最后再进行IDFT得到时域卷积序列y(n)。其中时域和频域之间的变换均用FFT实现，所以使卷积速度大大提高。框图如下:L点FETL点FETL点FET115
49、例7.19 用DFT求连续信号频谱n在计算机上用DFT对模拟信号进行谱分析时，只能以有限大的采样频率fs对模拟信号采样有限点样本序列（等价于截取模拟信号一段进行采样）作DFT变换，得到模拟信号的近似频谱。其误差主要来自以下因素：n 截断效应（频谱泄露和谱间干扰）n 频谱混叠失真n因素使谱分辨率（能分辨开的两根谱线间的最小间距）降低，并产生谱间干扰；n因素使折叠频率（fs/2）附近的频谱产生较大失真。116例7.19 用DFT求连续信号频谱n加大截取长度Tp可提高频率分辨率；选择合适的窗函数可降低谱间干扰；而频谱混叠失真要通过提高采样频率fs和（或）预滤波（在采样之前滤除折叠频率以外的频率成分）
50、来改善。n编写程序q719.m验证截断效应及加窗的改善作用，先选取以下参数：采样频率fs=400Hz，T=1/fs 采样信号序列 对x(n)作4096点DFT作为的近似频谱Xa(jf)。a()()(),n=0,1,N-1x nx nT w n117例7.19 用DFT求连续信号频谱n如图7.19所示。图中X1(jf)，X4(jf)和X8(jf)分别表示Tp=0.04s，0.04*4s和0.04*8s时的谱分析结果。由图可见，由于截断使原频谱中的单频谱线展宽（又称之为泄漏），Tp越大泄漏越小，频率分辨率越高。Tp=0.04s时，25Hz与50Hz两根谱线已分辨不清了。所以实际谱分析的截取时间Tp
