第七讲MATLAB在电路中的应用

利用MATLAB进行信号与系统分析

信号的MATLAB表示信号的基本运算信号的能量和功率利用MATLAB对信号进行时域分析利用MATLAB对信号进行频域分析傅里叶(Fourier)变换基本信号的MATLAB表示

指数信号Aeat、指数序列ak、抽样函数Sa(t)、正弦型信号、矩形脉冲信号、三角脉冲信号一、信号的MATLAB表示一、信号的MATLAB表示1、连续时间信号的MATLAB表示连续时间信号：时间变化连续如y=x(t)离散时间信号(序列)：时间离散如x(nT)=x(t)|t=nT例：用MATLAB绘制连续时间信号关于t的曲线，t的取值范围为0-30s，并以0.1s递增。%file65.mt=0:0.1:30;x=exp(-0.707*t).*sin(2*t/3);plot(t,x);gridon;xlabel('Time(sec)');ylabel('x(t)');函数名功能函数名功能sawtooth(t,width)产生锯齿波或三角波信号(周期为2π，幅值－1~1)pulstran(t,d,’func’)产生冲激串square(t,duty)产生方波信号(周期为2π，幅值－1~1)rectpuls(t,width)产生非周期的方波信号sinc(t)产生抽样函数波形tripuls(t,width,skew)产生非周期的三角波信号chirp(t,f0,t1,f1,’method’,phi)产生调频余弦扫频信号diric(x,n)产生Dirichlet或周期sinc函数gauspuls(T,FC,BW,BWR)产生高斯正弦脉冲信号gmonopuls(t,fc)产生高斯单脉冲信号voc(x,fc,fs)电压控制振荡器2、工具箱中的信号产生函数

举例（1）%产生周期为0.02的三角波file66.mt=0:0.0001:1;x1=sawtooth(2*pi*50*t,0.8);plot(t,x1);axis([0,0.2,-1,1])%产生频率为50HZ，占空比为50%的方波file67.mt=0:0.0001:0.08;x1=square(2*pi*50*t,50);plot(t,x1);axis([0,0.08,-1,1.5])举例（2）%产生线性调频信号：file68.mt=0:0.001:2;y=chirp(t,0,1,150);figure(1);plot(t,y);axis([0,0.5,0,1])%产生非周期三角波信号：file69.mt=0:0.01:1;y1=tripuls(t);y2=tripuls(t,0.6);subplot(2,1,1);plot(t,y1);subplot(2,1,2);plot(t,y2);3、离散时间信号的表示

在MATLAB中，离散时间信号x(n)的表示：需要用一个向量x表示序列幅值，用另一个等长的向量来定位时间变量n来完整地表示一个序列。例：绘制离散时间信号的棒状图。其中x(-1)=-1,x(0)=1,x(1)=2,x(2)=1,x(3)=0,x(4)=-1。%离散时间信号的表示file70.mn=-3:5;%定位时间变量x=[0,0,-1,1,2,1,-1,0,0];stem(n,x);grid;%绘制火柴杆图line([-3,5],[0,0]);%画x轴线xlabel('n');ylabel('x[n]')4、几种常用离散时间信号的表示

1.单位脉冲序列x=zeros(1,N);x(1,n0)=1;2.单位阶跃序列n=[ns:nf];x=[(n-n0)>=0]3.实指数序列n=[ns:nf];x=a.^n4.复指数序列n=[ns:nf];x=exp((sigema+jw)*n)5.正(余)弦序列n=[ns:nf];x=cos(w*n+sita)%decayingexponentialt=0:001:10;%file71.mA=1;a=-0.4;ft=A*exp(a*t);plot(t,ft)t=0:0.5:10;A=1;a=-0.4;ft=A*exp(a*t);stem(t,ft)举例（1）一、信号的MATLAB表示%rectpulsfile72.mt=0:0.001:4;T=1;ft=rectpuls(t-2*T,T);plot(t,ft)axis([0,4,-0.5,1.5])一、信号的MATLAB表示举例（2）%tripuls

file73.mt=-3:0.001:3;ft=tripuls(t,4,0.5);plot(t,ft)ft1=tripuls(t,4,1);plot(t,ft1)一、信号的MATLAB表示举例（3）%unitimpulssequence

file75.mk=-50:50;delta=[zeros(1,50),1,zeros(1,50)];stem(k,delta)function[f,k]=impseq(k0,k1,k2)%产生

f[k]=delta(k-k0)；k1<=k<=k2k=[k1:k2];f=[(k-k0)==0];k0=0;k1=-50;k2=50;%file74.m[f,k]=impseq(k0,k1,k2);stem(k,f)一、信号的MATLAB表示举例（4）%unitstepsequencek=-50:50;%file76.muk=[zeros(1,50),ones(1,51)];stem(k,uk)function[f,k]=stepseq(k0,k1,k2)%产生

f[k]=u(k-k0);k1<=k<=k2k=[k1:k2];f=[(k-k0)>=0];k0=0;k1=-50;k2=50;[f,k]=stepseq(k0,k1,k2);stem(k,f)function[f,k]=stepseq(k0,k1,k2)%产生

