matlab滤波器.ppt

第九章MATLAB在信号和系统中的应用,9.1基本信号的表示利用Matlab软件信号处理工具箱(SignalProcessingToolbox)中的专用函数产生信号并绘出波形。,9.1.1连续时间信号的表示RECTPULSSampledaperiodicrectanglegenerator.RECTPULS(T)generatessamplesofacontinuous,aperiodic,unity-heightrectangleatthepointsspecifiedinarrayT,centeredaboutT=0.Bydefault,therectanglehaswidth1.Notethattheintervalofnon-zeroamplitudeisdefinedtobeopenontheright,i.e.,RECTPULS(-0.5)=1whileRECTPULS(0.5)=0.,rectpuls函数：产生矩形波调用格式为：x=rectpuls(t,width)例1：生成幅度为2，宽度T=4，中心在t=0的矩形波x(t)以及x(t-T/2)Fs=10000;%采样频率t=-4:1/Fs:4;T=4;x1=2*rectpuls(t,T);subplot(121);plot(t,x1);title(x(t);axis(-4602.2);x2=2*rectpuls(t-T/2,T);subplot(122);plot(t,x2);title(x(t-T/2);axis(-4602.2);,方波发生函数：square()SQUARESquarewavegeneration.SQUARE(T)generatesasquarewavewithperiod2*PifortheelementsoftimevectorT.SQUARE(T)islikeSIN(T),onlyitcreatesasquarewavewithpeaksof+1to-1insteadofasinewave.SQUARE(T,DUTY)generatesasquarewavewithspecifieddutycycle.Thedutycycle,DUTY,isthepercentoftheperiodinwhichthesignalispositive.,例2：产生50Hz占空比分别为20和50的矩形波。Fs=10000;%采样频率t=0:1/Fs:0.1;%采样间隔1/Fsf=50;%50Hzx1=square(2*pi*50*t,20);x2=square(2*pi*50*t,50);subplot(211);plot(t,x1);subplot(212);plot(t,x2);,三角波或锯齿波发生函数：sawtooth()语法格式：sawtooth(t,width)。产生以2为周期幅值范围在-1，+1之间的三角波或锯齿波。参数t为时间向量；width是0，1之间的数，它决定函数在一个周期内上升部分和下降部分的比例。width=0.5产生三角波，width=1产生锯齿波，此时函数可简写为：sawtooth(t)。,例3：产生f=50Hz的锯齿波、三角波Fs=10000;%采样频率t=0:1/Fs:0.1;%采样间隔1/Fsf=50;%50Hzx1=sawtooth(2*pi*50*t,0);x2=sawtooth(2*pi*50*t,1);x3=sawtooth(2*pi*50*t,0.5);subplot(311);plot(t,x1);subplot(312);plot(t,x2);subplot(313);plot(t,x3);,sinc函数：sinc()语法格式：sinc(t)t=linspace(-10,+10,200);x=sinc(t);plot(t,x);,1、产生正弦波t=0:0.01:3*pi;y=sin(2*t);plot(t,y)2、产生矩形脉冲信号t=-3:0.01:3;y=rectpuls(t-1,2);%以t=1为中心，宽度为2plot(t,y),3、产生周期方波t=0:0.01:1;y=square(4*pi*t);plot(t,y)4、产生周期锯齿波t=0:0.001:2.5;y=sawtooth(2*pi*30*t);plot(t,y)axis(00.2-11),5、产生指数函数波形x=linspace(0,1);y=exp(-x);plot(x,y)6、产生抽样信号t=-10:0.1:10;%向量t时间范围t=t1:p:t2,p为时间间隔f=sinc(t);%sinc(t)=sin(pi*t)/pi*tplot(t,f);%显示该信号的时域波形title(f(t)=Sa(t);xlabel(t)axis(-10,10,-0.4,1.1),7、产生单位阶跃信号定义阶跃函数functionf=heaviside(t,t0)%产生一个阶跃序列%t0代表起始时刻f=(t-t0)0);调用阶跃函数t=-1:0.01:3;f=heaviside(t,0);plot(t,f);axis(-1,3,-0.2,1.2);,8、产生单位冲激信号定义冲激函数functionchongji(t1,t2,t0)dt=0.01;t=t1:dt:t2;n=length(t);x=zeros(1,n);x(1,(-t0-t1)/dt+1)=1/dt;%表示在t=t+t0时刻存在冲激stairs(t,x);axis(t1,t2,0,1.2/dt);title(单位冲击信号(t)调用chongji(-1,5,0);可以试着给别的t1,t2,t0.,9、均匀分布的随机信号x=rand(1,N);%产生0,1上均匀分布的随机信号N=200;x=rand(1,N);plot(x);,10、高斯分布的随机信号x=randn(1,N);%产生均值为0，方差为1的高斯分布随机信号（即白噪声信号）。N=200;x=randn(1,N);plot(x);,例4：信号产生举例t=0:0.0001:0.1;x1=sawtooth(2*pi*50*t);%在0，0.1之间产生5个周期的锯齿波subplot(2,2,1),plot(t,x1),gridon;x2=sawtooth(2*pi*50*t,0.5);%在0，0.1之间产生5个周期的三角波subplot(2,2,2),plot(t,x2),gridon;,x3=square(2*pi*50*t);%在0，0.1之间产生5个周期的方波subplot(2,2,3),plot(t,x3),gridon;axis(0,0.1,-1.2,1.2);t=-4:0.1:4;x4=sinc(t);%产生抽样函数subplot(2,2,4),plot(t,x4),gridon;,9.1.2离散时间信号的表示在MATLAB中，离散时间信号x(n)的表示：需用一个向量x表示序列幅值，用另一个等长的定位时间变量n，才能完整地表示一个序列。例5：绘制离散时间信号的棒状图。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=0,x(4)=-1。MATLAB源程序为：n=-3:5;%定位时间变量x=0,0,-1,1,2,1,-1,0,0;stem(n,x);gridon;%绘制棒状图line(-3,5,0,0);%画x轴线xlabel(n);ylabel(xn)运行结果如图所示。,几种常用离散时间信号的表示1单位脉冲序列直接实现：x=zeros(1,N);x(1,n0)=1;2单位阶跃序列直接实现：n=ns:nf;x=(n-n0)=0;,3实指数序列直接实现：n=ns:nf;x=a.n;4复指数序列直接实现：n=ns:nf;x=exp(sigema+jw)*n);5正(余)弦序列直接实现：n=ns:nf;直接实现：n=ns:nf;x=cos(w*n+sita);,生成上述三种信号t=-5:0.01:5;a=2;y1=a.t;subplot(2,2,1);plot(t,y1,r);a=2;theat=pi/3;y2=sin(2*pi*t+theat)subplot(2,2,2);plot(t,y2);w=4;y3=exp(a+j*w)*t);subplot(2,1,2);plot(t,y3,y);,9.2信号的时域运算与时域变换1、加(减)、乘运算要求二个信号序列长度相同.t=0:0.01:2;f1=exp(-3*t);f2=0.2*sin(4*pi*t);f3=f1+f2;f4=f1.*f2;subplot(2,2,1);plot(t,f1);title(f1(t);subplot(2,2,2);plot(t,f2);title(f2(t);subplot(2,2,3);plot(t,f3);title(f1+f2);subplot(2,2,4);plot(t,f4);title(f1*f2);,2、用matlab的符号函数实现信号的反折、移位、尺度变换、周期延拓.由f(t)到f(-at+b)(a0)步骤:,例6：已知f(t)=sin(t)/t,试通过反折、移位、尺度变换由f(t)得到f(-2t+3)的波形。symst;f=sin(t)/t;%定义符号函数f(t)=sin(t)/tf1=subs(f,t,t+3);%对f进行移位f2=subs(f1,t,2*t);%对f1进行尺度变换f3=subs(f2,t,-t);%对f2进行反褶subplot(2,2,1);ezplot(f,-8,8);gridon;%ezplot是符号函数绘图命令subplot(2,2,2);ezplot(f1,-8,8);gridon;subplot(2,2,3);ezplot(f2,-8,8);gridon;subplot(2,2,4);ezplot(f3,-8,8);gridon;注:也可用一条指令:subs(f,t,-2*t+3)实现f(t)到f(-2t+3)的变换,例7：已知信号用matlab求f(t+2),f(t-2),f(-t),f(2t),-f(t),并绘出时域波形。symstf=sym(t/2+1)*(heaviside(t+2)-heaviside(t-2);subplot(2,3,1);ezplot(f,-3,3);y1=subs(f,t,t+2);subplot(2,3,2),ezplot(y1,-5,1);y2=subs(f,t,t-2);subplot(2,3,3),ezplot(y2,-1,5);y3=subs(f,t,-t);subplot(2,3,4),ezplot(y3,-3,3);y4=subs(f,t,2*t);subplot(2,3,5),ezplot(y4,-2,2);y5=-f;subplot(2,3,6),ezplot(y5,-3,3),t=-5:0.01:5;y1=heaviside(t,-2);y2=heaviside(t,2);f=(1+t./2).*(y1-y2)plot(t,f)functionf=heaviside(t,t0)%产生一个阶跃序列%t0代表起始时刻f=(t-t0)0);,离散信号运算信号的相加与相乘y(n)=x1(n)+x2(n)y(n)=x1(n)x2(n)MATLAB实现：y=x1+x2;y=x1.*x2序列移位：y(n)=x(n-m)。MATLAB实现：y=x;ny=nx-m序列周期延拓：y(n)=x(n)MMATLAB实现：ny=nxs:nxf；y=x(mod(ny,M)+1)序列翻褶：y(n)=x(-n)。MATLAB可实现:y=fliplr(x)序列累加的数学描述为：MATLAB实现：y=cumsum(x)两序列的相关运算MATLAB实现：y=xcorr(x1,x2),例8：N=24;M=8;m=3;n=0:N-1;x1=(0.8).n;%生成指数序列x2=(n=0),3.线性卷积运算已知离散信号x(n)和h(n)，求y(n)=x(n)*h(n)，并用图形表示,Y=conv(x,h)实现x,h二个序列的卷积。若x(n)和y(n)的长度分别为M和N，则返回值是长度为M+N-1的矢量。例9:求二个信号的卷积t11=0;t12=1;t21=0;t22=2;t1=t11:0.001:t12;ft1=2.*rectpuls(t1-0.5,1);subplot(2,2,1),plot(t1,ft1),gridon;axis(0304)t2=t21:0.001:t22;ft2=t2;subplot(2,2,2),plot(t2,ft2),gridon;axis(0304)t3=t11+t21:0.001:t12+t22;ft3=conv(ft1,ft2)ft3=ft3*0.001subplot(2,2,3),plot(t3,ft3),gridon;axis(0304),例10：已知某线性移不变系统单位抽样响应为系统输入信号为求系统输出。n=0:5x1=(n-2)=0);%x1=zeros(1,6);%x1(1,3)=1;x2=(n=0)-(n=3);%x2=(n-2)Nerror(Nmustnotbelessthanlengthofx1);endiflength(x2)Nerror(Nmustnotbelessthanlengthofx2);end,%以上语句判断两个序列的长度是否小于Nx1=x1,zeros(1,N-length(x1);%填充序列x1(n)使其长度为N1+N2-1(序列x1(n)的长度为N1，序列x2(n)的长度为N2)x2=x2,zeros(1,N-length(x2);%填充序列x2(n)使其长度为N1+N2-1n=0:1:N-1;x2=x2(mod(-n,N)+1);%生成序列x2(-n)NH=zeros(N,N);forn=1:1:NH(n,:)=cirshiftd(x2,n-1,N);%该矩阵的k行为x2(k-1-n)Nendyc=x1*H;%计算循环卷积,functiony=cirshiftd(x,m,N)%directlyrealizecircularshiftforsequencex%y=cirshiftd(x,m,N);%x:inputsequencewhoselengthislessthanN%m:howmuchtoshift%N:circularlength%y:outputshiftedsequenceiflength(x)Nerror(thelengthofxmustbelessthanN);endx=x,zeros(1,N-length(x);n=0:1:N-1;y=x(mod(n-m,N)+1);,ft4=circonv(ft1,ft2,length(t2);ft4=ft4*0.001;subplot(2,2,4),plot(t2,ft4),gridon;ft44=circonv(ft1,ft2,length(t3);ft44=ft44*0.001;figure,plot(t3,ft44),gridon;,9.3信号的能量和功率1、信号能量数字定义：MATLAB实现:E=sum(x.*conj(x);E=sum(abs(x).2);2、信号功率数字定义：MATLAB实现:P=sum(x.*conj(x)/N;P=sum(abs(x).2)/N;,9.4傅里叶(Fourier)变换连续时间、连续频率傅里叶变换正变换：逆变换：连续时间、离散频率傅里叶级数正变换：逆变换：,时间离散、连续频率序列傅里叶变换正变换：逆变换：离散时间、离散频率离散傅里叶级数正变换：逆变换：,离散时间、离散频率离散傅里叶变换(DFT)正变换：逆变换：,各种形式的Fourier变换,DTFT,DFS,CFS,CTFT,对称关系时域周期性频域离散性时域离散性频域周期性时域非周期频域连续性时域连续性频域非周期傅立叶变换的对偶性,计算信号的傅里叶变换（1）符号解法调用函数：fourier、ifourier例11：求函数的傅立叶变换。symstF=fourier(exp(-2*abs(t)F=4/(4+w2)例12：求傅立叶逆变换symstwf=ifourier(1/(1+w2),t)f=1/2*exp(-t)*Heaviside(t)+1/2*exp(t)*Heaviside(-t),例13：设画出f(t)及其幅频谱。symst;f=1/2*exp(-2*t)*sym(heaviside(t);F=fourier(f);subplot(211);ezplot(f);subplot(212);ezplot(abs(F),符号解法的局限：运算结果是符号表达式，必须使用ezplot来作图。如果结果中含有冲激函数等项，那么ezplot无法作图。返回结果有可能包含一些不能直接表达的式子，甚至会出现“未被定义函数或变量”项，更加无法作图。很多场合，尽管连续，但却不可能表示成符号表达式。,（2）数值解法从定义出发：用MATLAB作数值计算，它不能计算无限区间。可以根据波形情况，确定积分上下限，将t分N等份，用求和代替积分。,同样，在频域内也是取一系列的样本点，然后计算出不同w处的F值，然后就得到了傅立叶变换的数值解。而不同w处的F值都用上面同一个公式计算，因此可以使用MATLAB的矩阵计算功能。将w设为一个行向量，则:其中，F是与w等长的行向量，t是列向量，w是行向量，t*w是一个矩阵，其行数与t相同，列数与w相同。,例14：求门信号的傅立叶变换。dt=0.02;tao=2;t=-tao:dt:tao;f=Heaviside(t+1)-Heaviside(t-1);w1=10*pi;N=500;k=0:N;w=k*w1/N;F=f*exp(-j*t*w)*dt;F=real(F);w=-fliplr(w),w(2:501);%双边谱F=fliplr(F),F(2:501);subplot(211),plot(t,f)subplot(212),plot(w,F),1、离散傅里叶变换（DFT）在离散信号与系统分析中，DFT有着频谱分析仪的特性。对连续信号与系统通过时域采样，也可以使用DFT进行谱分析。因此，DFT在频谱分析中起了核心的作用。另一方面，从技术上而言，DFT运算存在着高效的快速算法FFT，因此信号的频谱分析中将主要使用MATLAB的fft函数以及ifft函数。,使用DFT对连续信号进行谱分析工程实际应用中，常会遇到连续时间信号，其频谱也是连续函数。为了利用DFT对信号进行频谱分析，须先对进行时域采样，得到再对进行DFT，得到的则是的傅里叶变换在频域区间上的N点等间隔采样。对带限信号来说，X(k)与在范围内上的采样值除相差一个因子外完全相同。以上说明构成了用DFT对连续信号进行谱分析的理论基础。,用DFT进行连续信号频谱分析时注意点(1)若信号持续时间有限长，则其频谱无限宽。(2)若信号频谱有限，则其信号持续时间必定无限长。(3)实际应用中，所能观察到的信号长度总是有限的。(4)因此，用DFT对连续时间信号进行频谱分析时，采样后产生的频谱混叠失真是不可避的。(5)得到的谱分析结果总是近似的，近似的程度取决于信号带宽、采样率以及信号的截取长度。,2、一维快速离散Fourier变换：fft()语法格式：y=fft(x)。功能：采用FFT算法计算序列向量x的N点DFT变、换。当x为矩阵时，计算x中每一列信号的离散傅里叶变换。当N缺省时，fft函数自动按x的长度计算DFT。当x的长度为2的幂时，用基-2算法；否则，采用较慢的分裂基算法。y=fft(x,n)。功能：计算n点的FFT。当x的长度大于n时，截断x;否则补零。,3、一维快速离散Fourier逆变换：ifft()语法格式：y=ifft(x)。采用FFT算法计算序列向量X的N点IDFT变换，y是计算信号x的快速离散傅里叶变换的逆变换。y=ifft(x,n)。计算n点的快速离散傅里叶变换的逆变换。4、离散余弦变换（DCT）：dct()语法格式：y=dct(x)。计算信号x的离散余弦变换。y=dct(x,n)。计算n点的离散余弦变换。当x的长度大于n时，截断x；否则补零。离散余弦逆变换可由函数idct实现。,例15：计算信号,n=1,2,100的DCT。用MATLAB语言可实现如下：N=100;n=1:N;x=n+20*sin(2*pi*n/20);y=dct(x);z=idct(y);subplot(311);stem(x,.);ylabel(原始信号);gridon;subplot(312);stem(y,.);ylabel(DCT信号);gridon;subplot(313);stem(z,.);ylabel(IDCT信号);gridon;,利用离散傅里叶变换（DFT）分析信号的频谱例16：已知序列x(n)=2sin(0.48n)+cos(0.52n)0n100,试绘制x(n)及它的离散傅里叶变换|X(k)|图。MATLAB实现程序：N=100;n=0:N-1;xn=2*sin(0.48*pi*n)+cos(0.52*pi*n);XK=fft(xn,N);magXK=abs(XK);phaXK=angle(XK);,subplot(1,2,1);plot(n,xn);xlabel(n);ylabel(x(n);title(x(n)N=100);gridon;subplot(1,2,2);k=0:length(magXK)-1;stem(k,magXK,.);xlabel(k);ylabel(|X(k)|);title(X(k)N=100);gridon;,5、利用FFT实现线性卷积若序列x1(n)、x2(n)为长度分别为N1、N2的有限长序列，yc(n)=x1(n)x2(n)，yl(n)=x1(n)*x2(n)。由DFT的性质可知：当NN1+N2-1时有yl(n)=yc(n)=IDFTDFTx1(n)DFTx2(n)。序列较长时DFT运算通常用快速算法FFT实现。在MATLAB的信号处理工具箱中函数FFT和IFFT用于快速傅里叶变换和逆变换。,例17：已知两序列求它们的线性卷积yl(n)=h(n)*x(n)和N点的圆周卷积y(n)=h(n)Nx(n),并研究两者之间的关系。,n=0:1:11;m=0:1:5;N1=length(n);N2=length(m);xn=0.8.n;%生成x(n)hn=ones(1,N2);%生成h(n)yln=conv(xn,hn);%直接用函数conv计算线性卷积ycn=circonv(xn,hn,N1);%用函数circonv计算N1点圆周卷积ny1=0:1:length(yln)-1;ny2=0:1:length(ycn)-1;,subplot(2,1,1);%画图stem(ny1,yln);ylabel(线性卷积);subplot(2,1,2);stem(ny2,ycn);ylabel(圆周卷积);axis(0,16,0,4);,用FFT实现两序列的线性卷积。实现程序：n=0:1:11;m=0:1:5;N1=length(n);N2=length(m);xn=0.8.n;%生成x(n)hn=ones(1,N2);%生成h(n),N=N1+N2-1;XK=fft(xn,N);HK=fft(hn,N);YK=XK.*HK;yn=ifft(YK,N);ifall(imag(xn)=0)stem(x,yn,.),9.5线性时不变系统9.5.1系统的描述1常系数线性微分/差分方程2系统传递函数tf连续系统：离散系统：,3零极点增益模型zpk连续系统：离散系统：4极点留数模型连续系统：离散系统：,5二次分式模型sos连续系统：离散系统：6状态空间模型ss连续系统：离散系统：,1、传递函数表示法在Matlab中，传递函数用分子、分母两个多项式的系数表示，系数为降幂排列。分子（Numerator）：B=b(1)b(2)b(m+1)分母（Denominator）：A=a(1)a(2)b(n+1)sys=zpk(num,den)%获得传递函数模型表达式,例：num=10.21;den=10.51;sys1=tf(num,den)%连续系统的传递函数形式sys2=tf(num,den,0.1,variable,z-1)%离散系统的传递函数形式,Transferfunction:s2+0.2s+1-s2+0.5s+1Transferfunction:1+0.2z-1+z-2-1+0.5z-1+z-2Samplingtime:0.1,2、零极点模型表示法zplane在Matlab中，增益系数、零点向量、极点向量用三个列向量表示。增益系数（Gain）：k零点向量（Zero）：z=z1z2zn极点向量（Pole）：p=p1p2pnsys=zpk(z,p,k)%获得零-极点模型表达式,例：k=3;z=234;p=567;sys=zpk(z,p,k)Zero/pole/gain:3(s-2)(s-3)(s-4)-(s-5)(s-6)(s-7),3、状态空间模型表示法：ss连续系统状态空间方程：离散系统状态空间方程：状态向量：x输出向量：y激励向量（输入向量）：u在Matlab中，用矩阵A、B、C、D表示系统的状态空间模型。,9.5.2系统模型的转换函数在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递函数模型和零极点增益模型。其中sos表示二次分式，g为比例系数，sos为L6的矩阵，即,1ss2tf函数格式：num,den=ss2tf(A,B,C,D,iu)功能：将指定输入量iu的线性系统(A,B,C,D)转换为传递函数模型num,den。2zp2tf函数格式：num,den=zp2tf(z,p,k)功能：将给定系统的零极点增益模型转换为传递函数模型，z、p、k分别为零点列向量、极点列向量和增益系数。,线性系统模型的变换函数,例18：求离散时间系统的零、极点向量和增益系数。num=2,3;den=1,0.4,1;num,den=eqtflength(num,den);%使长度相等z,p,k=tf2zp(num,den)屏幕显示为z=0-1.5000p=-0.2000+0.9798i-0.2000-0.9798ik=2,例19：将以下系统转换成零-极点模型、部分分式模型b=23;a=10.41;z,p,k=tf2zp(b,a)z,p,k=residuez(b,a),z=-1.5000p=-0.2000+0.9798i-0.2000-0.9798ik=2z=1.0000-1.3268i1.0000+1.3268ip=-0.2000+0.9798i-0.2000-0.9798ik=,例20：将以下系统转换为状态空间模型b=023;121;a=10.41;A,B,C,D=tf2ss(b,a),A=-0.4000-1.00001.00000B=10C=2.00003.00001.60000D=01,9.5.3系统互联与系统结构1、系统的级联MATLAB实现函数series()格式：A,B,C,D=series(A1,B1,C1,D1,A2,B2,C2,D2)或num,den=series(num1,den1,num2,den2)将系统1、系统2级联，可得到级联连接的传递函数形式为：,2、系统的并联MATLAB实现函数parallel()格式：A,B,C,D=parallel(A1,B1,C1,D1,A2,B2,C2,D2)num,den=parallel(num1,den1,num2,den2)将系统1、系统2并联，可得到并联连接的传递函数形式为：,3、两个系统的反馈连接函数feedback格式：A,B,C,D=feedback(A1,B1,C1,D1,A2,B2,C2,D2,sign)num,den=feedback(num1,den1,num2,den2,sign)将系统1和系统2进行反馈连接，sign表示反馈方式(默认值为-1)；当sig=+1时表示正反馈；当sig=-1时表示负反馈。,例19：求两个单输入单输出子系统的级联、并联和反馈后系统的传递函数。MATLAB源程序为：num1=1;den1=1,1;%系统1num2=2;den2=1,2;%系统2nums,dens=series(num1,den1,num2,den2)%实现两个系统级联nump,denp=parallel(num1,den1,num2,den2)%实现两个系统并联numf,denf=feedback(num1,den1,num2,den2)%实现两个系统反馈,程序运行结果为：nums=002；dens=132nump=034；denp=132numf=012；denf=134因此，各系统的传递函数分别为：,9.5.4线性时不变系统的响应1、线性时不变系统的时域响应连续LTI系统的响应用MATLAB中的卷积函数conv()来实现。2、离散LTI系统的响应用MATLAB中的卷积函数conv()来实现。,3、时域响应函数(1)对任意输入的连续LTI系统响应函数lsim()格式：y,x=lsim(a,b,c,d,u,t)功能：返回连续LTI系统对任意输入时系统的输出响应y和状态记录x，其中u任意输入信号，一般情况下u为一个矩阵；t为等间隔的时间向量,指明要计算响应的时间点.,(2)对任意输入的离散LTI系统响应函数dlsim()格式：y,x=dlsim(a,b,c,d,u)功能：返回离散LTI系统对输入序列u的响应y和状态记录x。,4、求连续LTI系统的单位冲激响应函数impulse()格式：Y,T=impulse(sys)功能：返回系统的单位阶跃响应Y和仿真所用的时间向量T。自动选择仿真的时间范围。其中sys可为系统传递函数、零极增益模型或状态空间模型。impulse(sys,t)t:为等间隔的时间向量,指明要计算响应的时间点。绘出系统在0t时间范围内冲激响应的时域波形,num=23;den=10.41;impulse(num,den),求离散系统的单位冲激响应函数dimpulse()impz函数求离散系统（数字滤波器）的单位冲击响应。（注：Matlab7.0不再支持dimpulse函数）格式：y,x=dimpulse(num,den)功能：返回项式传递函数的单位冲激响应y向量和时间状态历史记录x向量。,num=0.18080.10470.31060.10470.1808;den=1-1.14801.5107-0.69910.2703;y,t=impz(num,den,50);impz(num,den,50);%50个采样点也可以用以下方法实现：n=50;imp=1zeros(1,n-1);y=filter(b,a,imp);stem(y);,5、求连续LTI系统的零输入响应函数initial()格式：y,t,x=initial(sys,x0)功能：计算出连续时间LTI系统由于初始状态x0所引起的零输入响应y。其中x为状态记录，t为仿真所用的采样时间向量。求离散系统的零输入响应函数dinitial()格式：y,x,n=dinitial(sys,x0)功能：计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x，取样点数由函数自动选取。n为仿真所用的点数。,6、求连续系统的单位阶跃响应函数step()格式：Y,T=step(sys)功能：返回系统的单位阶跃响应Y和仿真所用的时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数(TF)、零极增益模型(ZPK)或状态空间模型(SS)。绘出系统在0t时间范围内冲激响应的时域波形。step(num,den,T)绘出系统在0T时间范围内冲激响应的时域波形T:为等间隔的时间向量,指明要计算响应的时间点.,num=23;den=10.41;step(num,den),求离散系统的单位阶跃响应函数dstep()stepz()格式：y,x=dstep(num,den)功能：返回多项式传递函数G(z)=num(z)/den(z)表示的系统单位阶跃响应。,num=0.18080.10470.31060.10470.1808;den=1-1.14801.5107-0.69910.2703;stepz(num,den,50),例20：已知描述某连续系统的微分方程为：试用Matlab绘出该系统冲激响应和阶跃响应；输入信号为,该系统零状态响应y（t）。b=12;a=121;subplot(1,3,1);impulse(b,a),gridon;%冲激响应subplot(1,3,2);step(b,a),gridon;%阶跃响应p=0.5;%定义取样时间间隔t=0:p:5;%定义时间范围x=exp(-2*t);%定义输入信号subplot(1,3,3);lsim(b,a,x,t),gridon;%对系统的输出信号进行仿真,例21:对系统分别求脉冲响应、阶跃响应及对输入u(t)=sin(t)的响应.num=1,1;den=1,1.3,0.8;T=0:0.1:3;y1=impulse(num,den,T);y2=step(num,den,T);U=sin(T);y3=lsim(num,den,U,T);subplot(1,3,1);plot(T,y1);title(脉冲响应)subplot(1,3,2);plot(T,y2);title(阶跃响应)subplot(1,3,3);plot(T,y3);title(输入为u=sint的响应),9.6线性时不变系统的频率响应频响特性：系统在正弦激励下稳态响应随信号频率变化的特性.|H(j)|:幅频响应特性.():相频响应特性(或相移特性).,9.6.1连续系统的复频域分析1.系统函数H(s)定义系统零状态响应的拉氏变换与激励的拉氏变换之比.H(s)=R(s)/E(s)在matlab中,传递函数描述法是通过传递函数分子和分母关于s降幂排列的多项式系数来表示的.例如,某系统传递函数：则可用如下二个向量num(numerator)和den(denominator)来表示:num=1,1den=1,1.3,0.8,2、系统零、极点分布与系统稳定性关系系统函数H(s)集中表现了系统的性能,研究H(s)在S平面中极点分布的位置,可很方面地判断系统稳定性.(1)稳定系统：H(s)全部极点落于S左半平面(不包括虚轴),即若满足，则系统是稳定的.(2)不稳定系统：H(s)极点落于S右半平面,或在虚轴上具有二阶以上极点,则在足够长时间后,h(t)仍继续增长,系统是不稳定(3)临界稳定系统:H(s)极点落于S平面虚轴上,且只有一阶,则在足够长时间后,h(t)趋于一个非零数值或形成一个等幅振荡.H(s)零、极点可用matlab的多项式求根函数roots求得：极点:p=roots(den)零点:z=roots(num)根据p和z，用plot命令即可画出系统零极点分布图,进而分析判断系统稳定性。,例22:系统函数画出系统零、极点分布图,判断该系统稳定性。num=1,0,-4;den=1,2,-3,2,1;p=roots(den)z=roots(num)plot(real(p),imag(p),*);holdon;plot(real(z),imag(z),o);gridon由系统零极点分布图可知,该系统有一对极点位于s右半平面,故系统是不稳定的.,1Matlab求系统频响特性函数freqs的调用格式:h=freqs(num,den,)即：求模拟滤波器Ha(s)的频率响应函数功能：计算由向量(rad/s)指定的频率点上模拟滤器系统函数Ha(s)的频率响应Ha(j)，结果存于H向量中。例23：已知某模拟滤波器的系统函数求该模拟滤波器的频率响应。MATLAB源程序如下。B=1;A=12.61313.41422.61311;W=0:0.1:2*pi*5;freqs(B,A,W),例24:求系统的频响特性.num=1,1;den=1,1.3,0.8;W=0:0.1:100;h=freqs(num,den,W);subplot(2,1,1);plot(W,abs(h);title(幅频特性)axis(0,20,0,1.5);set(gca,xtick,0,10,20);set(gca,ytick,0,1/sqrt(2),1.25);gridon;subplot(2,1,2);plot(W,angle(h);title(相频特性)axis(0,20,-pi/2,0.2);set(gca,xtick,0,10,20);set(gca,ytick,-pi/2,-pi/4,0);gridon;,2求数字滤波器H(z)的频率响应函数freqz()格式：H=freqz(B,A,W)功能：计算由向量W(rad)指定的数字频率点上(通常指在H(z)的频率响应H(ejw)。例25：已知某滤波器的系统函数为求该滤波器的频率响应。MATLAB源程序为：B=100000001;A=1;freqz(B,A)该程序运行所绘出的幅频与相频性曲线如图所示。,3滤波函数filter格式：y=filter(B,A,x)功能：对向量x中的数据进行滤波处理，即差分方程求解，产生输出序列向量y。B和A分别为数字滤波器系统函数H(z)的分子和分母多项式系数向量。例：设系统差分方程求该系统对信号的响应。MATLAB源程序为：B=1;A=1,-0.8;n=0:31;x=0.8.n;y=filter(B,A,x);subplot(2,1,1);stem(x)subplot(2,1,2);stem(y)该程序运行所得结果如图所示。,9.7信号采样与重构,例26：信号的采样与重构。t1=-15:0.5:15;f1=sinc(t1/pi);wm=1;%信号带宽wc=wm;%滤波器截至频率ts=pi/wm%采样间隔，临界采样ws=2*pi/ts;%采样角频率n=-100:100;%时域采样点数目nts=n*ts;%时域采样点f=sinc(nts/pi);%采样信号dt=0.005;t=-15:dt:15;fa=ts*f*wc/pi*sinc(wc/pi)*(ones(length(nts),1)*t-nts*ones(1,length(t);%信号重构error=abs(fa-sinc(t/pi);%重构信号与原信号的误差subplot(311);stem(t1,f1);title(Sa(t)=sinc(t/pi);subplot(312);plot(t,fa)title(Sa(t)信号的重构);grid;subplot(313);plot(t,error);,1、利用MATLAB画出信号的波形。取t04；,2、设p1=-2，p2=-30；p1=-2，p2=3针对极点参数,画出系统零极点分布图,判断系统稳定性.针对极点参数,绘出系统脉冲响应曲线,并观察t时,脉冲响应变化趋势.针对极点参数,绘出系统的频响曲线.,3、设信号用MATLAB绘出f(t+4),f(-2t+3),-f(t)的时域波形。4、已知连续时间信号利用MATLAB求的时域波形图。,5、设系统传递函数为试求在输入激励下系统的稳态响应。6、计算下列函数的拉普拉斯正、反变换：（1）（2）,试画出系统的零极点图分布图、求系统的单位脉冲与阶跃响应、频率响应判断系统是否稳定,9.8数字滤波器设计滤波：从复杂的信号中提取所需要的信号，抑制不需要的部分。保留信号中的有用频率成分，去除信号中的无用频率成分。数字滤波器的频率响应函数幅度响应：相位响应：,数字滤波器设计包括了无限冲激响应(IIR)和有限冲激响应（FIR）滤波器设计。MATLAB为此提供了多种成熟算法的相应函数以及丰富的设计工具。,数字滤波器的系统函数H(z)用如下式表示:在MATLAB中，用向量b，a来表示滤波器的系数b(i)和a(i)。,数字滤波器分类:递归系统IIR非递归系统FIR,数字滤波器的单位冲激响应定义为输入为单位样本序列时数字滤波器的响应，即：h(n)=T(n)其中：,为了形象地看到滤波器的性能，需要绘出滤波器的幅频和相频特性曲线。时域分析（对任意输入的响应）频域分析（幅值响应、相位响应）。滤波器时域和频域分析是设计各类滤波器、评价滤波器性能以及滤波器应用的基础，,9.8.1滤波器特性及使用函数1、频率响应MATLAB数字信号处理工具箱有很多函数提供对模拟和数字滤波器的频率响应分析。其中，freqz函数和freqs函数分别返回数字和模拟滤波器的频率响应。工具箱中通常使用的单位频率是Nyquist频率，即采样频率的1/2。注意：就数字滤波器函数来说，其频域指标中的所有频率都以Nyquist频率进行归一化。因此Nyquist频率也称归一化频率。,2、关于Nyquist频率的说明例如：系统采样频率为1000Hz，则若数字滤波器的截止频率等于300Hz，经Nyquist频率归一化后，其归一化频率就是300/500=0.6。若将归一化频率转换成数字信号处理教科书中所使用的数字频率(rad)，需乘以；反之，若乘以采样频率fs的一半，则将归一化频率转换回了模拟域频率(Hz)。归一化频率f应满足0b=1;%分子系数向量b(i)a=1-0.9;%分母系数向量a(i)如果用filter函数实现对信号x滤波，只要调用：y=filter(b,a,x);就可给出输入x经过滤波以后的输出y。,filtfilt函数实现零相位前向与后向结合的滤波。调用格式为：y=filtfilt(b,a,x)式中，b,a分别为滤波器传递函数H(z)的分子和分母多项式系数。x为滤波器的输入，为数值向量。y为滤波器的输出。该函数对序列x进行正常的正向滤波后，将滤波后的输出翻转重新用该滤波器进行滤波，第二次滤波后的输出序列的翻转即得到零相位