资源描述
单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第 7章 MATLAB在数字信号处理中的应用,1,7.1 时域离散信号的产生及时域处理,时域离散信号用,x,(,n,)表示，时间变量,n,（表示采样位置）只能取整数。因此，,x,(,n,)是一个离散序列，以后简称序列。,用一个向量,x,不足以表示序列值,x,(,n,)。必须再用另一个等长的定位时间变量,n,。,x,和,n,同时使用才能完整地表示一个序列,，由于n序列是按整数递增的，可简单地用其初值ns决定，因为它的终值nf取决于ns 和x的长度length(x)，故可写成：,n = ns:nf,或n = ns: ns,length(x),1,2,单位脉冲序列,(,n,-,n,0,),的生成函数impseq,function x,n=impseq(n0,ns,nf),n=ns:nf;x=(n-n0)=0;,单位阶跃序列,u,(,n,-,n,0,),的生成函数stepseq,function x,n=stepseq(n0,ns,nf),n=ns:nf; x=(n-n0)=0;,P173表7.1,3,例7.1 序列的相加和相乘,给出两个序列,x,1,(,n,)和,x,2,(,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及乘积yp。,解：编程的思路是把序列长度延拓到覆盖n1和n2的范围，这样才能把两序列的时间变量对应起来，然后进行对应元素的运算。,4,%第七章例7.1q701离散信号序列的相加和相乘,x1=0,1,2,3,4,3,2,1,0;ns1=-2;% 给定x1及ns1,x2=2,2,0,0,0,-2,-2; ns2=2;% 给定x2及ns2,nf1=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)&(ny=ns2)&(ny=0; % 产生单位阶跃序列（u(n-n1)）,x1=(n-n1)=0-(n-n1-N)=0;% 用阶跃序列差产生矩形序列,x2=(n=n1)% 用逻辑式产生矩形序列,x3=exp(j*n*pi/8).*x2; % 对复正弦序列加矩形窗（元素群乘）,subplot(2,2,1);stem(n,x1,.);,xlabel(n);ylabel(x1(n);%,标注,axis(n0,max(n),0,1);%,定坐标范围,subplot(2,2,3);stem(n,x2,.);,xlabel(n);ylabel(x2(n);%,标注,axis(n0,max(n),0,1);%,定坐标范围,8,subplot(2,2,2);stem(n,real(x3),.);,xlabel(n);ylabel(x3(n),的实部,);%,标注,line(n0,max(n),0,0);%,画横轴,axis(n0,max(n),-1,1);%,定坐标范围,subplot(2,2,4);stem(n,imag(x3),.);,xlabel(n);ylabel(x3(n)的虚部);% 标注,line(n0,max(n),0,0);% 画横轴,axis(n0,max(n),-1,1);%定坐标范围,结果,输入序列起点:n0=-6,输入序列长度:N=15,输入位移:n1=3,9,10,例7.3 序列的移位和周期延拓,已知，利用MATLAB生成并图示 表示,x,(,n,)以8为周期的延拓）和 ，其中N=24。,解：周期延拓,方法1：利用矩阵乘法和冒号运算,x=1 2 3 4;y=x*ones(1,3);y1=(y(:,),P176,方法2：采用求余函数mod,y = x(mod(n, M)+1)可实现对x(n) 以M为周期的周期延拓。加1是因为MATLAB向量下标只能从1开始。,11,%,第七章例,7.3,程序,q703,序列,x(n),的移位与周期延拓,%N:,观察窗口长度,%M:,序列,x(n),长度,%m:,移位样点数,clear;close all,N=24;M=8;m=input(,输入移位值：,m=);,if (m=N-M+1)%,检验输入参数,m,是否合理,error,(,输入数据不在规定范围内！,);,end,n=0:N-1;,x1=(0.8).n; x2=(n=0),xn=x1.*x2; %,产生,x(n),xm=zeros(1,N);%,设定,xm,的长度,for k=m+1:m+M,xm(k)=xn(k-m);,end,%xn移位m位，变为xm,xc=xn(mod(n,8)+1); %,产生 x(n)的周期延拓,xcm=xn(mod(n-m,8)+1); % 产生 x(n) 周期延拓后的移位,12,subplot(2,2,1); stem(n,xn,.);,title(x(n);axis(0,length(n),0,1);,ylabel(x(n);xlabel(n);,subplot(2,2,2);stem(n,xm,.);,title(x(n-m);axis(0,length(n),0,1);,ylabel(x(n-m);xlabel(n);,subplot(2,2,3); stem(n,xc,.);,title(x(n),的周期延拓,);axis(0,length(n),0,1);,ylabel(x(mod(n,8);xlabel(n);,subplot(2,2,4);stem(n,xcm,.);,title(x(n),的循环移位,);axis(0,length(n),0,1);,ylabel(x(mod(n-m,8);xlabel(n);,结果：,输入移位值：,m=3,13,14,如果，xcm=xn(mod(n-m,8)+1);改为： xcm=xm(mod(n,8)+1); %,产生 x(n)移位后的周期延拓,15,例7.4 离散系统对信号的响应,本题给定6阶低通数字滤波器的系统函数，求它在下列输入序列,x,(,n,)下的输出序列,y,(,n,)。,解：本题的计算原理见例6.14，在这里用工具箱函数filter来解。如果已知系统函数,H,(,z,)=,B,(,z,)/,A,(,z,)，则filter函数可求出系统对输入信号,x,(,n,)的响应,y,(,n,)。,y,= filter(B, A, x),P177例题7.4,由差分方程可得到,H,(,z,)的分子和分母多项式系数向量A和B，再给出输入向量x即可。,16,%,第七章例,7.4,程序,q704,求时域离散系统对常用序列的响应,%N:,输入数据长度,%B:,差分方程系数,%A:,差分方程系数,clear;close all;,N=64;n=0:N-1;m=10; %,下面设定系统参数,A,B,B=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),的响应,x3=n=0; %,产生输入信号,x3(n),y3=filter(B,A,x3); %,对,x3(n),的响应,x4=(n=0)%,产生输入信号,x4(n),y4=filter(B,A,x4); %,对,x4(n),的响应,x5=exp(j*pi*n/8),.*x4,;%,产生输入信号,x5(n),y5=filter(B,A,x5); %,对,x5(n),的响应,17,subplot(3,2,1);stem(n,y1,.);line(0,N,0,0),axis(0,N,min(y1),max(y1);ylabel(y1(n),subplot(3,2,2);stem(n,y2,.);line(0,N,0,0),axis(0,N,min(y2),max(y2);ylabel(y2(n),subplot(3,2,3);stem(n,y3,.);line(0,N,0,0),axis(0,N,min(y3),max(y3);ylabel(y3(n),subplot(3,2,4);stem(n,y4,.);line(0,N,0,0),axis(0,N,min(y4),max(y4);ylabel(y4(n),subplot(3,2,5);stem(n,real(y5),.);line(0,N,0,0),axis(0,N,-1,1);ylabel(Rey5(n),subplot(3,2,6);stem(n,imag(y5),.);line(0,N,0,0),axis(0,N,-1,1);ylabel(Imy5(n),18,运行结果：,B = 0.0004 0.0022 0.0056 0.0075 0.0056 0.0022 0.0004A = 1.0000 -3.1836 4.6223 -3.7795 1.8136 -0.4800 0.0544,19,例7.5 系统线性性质验证,设系统差分方程为,y,(,n,) =,x,(,n,) + 0.8,y,(,n,-1) ，,要求用程序验证系统的线性性质。,解：产生两种输入序列，分别乘以常数后,1. 分别激励系统，再求输出之和；,2. 先相加，再激励系统求输出；,对两个结果进行比较，方法是求它们之差，按,误差的绝对值是否极小,进行判断。,20,% Q705,时域离散系统线性性质验证,%N:,输入信号长度,%B:,差分方程系数,%A:,差分方程系数,clear;close all;,N=64;n=0:N-1;m=20;,B=1;A=1,-0.8;,%,设定系统参数,A,B,x1=0.8.n; %,产生输入信号,x1(n),x=(n=0),x1=x1.*x;,y1=filter(B,A,x1); %,对,x1(n),的响应,y1(n),x2=(n-m)=0;,y2=filter(B,A,x2); %,对,x2(n),的响应,y2(n),x3=5*x1+3*x2;,y3=filter(B,A,x3); %,对,5x1(n)+3x2(n),的响应,y3(n),y=5*y1+3*y2; % y(n)=5y1(n)+3y2(n),21,subplot(2,2,1);stem(n,y1,.);line(0,N,0,0),axis(0,N,min(y1),max(y1);ylabel(y1(n),subplot(2,2,2);stem(n,y2,.);line(0,N,0,0),axis(0,N,min(y2),max(y2);ylabel(y2(n),subplot(2,2,3);stem(n,y3,.);line(0,N,0,0),axis(0,N,min(y3),max(y3);xlabel(n);ylabel(y3(n),subplot(2,2,4);stem(n,y,.);line(0,N,0,0),axis(0,N,min(y),max(y);xlabel(n);ylabel(y(n),22,23,e=y-y3,figure;stem(e),24,例7.6 离散序列的卷积计算,给出两个序列 和 ,计算其卷积y(n)，并图示各输入输出序列。,解：在例6.4中，已经给出了直接调用MATLAB的卷积函数conv的方法，也给出了自编卷积计算程序的方法，要注意的是本例时间变量的设定和移位方法。在本例中，设定n为从零开始，向量x和h的长度分别为Nx=20和Nh=10；结果向量y的长度为length(y)=Nx+Nh,-,1。,25,% q706,序列卷积运算,p180,% Nx,：输入序列,x(n),长度,% Nh,：,h(n),长度,% m,：,x(n),的移位样点数,clear;close all,Nx=20;Nh=10;m=5;,n=0:Nx-1;,x1=(0.9).n; %,产生,x1(n),x2=zeros(1,Nx+m);,for k=m+1:m+Nx %,产生,x2(n)=x1(n-m),x2(k)=x1(k-m);,end,nh=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),26,Ny1=length(y1); %y1(n),长度,Ny2=length(y2); %y2(n),长度,subplot(3,2,1);,stem(n,x1,.);,xlabel(n);ylabel(x1(n);axis(0,Ny1,0,1);,subplot(3,2,3);,stem(nh,h1,.);,xlabel(n);ylabel(h1(n);axis(0,Ny1,0,1);,n1=0:Ny1-1;,subplot(3,2,5);,stem(n1,y1,.);,xlabel(n);ylabel(y1(n);axis(0,Ny1,0,max(y1);,nx=0:length(x2)-1;,subplot(3,2,2);,stem(nx,x2,.);,xlabel(n);ylabel(x2(n);axis(0,Ny2,0,1);,subplot(3,2,4);,stem(nh,h2,.);,xlabel(n);ylabel(h2(n);axis(0,Ny2,0,1);,n2=0:Ny2-1;,subplot(3,2,6);,stem(n2,y2,.);,xlabel(n);ylabel(y2(n);axis(0,Ny2,0,max(y2),27,28,7.2 Z变换和傅里叶变换,Z变换的定义：,有限长序列x(n)，n=ns:nf，其Z变换为：,29,7.2 Z变换和傅里叶变换,求z的逆变换的方法：,对于z变换分式,可以用,部分分式法,或,长除法,求其反变换。,用函数,residuez,可以求出它的极点留数分解,其中,r, p, k = residuez (B, A),其反变换为：,30,例7.7 有限长度序列的Z和逆Z变换,已知两序列,x,1,= 1,2,3，,n,1,= -1:1 及,x,2,= 2,4,3,5，,n,2,= -2:1，求出,x,1,与,x,2,及其卷积,x,的Z变换。,解： 序列,x,1,、,x,2,的Z变换可写成,两个多项式乘积 ，可用,conv,函数来求得序列x(n)。,conv,函数只能给出x数组，n数组要自己判别。n的起点ns = ns1 + ns2 =,3，终点nf = nf1 + nf2 = 2。 n=ns:nf。由x和n即可得出X(z)。,31,%,第七章例,7.7,程序,q707,%,计算,Z,变换,x1=1,2,3;ns1=-1; %,设定,x1,和,ns1,nf1=ns1+length(x1)-1;% nf1,可以算出,x2=2,4,3,5;ns2=-2;%,设定,x2,和,ns2,nf2=ns1+length(x1)-1;% nf2,可以算出,x=conv(x1,x2)%,求出,x,n=(ns1+ns2):(nf1+nf2)%,求出,n,程序运行结果：,x = 2 8 17 23 19 15,n = -3 -2 -1 0 1 2,由此可得,32,例7.8 求Z多项式分式的逆变换,设系统函数为，,输入例7.7中的,x,2,信号，用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变换，得,33,%第七章例7.8程序q708,% 计算逆Z变换,clear,close all,x=2,4,3,5;nsx=-2; % 输入序列及初始时间,nfx=nsx+length(x)-1; % 计算序列终止时间,Bw=-3;nsbw=-1; % 系统函数的分子系数，及z的最高次数,Aw=2,-2.2,.5;nsaw=0; %系统函数的分母系数，及z的最高次数,B=conv(-3,x); % 输入与分子z变换的多项式乘积,A=Aw; % 分母不变,nsy=,nsaw,-(,nsbw-nsx,) %,分母,比,分子,z的次数高nsy,r,p,k=residuez(B,A) %求留数r，极点p及直接项k,nf=input(终点时间nf=) %要求键入终点时间,n=min(nsx,nsy):nf ; % 生成总时间数组,34,% 求无限序列yi和直接序列yd,yi=(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; % 合成输出,xe=zeros(1,length(n); %初始化，将x延拓为xe,xe(find(n=nsx)&(n=nfx),=1,)=x; %在对应的n处把xe赋值x,subplot(2,1,1),stem(n,xe,.),line(min(n(1),0),nf,0,0) %绘图,subplot(2,1,2),stem(n,y,.),line(min(n(1),0),nf,0,0),35,由程序算出nsy =-1,留数：r = -57.7581 204.7581,极点：p = 0.7791 0.3209,直接项：k = -150 -30,代入,得,如果不需要求出表达式，可用filter函数求解。,36,x,(,n,),y,(,n,),37,例7.9 离散时间傅里叶变换,取一个周期的正弦信号，作8点采样，求它的连续频谱。然后对该信号进行,N,个周期延拓，再求它的连续频谱。把,N,无限增大，比较分析其结果。,解：离散时间傅里叶变换的定义：,MATLAB实现：设置一系列较密频率w,i,，求出一系列X(w,i,),X = x*exp(-j*n,*w),。,先把离散时间傅立叶变换写成一个子程序dtft.m，有了子程序，本例就没有什么难度了。,38,计算离散时间傅里叶变换的函数文件dtft.m：,function X=dtft(x,w),% 计算离散时间傅立叶变换,% X 为在w频率点上的DTFT数组,% x 为输入离散序列，时间数组为增序整数，可用下标表示,% w 为频率数组，由于exp(-j*n *w)以2为周期，,% 所以其范围通常选为-, 或0, 2,n=1:length(x);,X=x*exp(-j*n*w);,*注意：要恢复原来量纲，应把无量纲时间n乘以采样周期Ts，把无量纲频率w（数字频率）乘以采样频率Fs。,39,%第七章例7.9程序q709 % 离散时间傅立叶变换,clear all; close all;,disp(八点时间信号的离散时间傅立叶变换),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,subplot(3,1,1),plot(w,abs(X0),grid,shg% 画图,disp(重复N次的八点时间信号的离散时间傅立叶变换),N=input(N= );% 用键盘输入延拓周期数,x1=,reshape,(x0*ones(1,N),1,N*length(x0);% 延拓后的时域信号x1,X1 = dtft(x1,w)*dt;% 求x1的频率响应X1,subplot(3,1,2),plot(w,abs(X1),grid,shg% 画图,disp(重复无穷次的八点信号的离散时间傅立叶变换-傅立叶级数),pause,X2=fft(x0) *dt;% 离散傅立叶变换,w1=2*pi* 0:length(x0)-1/length(x0); % 离散频点向量,subplot(3,1,3),stem(-w1,w1,abs(X2),abs(X2),grid,shg,axis(min(w),max(w),0,max(abs(X2),grid,shg,40,程序运行结果,执行程序q709并按提示键入N = 4，所得图形如图7.10所示。N取得愈大，其峰值愈大，宽度愈窄。当N取得很大时，会出现内存不足的问题，这是用矩阵乘法做傅里叶变换的缺点。另外，因为那时峰值点处的宽度很窄，也会出现所选频点对不上峰值点的问题。所以对于N无限增大的情况，必须用fft函数来求。这时用连续频谱也没有意义了。这里用同样的横坐标把几种频谱进行对比，使读者更好地理解其关系。,41,图7.10,42,例7.10 时域采样频率与频谱混叠,分别以采样频率f,s,=1000Hz，400Hz和200Hz对x,a,(t)进行等间隔采样，计算并图示三种采样频率下的采样信号及其幅频特性。,解：程序分别设定4种采样频率f,s,=10kHz，1kHz，400Hz和200Hz，对x,a,(t)进行采样，得到采样序列xa(t)，xa1(n)，xa2(n)，xa3(n)，画出其幅度频谱。采样时间区间均为0.1秒。为了便于比较，画出了幅度归一化的幅频曲线，如图7.11所示。,43,例7.10 采样频率与频谱混叠(续),由于,由以上关系式可见，,采样信号的频谱函数是原模拟信号频谱函数的周期延拓，延拓周期为2,/T,。,如果以频率f为自变量（,= 2,f），则以采样频率f,s,= 1/T为延拓周期。,对频带限于f,c,的模拟信号x,a,(t)，只有当f,s,2f,c,时，采样后 才不会发生频谱混叠失真。这就是著名的,采样定理,44,%第七章例7.10程序q710% 时域采样及其频谱,clear;close all;,fs=10000; fs1=1000; fs2=400; fs3=200;% 设置四种采样频率,t=0:1/fs:0.1;% 采集信号长度为0.1秒,A=444.128; a=50*sqrt(2)*pi; b=a;,xa=A*exp(-a*t).*sin(b*t);,k=0:511; f=fs*k/512;%由wk=2k/512=2fT求得模拟频率f,Xa=dtft(xa,2*pi*k/512); % 近似模拟信号频谱,T1=1/fs1;t1=0:T1:0.1;% 采集信号长度为0.1秒,x1=A*exp(-a.*t1).*sin(b*t1);%1kHz采样序列x1(n),X1=dtft(x1,2*pi*k/512); % x1(n)的512点dtft,T2=1/fs2;t2=0:T2:0.1;% 采集信号长度为0.1秒,x2=A*exp(-a.*t2).*sin(b.*t2);%400Hz采样序列x2(n),X2=dtft(x2,2*pi*k/512); % x2(n)的512点dtft,T3=1/fs3;t3=0:T3:0.1;% 采集信号长度为0.1秒,x3=A*exp(-a.*t3).*sin(b.*t3); %200Hz采样序列x3(n) X3=dtft(x3,2*pi*k/512); % x3(n)的512点dtft,45,figure(1);,subplot(2,2,1);plot(t,xa);,axis(0,max(t),min(xa),max(xa);title(模拟信号);,xlabel(t(s);ylabel(Xa(t);line(0,max(t),0,0),subplot(2,2,2);plot(f,abs(Xa)/max(abs(Xa);,title(模拟信号的幅度频谱);axis(0,500,0,1),xlabel(f (Hz);ylabel(|Xa(jf)|);,subplot(2,2,3);stem(t1,x1,.);,line(0,max(t1),0,0);axis(0,max(t1),min(x1),max(x1),title(采样序列x1(n)（fs1=1kHz）);xlabel(n);ylabel(X1(n);,f1=fs1*k/512;,subplot(2,2,4);plot(f1,abs(X1)/max(abs(X1);,title(x1(n)的幅度谱);xlabel(f (Hz);ylabel(|X1(jf)|);,46,figure(2);,subplot(2,2,1);stem(t2,x2,.);,line(0,max(t2),0,0);axis(0,max(t2),min(x2),max(x2);,title(采样序列x2(n)（fs2=400Hz）);xlabel(n);ylabel(X2(n);,f=fs2*k/512;,subplot(2,2,2);plot(f,abs(X2)/max(abs(X2);,title(x2(n)的幅度谱);xlabel(f (Hz);ylabel(|X2(jf)|);,subplot(2,2,3);stem(t3,x3,.);,line(0,max(t3),0,0);axis(0,max(t3),min(x3),max(x3);,title(采样序列x3(n)（fs3=200Hz）);xlabel(n);ylabel(X3(n);,f=fs3*k/512;,subplot(2,2,4);plot(f,abs(X3)/max(abs(X3);,title(x3(n)的幅度谱);xlabel(f (Hz);ylabel(|X3(jf)|),47,48,例7.11 由离散序列恢复模拟信号,用时域内插公式,其中,模拟用理想低通滤波器恢复的过程，观察恢复波形，计算出最大恢复误差。,解：这个公式与卷积公式相像，可以用向量和矩阵乘法来解决。,49,例7.11 由离散序列恢复模拟信号,xa = x *,sinc(Fs * TNM),TNM=ones(length(n),1)*t- n*T*ones(1,length(t);,50,%第七章例7.11程序q711 % 时域采样与重构,clear; close all;,A=444.128; a=50*sqrt(2)*pi; b=a;,for k=1:2,if k=1 Fs=400;,elseif k=2 Fs=1000;end,T=1/Fs; dt=T/3; %每个采样间隔T上g(t)取三个样点,Tp=0.03; % 重构时间区间为0,0.03s,t=0:dt:Tp;% 生成序列t,n=0:Tp/T;% 生成序列n,TNM=ones(length(n),1)*t- n*T*ones(1,length(t);% 生成TNM矩阵,x=A*exp(-a.*n*T).*sin(b*n*T);% 生成模拟信号采样序列x(n),xa=x*sinc(Fs*TNM); %内插公式,subplot(2,1,k);plot(t,xa);hold on,axis(0,max(t),min(xa)-10,max(xa)+10);,st1=sprintf(由Fs= %d,Fs);st2=Hz的采样序列x(n)重构的信号;,st=st1,st2;title(st),ylabel(xa(t);,xo=A*exp(-a.*t).*sin(b*t);%以3Fs对原始模拟信号采样,stem(t,xo,.);line(0,max(t),0,0),emax2=max(abs(xa-xo),end,51,最大重构误差：emax2 = 27.7015,emax2 = 9.9436,52,例7.12 梳状滤波器零极点和幅频特性,梳状滤波器系统函数有如下两种类型。,FIR型：,IIR型：,freqz：数字滤波器频率特性计算和绘制函数,zplane：H(z)的零-极点图绘制。,zplane(B,A),解：调用函数freqz和zplane 很容易写出程序q712.m。,53,%第七章例7.12程序q712,% 梳状滤波器零极点和幅频特性,clear;close all,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);,54,figure(1);,subplot(2,2,1);,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梳状滤波器零、极点图,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(/),55,56,57,例7.13 低通滤波及时域卷积定理,输入信号,x,(,n,) = cos(0.04,n,) + cos(0.08,n,) + cos(0.4,n,) + 0.3,(,n,)，0,n,63 通过低通滤波器，计算滤波器对,x,(,n,)的响应输出,y,(,n,)，并图示,x,(,n,)和,y,(,n,)，观察滤波效果。,解：如前所述，只要求出,H,(,z,)=,B,(,z,)/,A,(,z,)的分子和分母多项式系数向量,B,和,A,，则可调用滤波器直接型实现函数,filter,对输入信号,x,(,n,)进行滤波。,y = filter(B, A, x),58,%,第七章例,7.13,程序,q713,% IIR,滤波器实现及,FT,的时域卷积定理验证,clear;close,all,%,产生输入信号,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;,%,求,H(z,),分子分母多项式系数向量,B,和,A,b=1,2,1; %(1+z-1)2,的展开系数,B=0.0003738*,conv(conv(b,b),b,);%,嵌套调用卷积函数,conv,a1=1,-1.2686,0.7051;,a2=1,-1.0106,0.3583;,a3=1,-0.9044,0.2155;,A=conv(conv(a1,a2),a3);,59,%,对,x(n,),滤波,y=,filter(B,A,x,);,%,绘图,subplot(3,2,1);stem(x,.),axis(0,max(n)/4,min(x),max(x);line(0,max(n),0,0),title(,输入信号,x(n);xlabel(n);ylabel(x(n,),subplot(3,2,3);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,),60,X=,fft(x,N,); %,计算,x(n,),的,N,点离散,DFT,Y=,fft(y,N,); %,计算,y(n,),的,N,点离散,DFT,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,4);plot(f,abs(Y),title(,输出信号,y(n,),的幅频曲线,);,xlabel(/);ylabel(|FTy(n,)|),axis(0,0.5,0,max(abs(Y);,61,H,f,=,freqz(B,A,N,whole,);,subplot(3,2,5);plot(f/pi,abs(H),title(,滤波器幅频响应曲线,);,xlabel(/);ylabel(H,幅度,),axis(0,0.5,0,max(abs(H);,Ym,=H.*X;,subplot(3,2,6);plot(f/pi,abs(Ym),title(|FTx(n)FTh(n)|);xlabel(/);ylabel(Ym,幅度,),axis(0,0.5,0,max(abs(Ym);,62,63,例7.14 用符号运算工具箱解z变换问题,解：无限长度时间序列的z变换和逆z变换都属于符号运算的范围。MATLAB的symbolic（符号运算）工具箱已提供了这种函数。如果读者已在计算机上安装了这个工具箱，可以键入以下程序。,MATLAB程序q714.m,其特点是程序的开始要指定符号自变量,syms z n a N w0% 规定z,n,a为符号变量,64,syms z n a N w0,y1=an,Y1=simplify(ztrans(y1),y2=n,Y2=simplify(ztrans(y2),y3=n*an,Y3=simplify(ztrans(y3),y4=n*(n-1)/2,Y4=simplify(ztrans(y4),y5=exp(j*w0*n),Y5=simplify(ztrans(y5),y6=sin(w0*n),Y6=simplify(ztrans(y6),pause,X1=z/(z-1),x1=iztrans(X1),X2=-3*z-1/(2-5*z-1+2*z-2),x2=iztrans(X2),X3=z/(z-a),x3=iztrans(X3),X4=z/(z-1)2),x4=iztrans(X4),X5=z/(z-1)3),x5=iztrans(X5),X6=(1-z-N)/(1-z-1),x6=iztrans(X6),X7=z/(z-exp(j*w0),x7=iztrans(X7),65,
展开阅读全文