f[k]=u(k-k0);k1<=k<=k2k=[k1:k2];f=[(k-k0)>=0];k0=0;k1=-50;k2=50;%file77.m[f,k]=stepseq(k0,k1,k2);stem(k,f)一、信号的MATLAB表示举例（5）求的实部序列，虚部序列，幅度序列，相位序列

fsxl(0.1,pi/3,-10,10)file78.mfunctionfsxl(delta,omig,n1,n2)e=delta+omig.*j;n=[n1:n2];x=exp(e*n);x_real=real(x);x_imag=imag(x);x_magnitude=abs(x);x_phase=(180/pi)*angle(x)subplot(2,2,1);stem(n,x_real);xlabel('n');ylabel('x_real');title('realpart');subplot(2,2,2);stem(n,x_imag);xlabel('n');ylabel('x_imag');title('imaginarypart');subplot(2,2,3);stem(n,x_magnitude);xlabel('n');ylabel('x_magnitude');title('amplitudepart');subplot(2,2,4);stem(n,x_phase);xlabel('n');ylabel('x_phase');title('phasepart');二、信号的基本运算%信号的翻转和时移file79.mt=-3:0.001:3;ft1=tripuls(2*t,4,0.5);subplot(2,1,1)plot(t,ft1)title('f(2t)')ft2=tripuls((2-2*t),4,0.5);subplot(2,1,2)plot(t,ft2)title('f(2-2t)')1.

信号的尺度变换、翻转、时移（平移）

已知三角波f(t)，用MATLAB画出的f(2t)和f(2-2t)波形

二、信号基本运算2.

信号的相加与相乘相加用算术运算符“+”实现相乘用数组运算符“.*”实现

例:画信号Aeatcos(w0t+f)的波形

t=0:0.001:8;%file80.mA=1;a=-0.4;w0=2*pi;phi=0;ft1=A*exp(a*t).*sin(w0*t+phi);plot(t,ft1)二、信号基本运算3.

序列移位与周期延拓运算序列移位：y(n)=x(n-m)MATLAB实现：y=x;ny=nx-m

file85.m序列周期延拓：y(n)=x((n))MMATLAB实现：ny=nxs:nxf；y=x(mod(ny,M)+1)file86.mn1=[-5:4];%序列x1(n)的起始及终止位置file81.mn1s=-5;n1f=4;x1=[231-13421-5-3];%序列x1(n)不同时间的幅度n2=[0:9];%序列x2(n)的起始及终止位置n2s=0;n2f=9;x2=[1111111111];%序列x2(n)不同时间的幅度ns=min(n1s,n2s);nf=max(n1f,n2f);%新信号的时间起始终止位置n=ns:nf;y1=zeros(1,length(n));%延拓序列初始化y2=zeros(1,length(n));y1(find((n>=n1s)&(n<=n1f)==1))=x1;%给延拓序列y1赋值x1y2(find((n>=n2s)&(n<=n2f)==1))=x2;%给延拓序列y2赋值x2ya=y1+y2;%逐点相加yp=y1.*y2;%逐点相乘subplot(4,1,1);stem(n,y1,'.');subplot(4,1,2);stem(n,y2,'.');subplot(4,1,3);stem(n,ya,'.');subplot(4,1,4);stem(n,yp,'.');%file82.m周期延拓N=24;M=8;m=3;n=0:N-1;x1=(0.8).^n;%产生指数序列x2=[(n>=0)&(n<=M)];%形成矩形序列Rm(n)x=x1.*x2;%截取操作形成新序列x(n)xm=zeros(1,N);fork=m+1:m+Mxm(k)=x(k-m);%产生移位序列x(n-3)Endxc=x(mod(n,M)+1);%产生x(n)的周期延拓x(n)8xcm=x(mod(n-m,M)+1);%产生x(n)的周期延拓x(n-3)8subplot(4,1,1);stem(n,x,'.');ylabel('x(n)');subplot(4,1,2);stem(n,xm,'.');ylabel('x(n-3)');subplot(4,1,3);stem(n,xc,'.');ylabel('x(n)_8');subplot(4,1,4);stem(n,xcm,'.');ylabel('x(n-3)_8')4.序列翻褶与序列累加运算序列翻褶数学描述：y(n)=x(-n)。MATLAB可实现:y=fliplr(x)序列累加的数学描述为：MATLAB实现：y=cumsum(x)5.两序列的卷积运算两序列卷积运算：MATLAB实现：y=conv(x1,x2)。序列x1(n)和x2(n)必须长度有限。二、信号基本运算6.两序列的相关运算两序列相关运算：MATLAB实现：y=xcorr(x1,x2)。例：求序列的翻褶序列y(n)=x(-n)和y(n)=x(-n+3)及累加和。file83.m二、信号基本运算n=0:10;x=3*exp(-0.2*n);y=fliplr(x);n1=-fliplr(n);n2=-fliplr(n-3);%在指定位置m=3处的时间序列翻褶s=cumsum(x);subplot(2,2,1);stem(n,x);xlabel('n');ylabel('x(n)');subplot(2,2,2);stem(n2,y);xlabel('n');ylabel('y(n)=x(-3+n)');subplot(2,2,3);stem(n1,y);xlabel('n');ylabel('y(n=)x(-n)');subplot(2,2,4);stem(n,s);xlabel('n');ylabel('s(n)');1.信号能量数字定义：MATLAB实现:E=sum(x.*conj(x));

