《理论课第7章yyxMATLAB在数字信号处理中的应用精讲课件.ppt》由会员分享，可在线阅读，更多相关《理论课第7章yyxMATLAB在数字信号处理中的应用精讲课件.ppt（156页珍藏版）》请在三一办公上搜索。
1、1,第 7章 MATLAB在数字信号处理中的应用,2,书山有路勤为径，学海无涯苦作舟。,好好学习，天天向上,3,引言,数字信号处理中包含离散时间信号与系统的一部分知识。本章主要研究内容分为两部分：离散时间信号与系统时域分析离散时间信号与系统频域分析在数字信号处理中还有一个非常大的知识点：数字滤波器设计，留待下学期开设本课程时，可以学习。,4,一、离散时间信号与系统时域分析,5,包含内容,信号分析典型信号表示impseq(),stepseq()7.1 序列相加，相乘7.2 序列合成和截取7.3 序列移位和周期延拓7.6离散序列卷积系统分析7.4系统响应7.5系统线性判定,6,7.1 离散信号的产
2、生及时域处理,时域离散信号用x(n)表示（离散序列）。x和n同时使用才能完整地表示一个序列：时间变量n（表示采样位置）只能取整数；向量x（表示采样点的幅度）与n一一对应。由于n序列是按整数递增的,可简单地用其初值ns决定,因为它的终值nf取决于ns 和x的长度length(x),故可写成：n=ns:nf 或n=ns:nslength(x)1,7,典型信号表示impseq(),stepseq(),8,例7.1 序列的相加和相乘,给出两个序列x1(n)和x2(n)。x1=0,1,2,3,4,3,2,1,0;n1=-2:6;x2=2,2,0,0,0,-2,-2;n2=2:8;要求它们的和ya及乘积y
3、p。解：编程的思路是把序列长度延拓到覆盖n1和n2的范围，这样才能把两序列的时间变量对应起来，然后进行对应元素的运算。,9,例程,x1=0,1,2,3,4,3,2,1,0;ns1=-2;%给定x1及ns1x2=2,2,0,0,0,-2,-2;ns2=2;%给定x2及ns2nf1=ns1+length(x1)-1;nf2=ns2+length(x2)-1;ny=min(ns1,ns2):max(nf1,nf2);%y(n)的时间变量xa1=zeros(1,length(ny);xa2=xa1;%延拓序列初始化xa1(find(ny=ns1)%给xa2 赋值x2ya=xa1+xa2%序列相加yp=
4、xa1.*xa2,10,结果显示,11,例7.2 序列的合成和截取,用例6.13的结果编写产生矩形序列RN(n)的程序。序列起点为n0，矩形序列起点为n1，长度为N（n0，n1，N由键盘输入）。并用它截取一个复正弦序列exp(jn/8).解：建模:矩形序列可看成两个阶跃序列之差。用MATLAB逻辑关系产生矩形序列x2(n)。而用它截取任何序列相当于元素群相乘x2.*x，也称为加窗运算。序列的合成和截取就是相加和相乘。,12,结果显示,13,例7.2 序列的合成和截取,用例6.13的结果编写产生矩形序列RN(n)的程序。序列起点为n0，矩形序列起点为n1，长度为N（n0，n1，N由键盘输入）。并
5、用它截取一个复正弦序列exp(jn/8).解：建模:矩形序列可看成两个阶跃序列之差。用MATLAB逻辑关系产生矩形序列x2(n)。而用它截取任何序列相当于元素群相乘x2.*x，也称为加窗运算。序列的合成和截取就是相加和相乘。,14,程序要点,n=n0:n1+N+5;%生成自变量数组u=(n-n1)=0;%产生单位阶跃序列（u(n-n1)）x1=(n-n1)=0-(n-n1-N)=0%用阶跃序列之差产生矩形序列x2=(n=n1)%对复正弦序列加矩形窗（元素群乘）,15,程序显示结果,16,例7.3 序列的移位和周期延拓,已知，利用MATLAB生成并图示 表示x(n)以8为周期的延拓）和解：方法1
6、，利用矩阵乘法和冒号运算 x=1 2 3 4;y=x*ones(1，3);方法2，采用求余函数mod，y=x(mod(n,M)+1)可实现对x(n)以M为周期的周期延拓。加1是因为MATLAB向量下标只能从1开始，,17,程序实现分析,构造x(n)x1=(0.8).n;x2=(n=0)end,18,周期化实现,方法1：xc=xn(mod(n,8)+1);%产生 x(n)的周期延拓 xcm=xn(mod(n-m,8)+1);%产生 x(n)移位后的周期延拓,19,周期化实现,方法2x=1 2 3 4;y=x*ones(1,3);y=1 1 1 2 2 2 3 3 3 4 4 4y1=y(:);%
7、1 2 3 4 1 2 3 4 1 2 3 4,20,程序结果,21,例7.6 离散序列的卷积计算,给出两个序列 和,计算其卷积y(n)，并图示各输入输出序列。解：在例6.4中，已经给出了直接调用MATLAB的卷积函数conv的方法，也给出了自编卷积计算程序的方法，要注意的是本例时间变量的设定和移位方法。在本例中，设定n为从零开始，向量x和h的长度分别为Nx=20和Nh=10；结果向量y的长度为length(y)=Nx+Nh-1。,22,程序要点,x1=(0.9).n;%x1(n)%产生x2(n)=x1(n-m)x2=zeros(1,Nx+m);for k=m+1:m+Nx x2(k)=x1(
8、k-m);endnh=0:Nh-1;h1=ones(1,Nh);%h1(n)h2=h1;%h2(n)y1=conv(x1,h1);%y1(n)y2=conv(x2,h2);%y2(n)Ny1=length(y1);%y1(n)长度Ny2=length(y2);%y2(n)长度,23,程序显示结果,24,例7.4 离散系统对信号的响应,本题给定6阶低通数字滤波器的系统函数，求它在下列输入序列x(n)下的输出序列y(n)。解：本题的计算原理见例6.14，在这里用工具箱函数filter来解。如果已知系统函数H(z)=B(z)/A(z)，则filter函数可求出系统对输入信号x(n)的响应y(n)。y
9、=filter(B,A,x)由差分方程可得到H(z)的分子和分母多项式系数向量A和B,再给出输入向量x即可。,系统分析,25,程序要点,%设定系统参数A,BB=0.0003738*conv(1,1,conv(1,1,conv(1,1,conv(1,1,conv(1,1,1,1)A=conv(1,-1.2686,0.7051,conv(1,-1.0106,0.3583,1,-0.9044,0.2155)x1=n=0;%产生输入信号x1(n)y1=filter(B,A,x1);%对x1(n)的响应x2=(n-m)=0;%产生输入信号x2(n)y2=filter(B,A,x2);%对x2(n)的响应
10、x3=n=0;%产生输入信号x3(n)y3=filter(B,A,x3);%对x3(n)的响应x4=(n=0)%对x5(n)的响应,26,显示结果,27,例7.5 系统线性性质验证,设系统差分方程为y(n)=x(n)+0.8y(n-1)要求用程序验证系统的线性性质。解：产生两种输入序列，分别乘以常数后1。分别激励系统，再求输出之和；2。先相加，再激励系统求输出；对两个结果进行比较，方法是求它们之差，按误差的绝对值是否极小进行判断。,28,程序要点,B=1;A=1,-0.8;x1=0.8.n;%产生输入信号x1(n)x=(n=0)%y(n)=5y1(n)+3y2(n),29,程序显示结果,30,
11、二、离散时间信号与系统频域分析,31,包含内容,Z变换序列z变换与逆变换7.7Z反变换法求解系统响应7.8离散时间傅里叶变换7.9离散傅里叶变换7.15系统频域响应7.12,7.13,32,(1)序列的正反Z变换,其中，符号表示取z变换，z是复变量。相应地，单边z变换定义为：,Z变换与其逆变换,33,MATLAB符号数学工具箱提供了计算离散时间信号单边z变换的函数ztrans和z反变换函数iztrans，其语句格式分别为,Z=ztrans(x)x=iztrans(z),上式中的x和Z分别为时域表达式和z域表达式的符号表示，可通过sym函数来定义。,MATLAB实现-方法1（符号数学工具箱）,3
12、4,【例1】试用ztrans函数求下列函数的z变换。,x=sym(an*cos(pi*n);Z=ztrans(x);simplify(Z)ans=z/(z+a),35,【例2】试用iztrans函数求下列函数的z反变换。,Z=sym(8*z-19)/(z2-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)是函数在MATLAB符号工具箱中的表示，反变换后的函数形式为：,36,例7.14 用符号运算工具箱解z变换,解
13、：无限长度时间序列的z变换和逆z变换都属于符号运算的范围。MATLAB的symbolic（符号运算）工具箱已提供了这种函数。如果读者已在计算机上安装了这个工具箱，可以键入以下程序。MATLAB程序q714.m 其特点是程序的开始要指定符号自变量syms z n a N w0%规定z,n,a为符号变量,37,程序分析：,syms z n aN w0%规定z,n,a为符号变量y1=an;%给出y的表示式Y1=ztrans(y1)%求y的z变换X1=-3*z-1/(2-5*z-1+2*z-2);%给出X的表示式x1=iztrans(X1)%求X的逆Z变换simplify(x);%可显示变换后表达式的
14、形式,38,如果信号的z域表示式是有理函数，进行z反变换的另一个方法是对进行部分分式展开，然后求各简单分式的z反变换.如果X(z)的有理分式表示为：,MATLAB实现-方法2（residuez函数实现反变换）,39,MATLAB信号处理工具箱提供了一个对进行部分分式展开的函数residuez，其语句格式为：R,P,K=residuez(B,A)其中，B，A分别表示X(z)的分子与分母多项式的系数向量，分子与分母多项式按照z-1升幂排列。R为部分分式的系数向量；P为极点向量；K为多项式的系数。若X(z)为有理真分式，则K为零。,40,【例1】试用MATLAB命令进行部分分式展开，并求出其z反变换
15、。,解：MATLAB源程序为B=18;A=18,3,-4,-1;R,P,K=residuez(B,A)R=0.3600 0.2400 0.4000,41,P=0.5000-0.3333-0.3333K=,从运行结果可知，,表示系统有一个二重极点。所以，X(z)的部分分式展开为,42,教材中求z的逆变换的方法,对于z变换分式可以用部分分式法或长除法求其反变换。用函数residuez可以求出它的极点留数分解其中r,p,k=residuez(B,A)其反变换为：,43,注意几点：,r,p,k=residuez(B,A),A首项不能为零，即分母z最高次项不能为0k只有当分式为假分式时会有值，即对应冲激
16、序列部分（FIR）；真分式，值为零还有一些细节需要处理：,44,例7.7 有限序列的z和逆z变换,两序列x1=1,2,3，n1=-1:1 及x2=2,4,3,5，n2=-2:1，求出x1与x2及其卷积x的z变换。解：其z变换可写成两个多项式乘积可用conv函数来求得。n数组要自己判别。n的起点ns=ns1+ns2=3，终点nf=nf1+nf2=2。n=ns:nf。由x和n即可得出X(z)。,45,程序分析,x1=1,2,3;ns1=-1;%设定x1和ns1nf1=ns1+length(x1)-1;%nf1可以算出x2=2,4,3,5;ns2=-2;%设定x2和ns2nf2=ns1+length
17、(x1)-1;%nf2可以算出x=conv(x1,x2)%求出xn=(ns1+ns2):(nf1+nf2)%求出n,46,47,例7.8 求z多项式分式的逆变换,设系统函数为，输入例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变换，得,48,例7.8 z多项式分式逆变换(续),由程序算出nsy=-1，留数、极点分别为r=-57.7581 和 204.7581p=0.7791 和 0.3209k=-150-30代入得,49,Z
18、反变换法求解系统输出(编程实现),Y(z)=X(z)W(z)x=2,4,3,5;nsx=-2;%输入序列及初始时间nfx=nsx+length(x)-1;%计算序列终止时间Bw=-3;nsbw=-1;%系统函数的分子系数,及z的最高次数Aw=2,-2.2,0.5;nsaw=0;%系统函数的分母系数,及z的最高次数B=conv(-3,x);%输入与分子z变换的多项式乘积A=Aw;%分母不变y(n)=IZTY(Z)r,p,k=residuez(B,A)%求留数r,极点p及直接项knf=input(终点时间nf=);%要求键入终点时间n=min(nsx,nsy):nf;%生成总时间数组%求无限序列y
19、i和直接序列ydyi=(r(1)*p(1).(n-nsy)+r(2)*p(2).(n-nsy).*stepseq(nsy,n(1),nf);yd=k(1)*impseq(nsy,n(1),nf)+k(2)*impseq(-1-nsy,n(1),nf);y=yi+yd;,50,另外一种方法(filter函数)：,filter(B,A,x),51,例7.9 离散时间傅里叶变换,取周期的正弦信号，作8点采样，求它的连续频谱。然后对该信号进行N个周期延拓，再求它的连续频谱。把N无限增大，比较分析其结果。解：先求离散傅立叶变换的MATLAB子程序最后得到X=x*exp(-j*w*n)。有了子程序，本例就
20、没有什么难度了。,52,DTFT程序：,function X=dtft(x,w)%计算离散时间傅立叶变换%X=dtft(x,n,w),%X=在w频率点上的DTFT数组,%x=沿n的有限长度序列,%n=样本位置向量%w=频率点位置向量n=1:length(x);ewn=exp(-n*w*i);X=x*ewn;,53,正弦信号x0的DTFT变换,x0=sin(2*pi*1:8/8)*5;%x0是8点行向量dt=2*pi/8;w=linspace(-2*pi,2*pi,1000)/dt;%w是1000点行向量X0=dtft(x0,w)*dt;%求得频率响应X0,54,正弦信号延拓N个周期后，DTFT
21、,N=4;%延拓周期数x1=reshape(x0*ones(1,N),1,N*length(x0);%延拓后的时域信号x1X1=dtft(x1,w)*dt;%求x1的频率响应X1%延拓100次后%x1=reshape(x0*ones(1,100),1,100*length(x0);%延拓后的时域信号x1%X1=dtft(x1,w)/100*dt;%求x1的频率响应X1,55,重复无穷次,disp(重复无穷次的八点信号的离散时间傅立叶变换-傅立叶级数)pause,X2=fft(x0*dt);%离散傅立叶变换w1=2*pi*0:length(x0)-1/length(x0);%离散频点向量subp
22、lot(3,1,3),stem(-w1,w1,abs(X2),abs(X2),grid,axis(min(w),max(w),0,max(abs(X2),grid,56,例7.9 离散时间傅里叶变换2,程序运行结果执行程序q709并按提示键入N=4，所得图形如图7.10所示。N取得愈大，其峰值愈大，宽度愈窄。当N取得很大时，会出现内存不足的问题，这是用矩阵乘法做傅里叶变换的缺点。另外，因为那时峰值点处的宽度很窄，也会出现所选频点对不上峰值点的问题。所以对于N无限增大的情况，必须用fft函数来求。这时用连续频谱也没有意义了。这里用同样的横坐标把几种频谱进行对比，使读者更好地理解其关系。,57,程
23、序显示结果：,58,例7.10 时域采样频率与频谱混叠,分别以采样频率fs=1000Hz，400Hz和200Hz对xa(t)进行等间隔采样，计算并图示三种采样频率下的采样信号及其幅频特性解：程序分别设定4种采样频率fs=10kHz，1kHz，400Hz和200Hz，对xa(t)进行采样，得到采样序列xa(t)，xa1(n)，xa2(n)，xa3(n)，画出其幅度频谱。采样时间区间均为0.1秒。为了便于比较，画出了幅度归一化的幅频曲线，如图7.11所示。,59,例7.10 采样频率与频谱混叠(续),由于由以上关系式可见，采样信号的频谱函数是原模拟信号频谱函数的周期延拓，延拓周期为2/T。如果以频
24、率f为自变量（=2f），则以采样频率fs=1/T为延拓周期。对频带限于fc的模拟信号xa(t)，只有当fs2fc时，采样后 才不会发生频谱混叠失真。这就是著名的采样定理,程序要点：,t=0:1/fs:0.1;%fs代不同的取样频率 xa=exp(-a*t).*sin(b*t);k=0:511;f=fs*k/512;%由wk=2k/512=2fT求得模拟频率fXa=dtft(xa,2*pi*k/512);%近似模拟信号频谱,60,结果1,61,结果2,62,63,例7.12 梳状滤波器零极点和幅特性,梳状滤波器系统函数有如下两种类型。FIR型：IIR型：freqz 数字滤波器频率特性计算和绘制函
25、数 zplane H(z)的零-极点图绘制。解：调用函数freqz和zplane 很容易写出程序q712.m。,64,程序分析,freqzfreqz(b,a,w)zplane,65,频率响应,b=1,0,0,0,0,0,0,0,-1;a0=1;a1=1,0,0,0,0,0,0,0,-(0.8)8;a2=1,0,0,0,0,0,0,0,-(0.9)8;a3=1,0,0,0,0,0,0,0,-(0.98)8;H,w=freqz(b,a0);H1,w1=freqz(b,a1);H2,w2=freqz(b,a2);H3,w3=freqz(b,a3);,66,零极点及绘图,subplot(2,2,1);
26、zplane(b,a0);title(FIR梳状滤波器零点图);subplot(2,2,2);zplane(b,a1);title(IIR梳状滤波器零、极点图,a=0.8);subplot(2,2,3);plot(w/pi,abs(H);title(FIR梳状滤波器幅频响应曲线);ylabel(幅度);xlabel(/);subplot(2,2,4);plot(w1/pi,abs(H1);title(IIR梳状滤波器幅频响应曲线,a=0.8);ylabel(幅度);xlabel(/);figure(2);subplot(2,2,1);zplane(b,a2);title(IIR梳状滤波器零、极
27、点图,a=0.9);subplot(2,2,2);zplane(b,a3);title(IIR梳状滤波器零、极点图,a=0.98);subplot(2,2,3);plot(w2/pi,abs(H2);title(IIR梳状滤波器幅频响应曲线,a=0.9);ylabel(幅度);xlabel(/);subplot(2,2,4);plot(w3/pi,abs(H3);title(IIR梳状滤波器幅频响应曲线,a=0.98);ylabel(幅度);xlabel(/),67,结果显示,68,结果显示,69,例7.13 低通滤波及时域卷积定理,输入信号 x(n)=cos(0.04n)+cos(0.08n
28、)+cos(0.4n)+0.3(n)，0n63 通过低通滤波器，计算滤波器对x(n)的响应输出y(n)，并图示x(n)和y(n)，观察滤波效果。解：如前所述，只要求出H(z)=B(z)/A(z)的分子和分母多项式系数向量B和A，则可调用滤波器直接型实现函数filter对输入信号x(n)进行滤波。y=filter(B,A,x),70,程序要点,输入信号构造%产生输入信号x(n)n=0:255;N=4096;x=cos(0.04*pi*n)+cos(0.08*pi*n)+cos(0.4*pi*n);w=randn(size(x);%产生正态零均值噪声x=x+0.3*w;,71,程序要点,%求H(z
29、)分子分母多项式系数向量B和Ab=1,2,1;%(1+z-1)2 的展开系数B=0.0003738*conv(conv(b,b),b);%嵌套调用卷积函数conva1=1,-1.2686,0.7051;a2=1,-1.0106,0.3583;a3=1,-0.9044,0.2155;A=conv(conv(a1,a2),a3);,72,程序要点,滤波及图形%对x(n)滤波y=filter(B,A,x);%绘图X=fft(x,N);Y=fft(y,N);subplot(3,2,1);stem(x,.)axis(0,max(n)/4,min(x),max(x);line(0,max(n),0,0)t
30、itle(输入信号x(n);xlabel(n);ylabel(x(n)subplot(3,2,5);stem(y,.)axis(0,max(n)/4,min(y),max(y);line(0,max(n),0,0)title(输出信号 y(n);xlabel(n);ylabel(y(n)k=0:N-1;f=2*k/N;subplot(3,2,2);plot(f,abs(X)title(输入信号x(n)的幅频曲线);xlabel(/);ylabel(|FTx(n)|)axis(0,0.5,0,max(abs(X);subplot(3,2,6);plot(f,abs(Y)title(输出信号y(n
31、)的幅频曲线);xlabel(/);ylabel(|FTy(n)|)axis(0,0.5,0,max(abs(Y);,73,结果显示：,74,程序要点,h,f=freqz(B,A,N,whole);figure(2)subplot(3,2,1);plot(f/pi,abs(h)title(滤波器幅频响应曲线);xlabel(/);ylabel(H幅度)axis(0,0.5,0,max(abs(h);Ym=h.*X;subplot(3,2,2);plot(f/pi,abs(Ym)title(|FTx(n)FTh(n)|);xlabel(/);ylabel(Ym幅度)axis(0,0.5,0,ma
32、x(abs(Ym);,75,显示结果,76,此处讲一下-离散系统时频域综合分析,77,1、系统函数的零极点分析,离散时间系统的系统函数定义为系统零状态响应的z变换与激励的z变换之比:,如果系统函数,的有理函数表示式为,离散系统时频域综合分析,78,在MATLAB中系统函数的零极点就可通过函数roots得到，也可借助DSP工具箱中的函数tf2zp得到，tf2zp的语句格式为：R,P,K=tf2zp(B,A)其中，B与A分别表示分子与分母多项式的系数向量。它的作用是将H(z)的有理分式表示式转换为零极点增益形式：,MATLAB实现,离散系统时频域综合分析,79,离散系统时频域综合分析,【例1】已知
33、一离散因果LTI系统的系统函数为：,试用MATLAB命令求该系统的零极点。,。,80,B=1,0.32;A=1,1,0.16;R,P,K=tf2zp(B,A)R=-0.3200P=-0.8000-0.2000K=1因此，零点为,，极点为：,。,求解：,81,离散系统时频域综合分析,若要获得系统函数的零极点分布图，可直接应用zplane函数，其语句格式为：zplane(B,A)其中，B与A分别表示的分子和分母多项式的系数向量。它的作用是在Z平面上画出单位圆、零点与极点。,。,82,离散系统时频域综合分析,【例2】已知一离散因果LTI系统的系统函数为：,试用MATLAB命令绘出该系统的零极点分布图
34、。,。,83,求解：,。,MATLAB源程序为：B=1,0,-0.36;A=1,-1.52,0.68;zplane(B,A),grid onlegend(零点,极点)title(零极点分布图)程序运行结果如图所示。,84,离散系统时频域综合分析,系统函数的零极点分布与其时域特性的关系：,在离散系统中，z变换建立了时域函数 与z域函数之间的对应关系。因此，z变换的函数 从形式可以反映 的部分内在性质。我们通过讨论H(z)的一阶极点情况，来说明系统函数的零极点分布与系统时域特性的关系。,85,离散系统时频域综合分析,MATLAB求解单位抽样响应 可利用函数filter，filter函数的常用语句格
35、式为：y=filter(b,a,x)表示由向量b和a组成的系统对输入x进行滤波，系统的输出为y;,2、系统时域响应分析,86,离散系统时频域综合分析,MATLAB另一种求单位抽样响应 的方法是利用控制系统工具箱提供的函数impz来实现。impz函数的常用语句格式为impz(b,a,N)其中，参数N通常为正整数，代表计算单位抽样响应的样值个数。,87,离散系统时频域综合分析,【例1】试用MATLAB命令画出系统函数的零极点分布图、以及对应的时域单位抽样响应 的波形。,b1=1,0;a1=1,-0.8;subplot(121)zplane(b1,a1)title(极点在单位圆内的正实数)subpl
36、ot(122)impz(b1,a1,30);grid on;,88,结果显示：,89,离散系统时频域综合分析,3、离散时间LTI系统的频率特性分析,离散时间系统的频率响应定义为：,其中，,称为离散时间系统的幅频特性；,称为离散时间系统的相频特性。,90,离散系统时频域综合分析,MATLAB提供了求离散时间系统频响特性的函数freqz，调用freqz的格式主要有两种。一种形式为H,w=freqz(B,A,N)其中，B与A分别表示的分子和分母多项式的系数向量；N为正整数，默认值为512；,返回值w包含 范围内的N个频率等分点；返回值H则是离散时间系统频率响应。,91,离散系统时频域综合分析,另一种
37、形式为：H,w=freqz(B,A,N,whole)与第一种方式不同之处在于角频率的范围扩展到,92,离散系统时频域综合分析,【例1】试用MATLAB命令绘制以下系统的频率响应曲线。,解：利用函数freqz计算出,然后利用函数abs和angle分别求出幅频特性与相频特性，最后利用plot命令绘出曲线。,93,离散系统时频域综合分析,MATLAB源程序为：b=1-0.96 0.9028;a=1-1.56 0.8109;H,w=freqz(b,a,400,whole);Hm=abs(H);Hp=angle(H);subplot(211)plot(w,Hm),grid onxlabel(omega(
38、rad/s),ylabel(Magnitude)title(离散系统幅频特性曲线)subplot(212)plot(w,Hp),grid onxlabel(omega(rad/s),ylabel(Phase)title(离散系统相频特性曲线),94,结果显示：,95,最后一个内容,DFT变换,96,7.3 离散傅里叶变换（DFT）,定义DFT：用类似于例7.9中的方法，可把（7.3）式写成矩阵乘法运算。其中，xn为序列行向量，Wnk是一NN阶方阵，而 称为旋转因子。,97,7.3 离散傅里叶变换（DFT）,用矩阵乘法计算N点DFT的程序如下。MATLAB程序q73a.m%用矩阵乘法计算N点DF
39、Tclear;close 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,98,例7.15 序列的离散傅立叶变换,求复正弦序列 余弦序列 正弦序列的离散傅立叶变换，分别按N=16和N=8进行计算。绘出幅频特性曲线，进行比较讨论。,99,程序分析,clear;close allN=16;N1=8;%产生序列x1(n),计算DFTx1(n)n=0:N-1;x1n=exp(j*pi*n/8);%产生x1(n)X1k=
40、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,N);%计算N点DFTx3(n)Xk3=fft(x3n,N1);%计算N1点DFTx1(n)%,100,例7.15 序列的离散傅立叶变换,在截取16点时，得到的是完整的余弦波形；而截取8点时，得到的是半截
41、的余弦波形，当然有大量的谐波成分。,101,本题分析,DFT求解结果会分析,102,例7.16 验证N点DFT的物理意义,（1），绘出幅频曲线和相频曲线。（2）计算并图示x(n)的8点DFT。（3）计算并图示x(n)的16点DFT。解:序列x(n)的点DFT的物理意义是 在0，2上进行点等间隔采样。程序先密集采样，绘制出幅频曲线图。然后再分别做8点和16点DFT来验证这个采样关系。,结果1：,103,结果2,104,结果3,105,106,对本章MATLAB在数字信号处理中的应用的内容讲解，告一段落。,107,例7.10 时域采样频率与频谱混叠,分别以采样频率fs=1000Hz，400Hz和2
42、00Hz对xa(t)进行等间隔采样，计算并图示三种采样频率下的采样信号及其幅频特性解：程序分别设定4种采样频率fs=10kHz，1kHz，400Hz和200Hz，对xa(t)进行采样，得到采样序列xa(t)，xa1(n)，xa2(n)，xa3(n)，画出其幅度频谱。采样时间区间均为0.1秒。为了便于比较，画出了幅度归一化的幅频曲线，如图7.11所示。,108,例7.10 采样频率与频谱混叠(续),由于由以上关系式可见，采样信号的频谱函数是原模拟信号频谱函数的周期延拓，延拓周期为2/T。如果以频率f为自变量（=2f），则以采样频率fs=1/T为延拓周期。对频带限于fc的模拟信号xa(t)，只有当
43、fs2fc时，采样后 才不会发生频谱混叠失真。这就是著名的采样定理,109,例7.11 由离散序列恢复模拟信号,用时域内插公式其中模拟用理想低通滤波器恢复的过程，观察恢复波形，计算出最大恢复误差。解：这个公式与卷积公式相像，可以用向量和矩阵乘法来解决。,110,例7.11 由离散序列恢复模拟信号,xa=x*g(TNM)=x*G 其中 G=sinc(Fs*TNM)M表示在两个采样点之间增加的间隔数，使输出更密，更接近模拟信号。,111,例7.16 验证N点DFT的物理意义,（1），绘出幅频曲线和相频曲线。（2）计算并图示x(n)的8点DFT。（3）计算并图示x(n)的16点DFT。解:序列x(n
44、)的点DFT的物理意义是 在0，2上进行点等间隔采样。程序先密集采样，绘制出幅频曲线图。然后再分别做8点和16点DFT来验证这个采样关系。,112,例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)的关系。,113,例7.17 频域与时域采样对偶性,由于频域在0，2上的采样点数N（N=32）小于x(n)的长度M（M=40），所以，产生时域混叠现象，不能由X1(k)
45、恢复原序列x(n)。只有满足NM时，可由频域采样X1(k)得到原序列x(n)。这就是频域采样定理。对NM的情况，请读者自己编程验证。,114,例7.18 快速卷积,快速卷积就是根据DFT的循环卷积性质，将时域卷积转换为频域相乘，最后再进行IDFT得到时域卷积序列y(n)。其中时域和频域之间的变换均用FFT实现，所以使卷积速度大大提高。框图如下:,115,例7.19 用DFT求连续信号频谱,在计算机上用DFT对模拟信号进行谱分析时，只能以有限大的采样频率fs对模拟信号采样有限点样本序列（等价于截取模拟信号一段进行采样）作DFT变换，得到模拟信号的近似频谱。其误差主要来自以下因素：截断效应（频谱泄
46、露和谱间干扰）频谱混叠失真因素使谱分辨率（能分辨开的两根谱线间的最小间距）降低，并产生谱间干扰；因素使折叠频率（fs/2）附近的频谱产生较大失真。,116,例7.19 用DFT求连续信号频谱,加大截取长度Tp可提高频率分辨率；选择合适的窗函数可降低谱间干扰；而频谱混叠失真要通过提高采样频率fs和（或）预滤波（在采样之前滤除折叠频率以外的频率成分）来改善。编写程序q719.m验证截断效应及加窗的改善作用，先选取以下参数：采样频率fs=400Hz，T=1/fs 采样信号序列 对x(n)作4096点DFT作为的近似频谱Xa(jf)。,117,例7.19 用DFT求连续信号频谱,如图7.19所示。图中
47、X1(jf)，X4(jf)和X8(jf)分别表示Tp=0.04s，0.04*4s和0.04*8s时的谱分析结果。由图可见，由于截断使原频谱中的单频谱线展宽（又称之为泄漏），Tp越大泄漏越小，频率分辨率越高。Tp=0.04s时，25Hz与50Hz两根谱线已分辨不清了。所以实际谱分析的截取时间Tp是由频率分辨率决定的。另外，在本应为零的频段上出现了一些参差不齐的小谱包（称为谱间干扰）。谱间干扰的大小取决于加窗的类型。用矩形窗比用Hamming窗的频率分辨率高（泄漏小），但谱间干扰刚好相反。,118,例7.20 IIR滤波器直接型的转换,程序调用了信号处理工具箱函数tf2sos和扩展函数dir2pa
48、r，sos,g=tf2sos(B,A)实现从直接型到级联型（二阶分割形式）的转换。g为式中的增益，sos为L6阶矩阵，表示式中的系数。,119,例7.20 IIR滤波器直接型的转换,Cp,Bp,Ap=dir2par(B,A)实现从直接型到并联型的转换。B为直接型H(z)的分子多项式系数向量，A为直接型H(z)的分母多项式系数向量；Cp，Bp，Ap的含义与扩展函数dir2par中的C，B，A相同。dir2par中又调用了复共轭对比较函数cplxcomp。由于dir2par和cplxcomp是文献7中开发的，不属于MATLAB工具箱函数，所以将其M文件清单附在程序q720.m之后。,120,例7.
49、20 IIR滤波器直接型的转换,根据计算结果,级联型H(z)表达式：级联型结构图。,121,例7.20 IIR滤波器直接型的转换,并联型结构H(z)表达式并联型结构图,122,例7.21 直接型结构到格型梯形结构,tf2latc函数实现直接型到格型转换 K,C=tf2latc(B,A)求出零-极点IIR系统格型梯形结构的格型参数向量K和梯形参数向量C（用A(1)归一化）。注意，当系统函数在单位圆上有极点时发生错误。K=tf2latc(1,A)求出全极点IIR系统的格型结构参数向量K。如果使用格式K,C=tf2latc(1,A)，返回的系数C为标量。K=tf2latc(B)求出FIR系统的格型梯
50、形结构参数（反射系数）向量K（用H(z)的常数项B(1)归一化）。,123,例7.21 直接型结构到格型梯形结构,直接型系统函数转换为格型梯形结构,124,例7.22 FIR滤波器直接型到其他型,系统函数为调用信号处理工具箱函数tf2sos和tf2latc,给变元A赋值1,B=2,13/12,5/4,2/3即可.级联型结构系数sos=1.0000 0.5360 0 1.0000 0 0 1.0000 0.0057 0.6219 1.0000 0 0g=2格型结构系数（反射系数）：K=0.2500 0.5000 0.3333,125,例7.22 FIR滤波器直接型到其他型,得出级联型为 格型结构