或E=sum(abs(x).^2);数字定义：2.信号功率MATLAB实现:P=sum(x.*conj(x))/N;

或E=sum(abs(x).^2)/N;非周期三角波信号能量的MATLAB计算：file84.mdt=0.0001;t=0:dt:1;x=tripuls(t);E=sum(abs(x).^2*dt);余弦信号的平均功率MATLAB计算：file85.mdt=0.0001;t=0:dt:2*pi;x=cos(t);P=sum(abs(x).^2*dt)./(2*pi);一个信号处理就是将输入信号变换成输出信号的运算过程，如图所示。在此过程中，输出的信号称为系统对输入信号作出的响应，输入信号称为系统的激励信号。当输入信号为连续时间信号x(t)时，该系统称为连续时间处理系统，可描述为：y(t)=T[x(t)]当输入信号为离散时间信号x[n]时，该系统称为离散时间处理系统，可描述为：y[n]=T[x[n]]当一个系统具有可加性和齐次性，则该系统称为线性系统。如果系统响应与激励加于系统的时刻无关，则该系统为时不变系统。1.系统的描述(1)．常系数线性微分/差分方程(2)．系统传递函数连续系统：离散系统：系统传递函数的MATLAB表示：注意：在MATLAB中，传递函数由分子、分母两个多项式的系数来表示，系数为降幂排列。例如：可以表示为：num=[1,0.2,1];den=[1,0.5,0.7]如果有空项的话，需要用零来补齐。sys=tf(num,den)%获得连续系统的传递函数形式sys=tf(num,den,0,1,’variable’,’z^-1’)%获得离散系统的传递函数形式(3)．零－极点增益模型zplane连续系统：离散系统：(4)．极点留数模型离散系统：连续系统：在MATLAB中，增益系数、零点向量、极点向量用三个列向量表示：增益系数(Gain)：k零点向量(Zero)：z=[z1z2

…zn]’极点向量(Pole)：p=[p1p2

…pn]’sys=zpk(z,p,k)%获得零-极点模型表达式(6)．状态空间模型ss连续系统：离散系统：(5)．二次分式模型连续系统：离散系统：状态向量：x输出向量：y激励向量(输入向量)：u在MATLAB中，用矩阵A，B，C，D表示系统的状态空间模型。系统模型的MATLAB表示系统模型表示名称参数传递函数模型tfnum-传递函数的分子系数向量den-传递函数的分母系统向量零极点增益模型zpz-系统零点向量;p-系统极点向量;k-增益系数二次分式模型sossos-二次分数的系数矩阵;g-比例系数状态空间模型ss[A,B,C,D]-状态空间模型参数部分分式模型residuezR-分子向量;P-分母系数向量;K-常数在MATLAB中，用sos、ss、tf、zp分别表示二次分式模型、状态空间模型、传递函数模型和零-极点增益模型。其中sos表示二次分式，g为比例系数，sos为L×6的矩阵。2.系统模型的转换函数(1)．ss2tf函数

格式：[num,den]=ss2tf(A,B,C,D,iu)功能：将指定输入量iu的线性系统(A,B,C,D)转换为传递函数模型[num,den]。(2)．zp2tf函数格式：[num,den]=zp2tf(z,p,k)功能：将给定系统的零－极点增益模型转换为传递函数模型，z、p、k分别为零点列向量、极点列向量和增益系数。(3).tf2sos函数格式：[sos,g]=tf2sos(num,den)功能：将给定系统的传递函数模型[num,den]转换为二次分式模型sos,g为增益系数。线性系统模型的转换函数函数名功能说明函数名功能说明ss2tf状态空间模型转换为传递函数模型tf2ss传递函数模型转换为状态空间模型ss2zp状态空间模型转换为零－极点增益模型zp2ss零－极点增益模型转换为状态空间模型zp2tf零－极点增益模型转换为传递函数模型tf2zp传递函数模型转换为零－极点增益模型residue传递函数模型转换为极点留数模型residue极点留数模型转换为传递函数模型ss2sos状态空间模型转换为二次分式模型sos2ss二次分式模型转换为状态空间模型tf2sos传递函数模型转换为二次分式模型sos2tf二次分式模型转换为传递函数模型zp2sos零－极点增益模型转换为二次分式模型sos2zp二次分式模型转换为零－极点增益模型例：将系统转换为状态空间模型[A,B,C,D]在命令窗口输入：>>z=[3];p=[1,2];k=2;>>[A,B,C,D]=zp2ss(z,p,k);

屏幕显示为:A=3.0000-1.41421.41420B=10C=2.0000-4.2426D=0例：求离散时间系统的零、极点向量和增益系数。在命令窗口输入：>>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=23.系统互联与系统结构MATLAB实现函数series()

格式：[A,B,C,D]=series(A1,B1,C1,D1,A2,B2,C2,D2)或

[num,den]=series(num1,den1,num2,den2)将系统1、系统2级联，可得到级联连接的传递函数形式为：(1).系统的级联(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时表示负反馈。例：求两个单输入单输出子系统的级联、并联和反馈后系统的传递函数。file90.mMATLAB源程序为：file90.mnum1=1;den1=[1,1];%系统1num2=2;den2=[1,2];%系统2[nums,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因此，各系统的传递函数分别为：

而在实际应用中，也可以把一个复杂的线性时不变系统分解为几个简单系统的组合结构，即直接型结构、级联型结构和并联型结构。MATLAB所提供的系统模型变换函数，实质就是给出了这几种系统结构的互换关系。系统的传递函数对应于系统的直接型结构，二次分式(sos)模型对应于级联型结构，系统传递函数的部分分式(residue或residuez)形式对应于并联型结构。例：已知某数字滤波器的传递函数为：请求出其级联型结构和格型结构。file91.m%file87.mclearall;b=[2,13/12,5/4,2/3];a=[1];%设定参数fprintf('级联型结构系数：');[sos,g]=tf2sos(b,a);%直接型到级联型转换fprintf('格型结构系数(反射系数)：');[K]=tf2latc(b,a);%直接型到格型转换sos=1.00000.536001.0000001.00000.00570.62191.000000g=2格型结构系数（反射系数）：K=0.25000.50000.3333根据级联结构系数，其传递函数可以改写为：例：已经某数字滤波器的传递函数为

请求出其级联型结构和并联型结构%file88.m

b=[1,-3,11,-27,18];a=[16,12,2,-4,-1];disp('级联型结构系数');[sos,g]=tf2sos(b,a);disp('并联型结构系数');[R,P,K]=residuez(b,a);级联型结构系数：sos=1.0000-3.00002.00001.0000-0.2500-0.12501.00000.00009.00001.00001.00000.5000g=0.0625并联型结构系数：R=-5.0250-1.0750i-5.0250+1.0750i0.925027.1875P=

-0.5000+0.5000i-0.5000-0.5000i0.5000-0.2500

K=

-18根据级联结构系数，其传递函数可以改写为：根据并联结构系数，其传递函数可改写为：模型转换函数示意图(1)．连续LTI系统的响应(2)．离散LTI系统的响应用MATLAB中的卷积函数conv()来实现。用MATLAB中的卷积函数conv()来实现。%file89.mdt=input('输入时间间隔dt')x=ones(1,fix(10/dt));h=exp(-0.1*[0:fix(10/dt)]*dt);y=conv(x,h);t=dt*([1:length(y)]-1);plot(t,y);grid;例：某LT1系统的单位冲激响应，输入初始条件为零，求系统的响应y(t)。%file90.mx=ones(1,10);n=[0:14];h=0.5.^n;y=conv(x,h);stem(y);xlabel('n');ylabel('y[n]');例：已知某LT1离散系统的单位冲激响应为：求输入信号序列的系统响应。nx=-5:4;x=ones(1,10);nh=[0:14];h=0.5.^n;y=conv(x,h);n0=nx(1)+nh(1);%求卷积序列y起始时间位置N=length(nx)+length(nh)-2;%求卷积序列y序列长度ny=n0:n0+N;%求卷积序列y的时间向量subplot(2,2,1);stem(nx,x);title('x[n]');xlabel('n');xlabel('x[n]');subplot(2,2,2);stem(nh,h);title('h[n]');xlabel('n');xlabel('h[n]');subplot(2,2,3);stem(ny,y);title('x[n]*h[n]');xlabel('n');xlabel('y[n]');h=get(gca,'position');h(3)=2.5*h(3)set(gca,'position',h)改写为：file91.m(3)．时域响应函数(a)连续系统的单位冲激响应 调用格式：[Y,T]＝impulse(sys,t)或impulse(sys)功能：返回系统的响应Y和时间向量T，自动选择仿真的时间范围。其中sys可为系统传递函数、零极增益模型或状态空间模型。%file92.mb=[23];a=[10.41];sys=tf(b,a);[y,t,x]=impulse(sys);plot(t,y)%file93.ma=[-0.55,-0.78;0.78,0];b=[1;0];c=[1.96,6.45];d=[0];impulse(a,b,c,d);title('LTI系统的冲激响应')例：有一个二阶系统：求系统的单位冲激响应。(b)连续系统的单位阶跃响应 调用格式：格式：[Y,T]=step(sys,t)功能：返回系统的单位阶跃响应Y和仿真时间T，仿真时间范围为t，或不设则可自动选择。其中sys可为系统传递函数(tf)、零极增益模型(zpk)或状态空间模型(ss)。%file94.mb=[23];a=[10.41];step(b,a)(c)连续系统的零输入响应 调用格式：[y,t,x]=initial(sys,t,x0)计算出连续时间系统LT1由于初始状态x0所引起的零输入响应y，其中x为状态记录，t为仿真所用的采样时间向量。(d)对任意输入的连续LTI系统响应函数lsim() 格式：[y,x]=lsim(a,b,c,d,u,t)功能：返回连续LTI系统对任意输入时，系统的输出响应y和状态记录x，其中u给出每个输入的时序列，一般情况下u为一个矩阵；t用于指定仿真的时间轴，它应为等间隔。%file95.mnum=[2,5,1];den=[1,2,3];t=0:0.1:10peiod=4;u=(rem(t,peiod)>=peiod/2)lsim(num,den,u,t);title('方波响应');

例：有一个二阶系统：求当输入是周期为4s的方波时的输出响应。(e)离散系统(数字滤波器)的单位冲激响应 调用格式：[y,x]＝impz(num,den,n)功能：返回单位冲激响应y向量和时间状态历史记录x向量。%file96.mnum=[0.18080.10470.31060.1047];den=[1.0000-1.14801.5107-0.6991]tf(num,den,0.1,’variable’,’z^-1’)impz(num,den,50);%取50个采样点%也可用以下方法实现%file97.mn=50;imp=[1zeros(1,n-1)];y=filter(num,den,imp);stem(y);%file98.mclearall;num=[2,-3.5,1.5];den=[1,-1.7,0.3];dimpulse(num,den,10);title('离散LTI系统的冲激响应')例：有一个二阶系统：求系统的单位冲激响应。(f)离散系统(数字滤波器)的单位阶跃响应 调用格式：[y,x]=stepz(num,den)/dstep功能：返回多项式传递函数G(z)=num(z)/den(z)表示的系统单位阶跃响应。%file99.mb=[12];a=[121];stepz(b,a)(g)离散系统(数字滤波器)的零输入响应 调用格式：[y,x,n]=dinitial(a,b,c,d,x0)功能：计算离散时间LTI系统由初始状态x0所引起的零输入响应y和状态响应响应x，取样点数由函数自动选取。n为仿真所用的点数。(h)对任意输入的离散LTI系统响应函数dlsim()调用格式：[y,x]=dlsim(a,b,c,d,u)功能：返回离散LTI系统对输入序列u的响应y和状态记录x。%file100.mnum=[2,-3.4,5.5];den=[1,-1.2,0.8];u=randn(1,100);dlsim(num,den,u);title('随机噪声响应')

例：有一个二阶系统：求系统对100点随机噪声的响应曲线。(i)滤波函数filter()格式：y=filter(B,A,x)功能：对向量x中的数据进行滤波处理，即差分方程求解，产生输出序列向量y。B和A分别为数字滤波器系统函数H(z)的分子和分母多项式系数向量。连续系统的频率响应H(jω)的定义与物理意义：

Yf(jω)=H(jω)F(jω)LT1系统把频谱为F(jω)的输入改变成频谱为H(jω)F(jω)的响应，改变的规律完全由H(jω)决定。H(jω)称为系统的频率响应，定义为：幅度响应相位响应H(jω)与h(t)的关系：H(jω)反映了系统对输入信号不同频率分量的传输特性。由H(jω)的定义，显然有

H(jω)=F[h(t)]即H(jω)等于系统冲激响应h(t)的Fourier变换。1．求模拟滤波器Ha(s)的频率响应函数freqs()格式一：H＝freqs(B,A,W)

功能：计算由向量W(rad/s)指定的频率点上模拟滤器系统函数Ha(s)的频率响应Ha(jΩ)，结果存于H向量中。格式二：[H,w]＝freqs(B,A,M)功能：计算出M个频率点上的频率响应存于向量H中,M个频率存放在向量w中。freqs自动将这M个频率点设置在适当的频率范围。默认w和M时,freqs自动选取200个频率点计算。求该模拟滤波器的频率响应。MATLAB源程序如下：%file101.mB=1;A=[12.61313.41422.61311];W=0:0.1:2*pi*5;freqs(B,A,W)已知某模拟滤波器的系统函数

已知某滤波器的系统函数为2．求数字滤波器H(z)的频率响应函数freqz()格式一：H=freqz(B,A,W)功能：计算由向量W(rad)指定的数字频率点上(通常指在H(z)的频率响应H(ejw)。

格式二：[H,W]=freqz(B,A,N)格式三：[H,W]=freqz(B,A,N,‘whole’)求该滤波器的频率响应。MATLAB源程序为：%file102.mB=[10000000-1];A=1;freqz(B,A)该程序运行所绘出的幅频与相频性曲线如图所示。例：已知一系统传递函数为：求系统的单位冲激响应h(n),以及h(n)的幅频和相频响应。%file103.mN=128;x=[1,zeros(1,N-1)];%产生单位冲激函数num=[0.01-0.030.06-0.030.01];den=[12.531.80.5];y=dlsim(num,den,x);%计算单位冲激响应figure(1);n=1:N;stem(n,y);gridon;title('单位冲激响应');figure(2)freqz(num,den,N);%绘制幅频和相频响应曲线gridon;dimpulse(num,den)单位冲激响应函数，不带输出变量时，直接绘出单位冲激响应幅频和相频响应冲激响应3．滤波函数filter从频域角度，无论是连续时间LTI系统还是离散时间LTI系统，系统对输入信号的响应实质上就是对输入信号频谱进行不同选择处理的过程，这个过程称为滤波。格式：y=filter(B,A,x)功能：对向量x中的数据进行滤波处理，即差分方程求解，产生输出序列向量y。B和A分别为数字滤波器系统函数H(z)的分子和分母多项式系数向量。设系统的差分方程为求该系统对信号的响应。MATLAB源程序为：file104.mB=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)该程序运行所得结果如图所示时域信号特性频域特性变换名称非周期连续信号非周期连续频谱傅里叶变换周期性连续信号非周期离散频谱傅里叶级数非周期离散信号周期连续频谱序列傅里叶变换周期离散信号周期离散频谱离散傅里叶级数离散信号（有限样本点）周期离散频谱离散傅里叶变换傅里叶变换：是建立以时间为自变量的“信号”与以频率为自变量的“频谱函数”之间的某种对应关系。遵循以下原则：时域连续周期离散非周期非周期离散周期连续频域1.连续时间、连续频率－傅里叶变换正变换：

逆变换：

这说明求和的问题可以用x(t)行向量乘以列向量来实现。式中是t的增量，在程序中用dt表示.例：求如图矩形脉冲信号f(t)在-40~40rad/s区间的频谱。%file105.mclc;clearall;N=input('取时间分隔的点数N=');dt=10/N;t=[1:N]*dt;%给出时间分割f=[ones(1,N/2),zeros(1,N/2)];%给出信号(方波)wf=input('需求的频谱宽度wf=');Nf=input('需求的频谱点数Nf=');w1=linspace(0,wf,Nf);F1=f*exp(-j*t'*w1)*dt;%求傅里叶变换w=[-fliplr(w1),w1(2:Nf)];%补上负频率F=[fliplr(F1),F1(2:Nf)];%补上负频率对应的频谱subplot(121)plot(t,f,'linewidth',1.5)gridon;axis([0,10,0,1.1]);subplot(122)plot(w,abs(F),'linewidth',1.5);gridon;取时间分隔的点数N=64需求的频谱宽度wf=40需求的频谱点数Nf=256取时间分隔的点数N=256需求的频谱宽度wf=40需求的频谱点数Nf=64采样较密采样稀，有频率泄漏2连续时间、离散频率－傅里叶级数正变换：

逆变换：

式中为离散频率相邻两谱线间的角频率间隔，k为谐波序号。3离散时间、离散频率－离散傅里叶级数正变换：

逆变换：

5离散时间、离散频率－离散傅里叶变换(DFT)正变换：

逆变换：

4时间离散、连续频率－序列傅里叶变换正变换：

逆变换：

注意，对于序列傅里叶变换，如果x(n)为无限长，那么就不能用MATLAB直接利用上式计算，只可以用它对表达式在频率点上求值，再画出它的幅度和相位。如果x(n)为有限长，就可直接用MATLAB计算5离散时间、离散频率－离散傅里叶变换(DFT)在MATLAB中，旋转因子矩阵可以由dftmtx函数实现。dftmtx函数:调用格式：W=dftmtx(n)功能：函数调用后，返回一个n*n的离散傅里叶变换DFT的旋转因子矩阵W。这样对于给定长度为n的行向量信号x(n),利用Xk=x*W就可以获得x(n)的傅里叶变换X(k)。其实，旋转因子矩阵也可以自己通过表达式构造。例如：编写函数，首先利用MATLAB计算序列x(n)的DFT和IDFT函数.functionxn=idft(xk)N=length(xk);n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);WNnk=WN.^(-n'*k);xn=xk*WNnk/N;IDFT:functionxk=dft(xn)N=length(xn);WNnk=dftmtx(N);xk=xn*WNnk;DFT:WNnk=conj(dftmtx(N))DFT参数对频率分辨率的影响：频率分辨率就是，在频率轴上能分辨出的两个频率点的最小间隔要能区分两个频率点f1,f2,，有效数据长度N必须满足：可以对序列进行补零操作，以增加数据的长度，这样做虽然不能提高频率分辨率，但补零后对原来的X(k)起到了插值作用，可以平滑频谱的包络。但是，如果增加有效数据的长度，则可以提高频率分辨率。利用上面编写的DFT函数，求的频谱特性，其中取fs=20Hz%file106.mN=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);k=(0:N/2-1)*fs/Nsubplot(211);plot(n,xn);gridon;title('x(n)(0=<n<128)')subplot(212);plot(k,abs(xk(1:N/2)));gridon;%file107.mN=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);N1=256;n1=0:N1-1;xn1=[xn,zeros(1,N1-N)];xk1=dft(xn1);m=(0:N1/2-1)*fs/N1subplot(211);plot(n1,xn1);gridon;title('x(n)补零后');subplot(212);plot(m,abs(xk1(1:N1/2)));gridon;title(‘x(n)补零后’)%file108.mN=402;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);k=(0:N/2-1)*fs/Nsubplot(211);plot(n,xn);gridon;title('x(n)(0=<n<401)')subplot(212);plot(k,abs(xk(1:N/2)));gridon;

%file109.mN=128;n=0:N-1;fs=20;xn=sin(3.1*2*pi*n/fs)+cos(3*2*pi*n/fs);xk=dft(xn);subplot(211)plot(abs(xk));gridon;subplot(212)k=(0:N-1)*fs/N;plot(k,abs(xk));gridon;(1)．一维快速正傅里叶变换函数fft格式：X=fft(x,N)功能：采用FFT算法计算序列向量x的N点DFT变换，当N缺省时，fft函数自动按x的长度计算DFT。当N为2整数次幂时，fft按基-2算法计算，否则用混合算法。(2)．一维快速逆傅里叶变换函数ifft格式：x=ifft(X,N)功能：采用FFT算法计算序列向量X的N点IDFT变换。(3).将零频分量移至频谱中心的函数格式：Y=fftshift(X)功能：用来重新排列X=fft(x)的输出，把X的左右两半进行交换，从而将零频分量移至频谱中心。%file110.mfs=1000;

%采样频率t=0:1/fs:0.5;%时间轴取值范围及间隔x=cos(160*pi*t);N=512;

%FFT点数subplot(311)；plot(t,x);gridon;y=fft(x,N);

%N点FFT[x(t)]xk=fftshift(abs(y));

%求幅度谱并将零频分量移至频谱中心位置f=(-N/2:N/2-1)*fs/512;

%频率范围subplot(312)plot(f,xk);%做频谱图gridon;%用IFFT恢复原始信号z=real(ifft(y,N));ti=[0:length(z)-1]/fs;subplot(313)plot(ti,z);gridon;例：求的FFT序列和IFFT序列用快速傅里叶变换FFT计算下面两个序列的卷积。MATLAB程序(部分)：file111.m%线性卷积xn=sin(0.4*[1:15]);

%对序列x(n)赋值,M=15hn=0.9.^(1:20);

%对序列h(n)赋值,N=20ticyn=conv(xn,hn); %直接调用函数conv计算卷积tocsubplot(2,2,1);stem(xn);subplot(2,2,2);stem(hn);subplot(2,2,3);stem(yn);h=get(gca,'position');h(3)=2.3*h(3)set(gca,'position',h)%园周卷积L=pow2(nextpow2(M+N-1));

ticXk=fft(xn,L);

Hk=fft(hn,L);

Yk=Xk.*Hk;

yn=ifft(Yk,L);

toc

x(n),h(n)及其线性卷积波形Elapsedtimeis0.005074seconds.Elapsedtimeis0.000072seconds.例题部分例1求系统

y"(t)+2y'(t)+100y(t)=10f(t)的零状态响应，已知f(t)=sin(2pt)u(t)。%file112.mts=0;te=5;dt=0.01;sys=tf([1],[12100]);t=ts:dt:te;f=10*sin(2*pi*t);y=lsim(sys,f,t);plot(t,y);xlabel('Time(sec)')ylabel('y(t)')例2求系统

y"(t)+2y'(t)+100y(t)=10f(t)的零状态响应，已知f(t)=d(t)。%连续时间系统的冲激响应file113.m

ts=0;te=5;dt=0.01;sys=tf([10],[12100]);t=ts:dt:te;y=impulse(sys,t);plot(t,y);xlabel('Time(sec)')ylabel('h(t)')例3

分析噪声干扰的信号f[k]=s[k]+d[k]通过M点滑动平均系统的响应，

其中s[k]=(2k)0.9k是原始信号，d[k]是噪声。File114.mR=51;d=rand(1,R)-0.5;k=0:R-1;s=2*k.*(0.9.^k);f=s+d;figure(1);plot(k,d,'r-.',k,s,'b--',k,f,'g-');M=5;b=ones(M,1)/M;a=1;y=filter(b,a,f);figure(2);plot(k,s,'b--',k,y,'r-');噪声干扰信号f[k]=s[k]+d[k]通过M点滑动平均系统的响应例4求系统y[k]+3y[k-1]+2y[k-2]=10f[k]的单位脉冲响应。%

离散系统的单位脉冲响应file115.m

k=0:10;a=[132];b=[1];h=impz(b,a,k);stem(k,h)根据卷积公式：因此编程的过程为：（1）写出h(t)＝e-0.1t的MATLAB表达式；（2）写出u(t)的MATLAB表达式；（3）利用MATLAB的卷积语句y=conv(u，h)求解（4）画曲线plot(t,y)。例5计算x[k]*y[k]并画出卷积结果，已知x[k]={1,2,3,4;k=0,1,2,3}，

y[k]={1,1,1,1,1;k=0,1,2,3,4}。%file116.mx=[1,2,3,4];y=[1,1,1,1,1];z=conv(x,y);N=length(z);stem(0:N-1,z);例6方波分解为多次正弦波之和图示的周期性方波，其傅里叶级数为分别计算直到9次谐波，并做图。file118.m输入周期性方波--x利用MATLAB进行信号的频域分析

周期信号频谱的MATLAB实现

用数值积分分析非周期信号频谱一、周期信号频谱的MATLAB实现

频谱Fn一般为复数，可分别利用abs和angle函数获得其幅频特性和相频特性。其调用格式分别为x=abs(Fn)y=angle(Fn)

周期信号的频谱Fn为离散信号，可以用stem例1试用MATLAB画出图示周期三角波信号的频谱。解：周期信号的频谱为画三角波信号频谱的MATLAB程序%file119.mN=8;n1=-N:-1;%计算n=-N到-1的Fourier系数

c1=-4*j*sin(n1*pi/2)/pi^2./n1.^2;c0=0;%计算n=0时的Fourier系数n2=1:N;%计算n=1到N的Fourier系数

c2=-4*j*sin(n2*pi/2)/pi^2./n2.^2;cn=[c1c0c2];n=-N:N;subplot(2,1,1);stem(n,abs(cn));ylabel('Cn的幅度');subplot(2,1,2);stem(n,angle(cn));ylabel('Cn的相位');xlabel('\omega/\omega0');程序运行结果例2求周期矩形脉冲的Fourier级数表示式。并用MATLAB求出由前N项Fourier级数系数得出的信号近似波形。取A=1,T=2,t=1,w0=p

%Gibbsphenomenon%file120.mt=-2:0.001:2;N=input('Numberofharmonics=');c0=0.5;cN=c0*ones(1,length(t));

%dccomponent

forn=0:2:N-1%evenharmonicsarezero

cN=cN+cos(pi*n*t).*sinc(n*pi/2);

endplot(t,cN);%GibbsphenomenonN=5N=15N=50N=500二、用数值积分分析非周期信号频谱数值函数积分quad8可用来计算非周期信号频谱F

是一个字符串，它表示被积函数的文件名。a,b

分别表示定积分的下限和上限quad8的返回是用自适应Simpson算法得出的积分值y=quad8('F',a,b)例3试用数值方法近似计算三角波信号的频谱F(jw)=Sa2(w

/2)解:图示三角波可表示为三角波信号频谱的理论值为例3试用数值方法近似计算三角波信号的频谱functiony=sf1(t,w);

y=(t>=-1&t<=1).*(1-abs(t)).*exp(-j*w*t);%file121.mw=linspace(-6*pi,6*pi,512);N=length(w);F=zeros(1,N);fork=1:NF(k)=quad('sf1',-1,1,[],[],w(k));end

figure(1);plot(w,real(F));title('')xlabel('\omega');ylabel('F(j\omega)');figure(2);plot(w,real(F)-sinc(w/2/pi).^2);xlabel('\omega');title('计算误差');例3试用数值方法近似计算三角波信号的频谱运行结果利用MATLAB进行系统频域分析连续系统频响特性的计算周期信号通过系统的响应离散系统频响特性的计算一、连续系统频响特性的计算b

分子多项式系数

a

分母多项式系数

w

需计算的H(jw)的抽样点

(数组w中少需包含两个w的抽样点)。计算频响的MATLAB函数H=freqs(b,a,w)一、连续系统频响特性的计算例1三阶归一化的Butterworth低通滤波器的系统函数为%file122.mw=linspace(0,5,200);b=[1];a=[1221];h=freqs(b,a,w);subplot(2,1,1);plot(w,abs(h));subplot(2,1,2);plot(w,angle(h));

试画出|H(jw)|和(w)。一、连续系统频响特性的计算三阶Butterworth低通滤波器的幅度响应和相位响应二、周期信号通过系统的响应例2周期方波通过RC系统的响应。

二、周期信号通过系统的响应例2周期方波通过RC系统的响应。%p5_2PeriodicsignalpassLTIsystemfile123.mT=4;w0=2*pi/T;RC=0.1;t=-6:0.01:6;N=51;c0=0.5;xN=c0*ones(1,length(t));

%dcforn=1:2:N%evenharmonicsarezero

H=abs(1/(1+j*RC*w0*n));phi=angle(1/(1+j*RC*w0*n));

xN=xN+H*cos(w0*n*t+phi)*sinc(n*0.5);

endplot(t,xN);xlabel(['timeRC=',num2str(RC)]);grid;set(gca,'xtick',[-5-3-10135]);二、周期信号通过系统的响应例2周期方波通过RC系统的响应。三、离散系统频率响应的计算计算频率响应的MATLAB函数b

分子的系数

a

分母系数w

抽样的频率点(至少2点)，

w在0~2p之间幅频特性:abs,相频特性:angle

h=freqz(b,a,w)

b=[1];w=linspace(0,2*pi,512);h2=freqz(b,a2,w);plot(w/pi,abs(h1),w/pi,abs(h2),':');a1=[1-0.9];a2=[10.9];h1=freqz(b,a1,w);legend('\alpha=0.9','\alpha=-0.9');%file124.m三、离散系统频响特性的计算利用MATLAB进行连续系统的复频域分析部分分式展开的MATLAB实现H(s)的零极点与系统特性的MATLAB计算一、部分分式展开的MATLAB实现[r,p,k]=residue(num,den)num,den分别为F(s)分子多项式和分母多项式的系数向量。r为部分分式的系数，p为极点，k为多项式二、H(s)的零极点与系统特性的MATLAB计算计算多项式根roots的函数可用于计算H(s)的零极点。r=roots(N)%计算多项式N的根H(s)零极点分布图可用pzmap函数画出，调用形式为pzmap(sys)表示画出sys所描述系统的零极点图。例1试画出系统

的零极点分布图，求其单位冲激响应h(t)和频率响应H(j)，并判断系统是否稳定。%file125.mnum=[1];den=[1221];sys=tf(num,den);poles=roots(den)figure(1);pzmap(sys);t=0:0.02:10;h=impulse(num,den,t);figure(2);plot(t,h)title('ImpulseRespone')[H,w]=freqs(num,den);figure(3);plot(w,abs(H))xlabel('\omega')title('Magnitu')

利用MATLAB进行离散系统的Z域分析部分分式展开的MATLAB实现H(z)的零极点与系统特性的MATLAB计算一、部分分式展开的MATLAB实现[r,p,k]=residuez(num,den)num,den分别为F(z)分子多项式和分母多项式的系数向量。r为部分分式的系数，p为极点，k为多项式的系数。若为真分式，则k为零。二、H(z)的零极点与系统特性的MATLAB计算利用t