第10章 MATLAB应用实例

1、2022年年4月月20日星期三日星期三1第第1010章章MATLABMATLAB应用实例应用实例2022年年4月月20日星期三日星期三2第10章 MATLAB应用实例 学习目标 学会利用MATLAB进行单位换算； 了解MATLAB在绘制曲线图、曲面方程图，模拟抛物运动，运算器放大电路等方面的应用。2022年年4月月20日星期三日星期三310.1单位换算 编写一个程序，能把用户输入的数值进行单位换算，在这里最起码能实现对长度，面积，体积，重量的单位换算。 这是一个非常基本的问题，如果我们使用其他编程工具，可能很容易实现。 在这里我们并不只是为了要完成这个程序，而是为了通过这个程序来了解如何进行m

2、atlab程序设计。2022年年4月月20日星期三日星期三4首先要进行建立模型 在这里分两步实现 第一步首先给出菜单选择，在这里让用户选择要实现什么单位换算 第二步把输入的转换为一个统一的单位，例如长度单位都统一转换为米，面积都转换为平方米；第三步再转换为我们所要的单位。 为什么要这样做，这是因为各个单位之间的转换的变换常数个各不相同，如果直接进行单位换算的话将会十分麻烦。另外，把变换常数直接表示为一个矩阵，选择单位的序号，也就是单位的序号的也就成了矩阵的下标了。这样程序就比较的简单。2022年年4月月20日星期三日星期三5MATLAB程序代码（ex101.m） clear all;Check

3、=0;while(Check=3) disp(单位换算程序) fprintf(1)长度 2)面积n) Check=input(选择进行何种单位换算：); %如果Check等于3的话就退出程序 switch Check case 1 fprintf(长度单位：n); fprintf(1)米 2)分米 3)厘米 4)毫米n); fprintf(5)英寸 6)英尺 7)英里 8)码n); InNum=input(选择输入单位); OutNum=input(选择输出单位); %设定各种单位对米的变换常数数组ToMeter ToMeter=1.00 0.1 0.01 0.001 0.0254 0.304

4、8 1609.3 0.914; 2022年年4月月20日星期三日星期三6FrmMeter=1./ToMeter;%反变换常数数组FrmMeter为ToMeter数组的倒数 Value=input(输入待变换的值：n); ValueinM=Value*ToMeter(InNum); %把输入值变换为米 NewValue=ValueinM*FrmMeter(OutNum);%把米变换为输出单位 fprintf(变换后的值是%gn,NewValue); %打印变换后的值 case 2 fprintf(面积单位：n); fprintf(1)平方米 2)公顷 3)亩 4)平方英尺n); fprintf(

5、5)平方码 6)英亩 7)平方英寸n); InNum=input(选择输入单位); 2022年年4月月20日星期三日星期三7OutNum=input(选择输出单位); %设定各种单位对平方米的变换常数数组ToArea ToArea=1.00 10000.0 2000/3 0.0929 0.836 4047 0.000645; FrmArea=1./ToArea;%反变换常数数组FrmArea为ToArea数组的倒数 Value=input(输入待变换的值：n); ValueinA=Value*ToArea(InNum); %把输入值变换为平方米 NewValue=ValueinA*FrmAre

6、a(OutNum); %把平方米变换为输出单位 fprintf(变换后的值是%gn,NewValue); %打印变换后的值 otherwise endend2022年年4月月20日星期三日星期三810.2曲线图的绘制曲线图的绘制 三次抛物线的方程为 试探讨参数a，b和c对其图形的影响。cxbxaxy232022年年4月月20日星期三日星期三9建模 对于这个函数是我们经常见到的一元三次函数，我们可以直接写入绘图语句中。 注意坐标的设定方法，以得到适于观察的图形。 同时了解如何使用绘图函数来实现图形，以及如何在图形上作标注字符。 给出的程序不是唯一的，例如可用plot函数编程，读者可自行探索其他编

7、程方法。2022年年4月月20日星期三日星期三10MATLAB程序代码（ex102.m） clear all; subplot(1,3,1) fplot(-2:2)*x.3+x.2+x,-2 2) %绘制变量a=-2,-1,0,1,2时的图形 grid,axis(equal),axis(-2 2 -4 4) %显示网格，纵横坐标轴保持一致图 %显示区间 subplot(1,3,2)2022年年4月月20日星期三日星期三11 fplot(x.3+(-2:2)*x.2+x,-2 2) %绘制变量b=-2,-1,0,1,2时的图形 grid,axis(equal),axis(-2 2 -4 4) s

8、ubplot(1,3,3) fplot(x3+x.2+(-2:2)*x,-2 2) %绘制变量c=-2,-1,0,1,2时的图形 grid,axis(equal),axis(-2 2 -4 4) gtext(a=2),gtext(a=-2),gtext(a=0) %用gtext命令在其他两个子图上标注字符2022年年4月月20日星期三日星期三12程序运行结果 2022年年4月月20日星期三日星期三1310.3极限判断极限判断 用MATLAB语言来表达推理过程是比较困难的，它必须与实际的数值联系起来，比如无法用无穷小和高阶无穷小的概念，只能用10-10，10-20等数值来表示。极限的定义恰恰是用

9、了和等数值的概念，因此不难用程序表达。2022年年4月月20日星期三日星期三14建模 用函数极限的定义，对于函数y=f(x)，当任意给定一个正数时，有一个对应的正数存在，使得 0时，A-f() 则A就是f(x)在时的极限，如果找不到这样的，A就不是它的极限。在考虑极限时，我们要找到左极限和右极限才能说极限存在，因此要检验 -。xxc xxcxxcxxcxxc2022年年4月月20日星期三日星期三15MATLAB程序代码（ex103.m） clear all; disp(A是否是f(xc)的极限？) A=input(A=(例如A=1); %输入极限值 xc=input(xc=(例如xc=0);

10、%输入对应的自变量值fxc=input(f(x)的表达式为,(例如sin(x)/x),s); %输入函数表达式 flag=1; while flag=1 2022年年4月月20日星期三日星期三16 epslon=input(任给一个小的数epsilon=); %任意给出 delta=1;n=1; %初始化 while abs(delta)eps %找，如果常数eps则失败eps= x=xc-delta; y1=abs(A-eval(fxc); %检验左极限是否存在 x=xc+delta; y2=abs(A-eval(fxc); %检验右极限是否存在 if y1epslon & y2ep

11、slon %如果左极限和右极限都存在则跳出 n=0;break 2022年年4月月20日星期三日星期三17 end delta=delta/2; %减小从新查找 end if n=1 disp(极限不正确或者不存在),delta,break,end disp(delta=),delta %找到了 disp(极限可能存在) flag=input(再试一个epslon吗?再试按1,不试按0或者任意数字键) %是否再试 end 2022年年4月月20日星期三日星期三18104曲面方程 二次曲面的方程如下 求讨论参数a，b，c，d对其形状的影响，并画出其图形。dczbyax2222222022年年4月

12、月20日星期三日星期三19建模 本题的数学模型很清楚，关键在于如何作出三维曲面图形。 特别要注意在给定了x，y求z时，若有开方运算，一是会出现虚数，二是对实数也有正负两个解。 为了使虚数不出现在绘图中，采用了一种技巧，那就是把虚数都换成非数（NaN）。（注，在这里a，b，c三个取值可以是虚数）2022年年4月月20日星期三日星期三20MATLAB程序代码（ex104.m） clear all; a=input(a=);b=input(b=);c=input(c=); d=input(d=);N=input(N=); %输入参数，N为网格线数目 xgrid=linspace(-abs(a*sqr

13、t(d),abs(a*sqrt(d),N); %建立x坐标 ygrid=linspace(-abs(b*sqrt(d),abs(b*sqrt(d),N); %建立y坐标 x,y=meshgrid(xgrid,ygrid); %确定NXN个点的x,y网格坐标 z=c*sqrt(d-y.*y/b/b-x.*x/a/a); z1=real(z); %取z的实部 for k=2:N-1 %以下几行的作用是取消z中的虚数点 2022年年4月月20日星期三日星期三21for j=2:N-1 if imag(z(k,j)=0 z1(k,j)=0;end if all(imag(z(k-1:k+1,j-1:j

14、+1)=0 z1(k,j)=NaN; end endendsubplot(2,2,1)surf(x,y,z1),hold on %画空间曲面z2=-z1;surf(x,y,z2); %画负半面axis(-abs(a*sqrt(d),abs(a*sqrt(d),-abs(b*sqrt(d),abs(b*sqrt(d),. -abs(c),abs(c*sqrt(d)xlabel(x),ylabel(y),zlabel(z)hold off2022年年4月月20日星期三日星期三22程序运行结果 2022年年4月月20日星期三日星期三23105 动态模拟球体表面切片穿过矩形体积的过程 本程序主要利用s

15、lice函数的绘制立体切片图的功能来做一个简单的动态图形，说明如何绘制动态图形2022年年4月月20日星期三日星期三24建模 首先要创建一个边界，同时创建一个球体，我们可以利用sphere函数来完成。然后通过擦除模式来绘制动态图形。 2022年年4月月20日星期三日星期三25MATLAB程序代码 clear all;x,y,z=meshgrid(-2:.2:2,-2:.25:2,-2:.16:2); %创建体积边界数据v=x.*exp(-x.2-y.2-z.2); %计算函数值xsp,ysp,zsp=sphere; %创建球体数据sphere函数绘制单位球体for i=-3:.2:3 hsp=

16、surf(xsp+2,ysp,zsp); rotate(hsp,1 0 0,90) %使球体沿x轴旋转90度 xd=get(hsp,XDATA); %获取x轴数据 yd=get(hsp,YDATA); zd=get(hsp,ZDATA); 2022年年4月月20日星期三日星期三26 delete(hsp) hold on slice(x,y,z,v,-2,2,2,-2) %绘制体积边界 hslicer=slice(x,y,z,v,xd,yd,zd); %动态绘制球面切片 axis tight xlim(-3 3) %限定x的范围 view(-10,35) %设置视点 drawnow delet

17、e(hslicer) hold off end2022年年4月月20日星期三日星期三27运行结果 2022年年4月月20日星期三日星期三28106 模拟抛物运动模拟抛物运动 设目标相对于射点的高度为，给定初速，试计算物体在真空中飞行的时间，距离以及最高点。最后动态的描述出小球的运动轨迹。2022年年4月月20日星期三日星期三29建模 无阻力抛物线的飞行是中学物理就解决了的问题，本题的不同点是目标和射点不在同一高度上，用MATLAB可使整个计算和绘图过程自动化。 其好处是可快速地计算物体在不同初速度和射角下飞行时间和距离。 关键在求落点时间时，需要解决一个二次线形代数方程。2022年年4月月20

18、日星期三日星期三30 由 解出t，它就是落点时间。会有两个解，我们只取其中一个有效解。再求 对于最高点，则如果则就是最高点，如果则最高点为，由于到达最高点时其垂直方向速度为0，所以我们有021sin200gttvyyfftvx00maxcosgvt001sin2022年年4月月20日星期三日星期三31 则最高点有 对于模拟运动轨迹则采用擦除方式，由运动轨迹方程即可绘出。gvygttvyyff2sin21sin202021100max2022年年4月月20日星期三日星期三32MATLAB程序代码 clear all;y0=0;x0=0; %初始位置vMag=input(输入初始速度（m/s）：)

19、; %输入初始速度的大小vDir=input(输入初速方向（度）（输入值大于-90小于90）：); %输入初始速度的方向yf=input(输入目标高度（m）：); vx0=vMag*cos(vDir*(pi/180); %计算x,y方向的初始速度vy0=vMag*sin(vDir*(pi/180);wy=-9.81;wx=0; %重力加速度tf=roots(wy/2,vy0,yf); %解二次线形代数方程，计算落点时间tftf=max(tf); %去除其中的错误的t=0:.1:tf;t=t tf; %添加上最后一个时间点 y=yf+vy0*t+wy*t.2/2; %计算轨迹x=y0+vx0*t

20、+wx*t.2/2;xf=max(x);end2022年年4月月20日星期三日星期三33if vDir-90 %计算飞行距离和最高点 ymax=yf;else ymax=yf+vy0*vy0/2/(-wy); endfprintf(飞行时间为%g秒n,tf);fprintf(最远飞行距离为%g米n,xf);fprintf(飞行最高点为%g米n,ymax);plot(x,y) %绘制抛物线轨迹axis(0 xf+5 0 ymax+5)h1=input(是否动态模拟抛物线运动轨迹？(1/0)（1为是，0为否）);if h1=1 hold on x1=0;y1=yf; h=plot(x1,y1,.)

21、; axis(0 xf+5 0 ymax+5) set(h,Erasemode,xor,MarkerSize,18) for t1=0:0.001:tf y1=yf+vy0*t1+wy*t1.2/2; x1=y0+vx0*t1+wx*t1.2/2; set(h,XDATA,x1,YDATA,y1) drawnow end hold off2022年年4月月20日星期三日星期三34程序运行结果 2022年年4月月20日星期三日星期三35107 弹性碰撞问题弹性碰撞问题 质量为m的小球以速度正面撞击质量为M的静止小球，假设碰撞完全弹性的，即没有能量损失，求碰撞后两小球的速度，及它们与两球质量比K=

22、M/m的关系。2022年年4月月20日星期三日星期三36建模 设碰撞后两小球速度都与v0相同，小球m的速度为v1，小球M的速度为v2，列出动量守恒和能量守恒方程，引入质量比K=M/m，和无纲量速度u1=v1/v0，u2=v2/v0后有 动量守恒 mv0=mv1+Mv2 (1)2022年年4月月20日星期三日星期三372022年年4月月20日星期三日星期三38MATLAB程序代码 clear all; K=logspace(-1,1,11); %设自变量数组K，从K=0.110，按等比取11个点 for i=1:length(K) %对各个K循环计算 ur1=roots(K(i)+1),-2,(

23、1-K(i); %二次方程有两个解 u1(i)=ur1(abs(ur1-1)0.01); %去掉在1临近的庸解 end u2=(1-u1)./K; %用（5）式求u2，用元素群运算 em=1-u1.*u1; %主动球损失的能量（相对值） K;u1;u2;em %显示输出数据 semilogx(K,u1;u2;em);grid2022年年4月月20日星期三日星期三39程序运行结果 数字结果为（省略了几行） K =/ =/ (m小球能量损失) 0.1000 0.8182 1.8182 0.3306 0.2512 0.5985 1.5985 0.6418 1.0000 0 1.0000 1.0000

24、 1.5849 -0.2263 0.7737 0.9488 3.9811 -0.5985 0.4015 0.6418 6.3096 -0.7264 0.2736 0.4724 10.0000 -0.8182 0.1818 0.33062022年年4月月20日星期三日星期三402022年年4月月20日星期三日星期三41108 气体分子运动规律 利用气体分子运动的麦克斯韦速度分布规律，求27度C下氧气分子运动的速度分布曲线，并求速度在400600m/s范围内的分子所占的比例，讨论温度T及分子量度速度分布曲线的影响。2022年年4月月20日星期三日星期三42建模 麦克斯韦速度分布律为 本例将说明如何

25、从复杂数学公式绘制曲线，并研究单个参数的影响。kTmvevkTmf22232)2(42022年年4月月20日星期三日星期三43MATLAB程序 R=8.31; %气体常数k=1.381*10(-23); %波尔茨曼常数NA=6.022*1023; %阿伏加特罗常数T=input(输入气体的绝对温度T（单位K）：);mu=input(输入气体的分子量mu：); m=mu.*1e-3/NA; %分子质量v=0:1500; %设置速度范围y=4*pi*(m/(2*pi*k*T).(3/2).*exp(-m*v.2/(2*k*T).*v.*v; %麦克斯韦分布律plot(v,y),hold on;v1

26、=input(给定速度范围（例如400：600）);if isempty(v1) %如果v1为空则跳出程序 break;endy1=4*pi*(m/(2*pi*k*T).(3/2).*exp(-m*v1.2/(2*k*T).*v1.*v1;fill(v1,v1(end),v1(1),y1,0,0,r) %绘出该范围fprintf(在该范围内的气体分子概率为%g,trapz(y1);2022年年4月月20日星期三日星期三44程序运行结果 2022年年4月月20日星期三日星期三45109 平行电流环之间截面磁场分布的计算 一对相同的共轴载流圆线圈，当它们的间距正好等于线圈半径时，称之为亥姆霍兹线圈

27、。 计算表明，亥姆霍兹线圈轴线附近的磁场的大小分布十分均匀，而且都沿x方向，本例要对这一论断进行验证。2022年年4月月20日星期三日星期三46建模 本题的计算模型与上例相同，只是把观测区域取在两线圈之间的小范围内。 B生成的线圈左边的磁场就等于A线圈的左边磁场，因此，A，B两线圈在中间部分的合成磁场等于A线圈的右磁场与左磁场平移Rh后的和。 因此，只要观测A线圈的左右区间x=-Rh,Rh内的磁场就够了。2022年年4月月20日星期三日星期三47MATLAB程序 clear all;mu0=4*pi*1e-7; %初始化（给定环半径，电流，图形）I0=5.0;Rh=1; %真空导磁率（T*m/

28、A）C0=mu0/(4*pi)*I0;%归并常数 NGx=21;NGy=21; %设定观测点网格数x=linspace(-Rh,Rh,NGx); %设定观测点范围及数组y=linspace(-Rh,Rh,NGy); Nh=20; %电流环分段theta0=linspace(0,2*pi,Nh+1); %环的圆周角分段theta1=theta0(1:Nh); y1=Rh*cos(theta1); %环各段向量的起点坐标y1,z1z1=Rh*sin(theta1);theta2=theta0(2:Nh+1); y2=Rh*cos(theta2); %环各段向量的终点坐标y2,z2z2=Rh*sin

29、(theta2);2022年年4月月20日星期三日星期三48dlx=0;dly=y2-y1;dlz=z2-z1; %计算环各段向量dl的三个长度分量xc=0;yc=(y1+y2)/2;zc=(z1+z2)/2;%计算环各段向量中点的三个坐标分量%循环计算各网格点上的B(x,y)值for i=1:NGy for j=1:NGx rx=x(j)-xc;ry=y(i)-yc;rz=0-zc; %观测点在z=0平面上 r3=sqrt(rx.2+ry.2+rz.2).3; %计算 dlXr_x=dly.*rz-dlz.*ry; %计算叉乘积dlXr的x和y的分量 dlXr_y=dlz.*rx-dlx.*

30、rz; Bx(i,j)=sum(C0*dlXr_x./r3);%把环各段产生的磁场分量累加 By(i,j)=sum(C0*dlXr_y./r3); endendBax=Bx(:,11:21)+Bx(:,1:11);Bay=By(:,11:21)+By(:,1:11);subplot(1,2,1);mesh(x(11:21),y,Bax);xlabel(x);ylabel(y); %画出其B分布三维图subplot(1,2,2); plot(y,Bax),grid,xlabel(y);ylabel(Bx);2022年年4月月20日星期三日星期三49运行结果 可以看出，在亥姆霍兹线圈的两个线圈之间

31、的轴线附近，有相当大的一个区域内，x方向的磁场强度Bx是十分均匀的。在区域内的y方向的磁场强度By近似为零，这留给读者自行证明。为了看出这个区域的大小和形状，可以用多种方法，这里我们用图形的方式，在MATLAB中输入如下代码： q=abs(Bax-Bax(11,6)/Bax(11,6)0.02; i,j=find(q);i=i-11; j=j-6;plot(j,i,.) axis(-6,6,-11,11)2022年年4月月20日星期三日星期三502022年年4月月20日星期三日星期三51图10-10 磁场相对误差点2022年年4月月20日星期三日星期三521010 两点（双缝）光干涉图案 单色

32、光通过两个窄缝射向屏幕，相当于位置不同的两个同频同相光源向屏幕照射的叠合，由于到达屏幕各点的距离（光程）不同引起相位差 叠合的结果是在右的点加强，在有的点抵消，造成干涉现象。 考虑到纯粹的单色光不易获得，通常都有一定的光谱宽度，这种光的非单色性对光的干涉会产生何种效应，要求用MATLAB计算并仿真这一问题。2022年年4月月20日星期三日星期三53建模 考虑到两个相干光源到屏幕上任意点的距离差引起的相位差 则光程差为221)2(zdyLs222)2(zdyLs21LLL2022年年4月月20日星期三日星期三542022年年4月月20日星期三日星期三55MATLAB程序 clear all;%初

33、始化，输入光波长，光缝距离，光栅到屏幕的距离Lambda=500*1e-9;d=2*1e-3;z=1; yMax=5*Lambda*z/d;xs=yMax; %设定图案的x,y向范围Ny=101;ys=linspace(-yMax,yMax,Ny); %y方向分成101点flag=input(是否考虑光的非单色性？（考虑输入1，不考虑输入2）：);for i=1:Ny %对屏上全部点进行循环计算 %计算第一个和第二个光源到屏幕各点的距离 L1=sqrt(ys(i)-d/2).2+z2); L2=sqrt(ys(i)+d/2).2+z2); switch flag case 1 %考虑光的非单色

34、性 N1=11;dL=linspace(-0.1,0.1,N1); %设光谱相对宽度10% 2022年年4月月20日星期三日星期三56Lambda1=Lambda*(1+dL); %分11根谱线，波长为一个数组 Phi1=2*pi*(L2-L1)./Lambda1;%从距离差计算各波长的相位差 B(i,:)=sum(4*cos(Phi1/2).2)/N1; %叠加各波长影响计算光强 case 2 %不考虑光的非单色性 Phi=2*pi*(L2-L1)/Lambda; %从距离差计算相位差 B(i,:)=4*cos(Phi/2).2; %计算该店光强（设两束光强相同） endendNCLevel

35、s=255; %确定用的灰度等级为255%定标；使最大光强（4.0）对应最大灰度级（白色）Br=(B/4.0)*NCLevels; subplot(1,2,1),image(xs,ys,Br); %画图像colormap(gray(NCLevels); %用灰度级颜色图subplot(1,2,2),plot(B(:),ys) %画出沿y向的光强变化曲线2022年年4月月20日星期三日星期三57程序运行结果 运行程序时分别输入1，2（输入1表示考虑光的非单色性，输入2表示不考虑光的非单色性）得到如图10-12，图10-13，可以看出，广的非单色性导致干涉现象的减弱。光谱很宽的光将不能形成干涉。2

36、022年年4月月20日星期三日星期三58图10-13双缝干涉条纹及光强分布 单色光时的结果非单色光时的结果图10-12双缝干涉条纹及光强分布非单色光时的结果2022年年4月月20日星期三日星期三591011弹簧震动系统模型 实时动画制作弹簧震动系统模型。当质量为M的物体连接在刚度为K的弹簧上，物体初始距平衡位置有一个位移。如果不考虑物体和平面的摩擦，则物体将做相对平衡位置的理想简谐运动，其运动方程为： 其解为： 其中x0=xt=0。参数选择K=M，x0=5。0 KxxMtMKxxcos02022年年4月月20日星期三日星期三60MATLAB程序（ex. 1010.m） offset=4; %滑

37、块的大小animinit(onecart Animation);onecart=findobj(Type,figure,Name,onecart Animation);axis(-10 20 -7 7);hold on;xySpr1=. 0 0; .4 0; .8 0.65; 1.6 -0.65; 2.4 0.65; 3.2 -0.65; 3.6 0; 4.0 0; %弹簧的横纵坐标xyBx11=. 0 1.2; 0 -1.2; 0 0; %滑块的数据2022年年4月月20日星期三日星期三61xyBx21=. 0 1.2; 2 1.2; 2 -1.2; 0 -1.2; 0 1.2; %滑块的数

38、据xBx11=xyBx11(:,1);yBx11=xyBx11(:,2);xBx21=xyBx21(:,1);yBx21=xyBx21(:,2);xSpr1=xySpr1(:,1);ySpr1=xySpr1(:,2);x=xBx11;xSpr1;xBx21(:,1)+offset;y=yBx11;ySpr1;yBx21;%画出滑动物体下的水平面plot(-10 20,-1.4 -1.4,blue,. -10:19;-9:20,-2 -1.4,blue,LineWidth,2);hold on;%画出弹簧和滑块2022年年4月月20日星期三日星期三62hndl=plot(x,y,y,color,

39、blue,EraseMode,background,LineWidth,3);set(gca,UserData,hndl);set(gcf,color,yellow);t=0;dt=0.005;u0=5;while t20 if any(get(0,Children)=onecart) if strcmp(get(onecart,Name),onecart Animation) set(0,currentfigure,onecart); t=t+dt; u=u0+5*cos(t); u=u+offset; distance=u;2022年年4月月20日星期三日星期三63 hndl=get(gc

40、a,UserData); x=xBx11; xSpr1/4*distance; xBx21+distance; set(hndl,xData,x); drawnow; if(t=0) for(p=0:0.01:10*p1); end end end endend2022年年4月月20日星期三日星期三64程序的运行结果 2022年年4月月20日星期三日星期三652022年年4月月20日星期三日星期三661012 运算器放大电路运算器放大电路 运算器放大电路如图10-16所示，试分析放大器开环增益和频率响应对整个电路闭环频率响应的影响，并绘出曲线。图10-16 运算放大器等效电路2022年年4月月

41、20日星期三日星期三67建模 设运算放大器的开环增益为A，它是频率的函数，则在图示的连接方法下，闭环输出与输入电压之比为 在增益A很大时，分母上的第二项可以忽略不计，因而得出理想运放的闭环传递函数AZZZZUUHio/ )/1 (1/1212)()()()()(12sZsZsUsUsHio2022年年4月月20日星期三日星期三68 通常，运算放大器的开环传递函数中包括三个实极点，即 其中-w2f)&(ww2f); %求出低通滤波器带宽内的下标w2=w(n2); %取出中段频率Y2=F(n2); %取出中段频谱figure(1),subpot(1, 2, 1),plot(w2,abs(Y

42、2),linewidth,1.5),grid %画出滤波后的频谱y2=Y2*exp(j*w2*t)/pi*dw; %对中段频谱求傅里叶逆变换2022年年4月月20日星期三日星期三100subplot(1,2,2),plot(t,f, t,y2,linewidth,1.5),grid %画出原波形及滤波后的波形%（2）1017b%滤波器2 ：截止频率为wc=10的三阶巴特沃斯低通滤波器转移函数为q611g%调用例6.11程序，设定数据为N=256,wf=20,Nf=128H=freqs(500,1, 20,200,1000,w);%如用信号处理工具箱，语句为%b,a=butter(wn,阶数)；

43、H=freqs（b,a,w）,figure(2),subplot(1,2,1),plot(w,abs(Y3),linewidth,1.5),grid %画出滤波后的频谱y3=Y3*exp(j*w*t)/pi*dw;%对中段频谱求傅里叶逆变换subplot(1,2,2)plot(t,f,t,y3,linewidth,1.5),grid %画出原波形及滤波后的波形2022年年4月月20日星期三日星期三101 （3）程序运行结果。 运行ex1017a.m，设定理想滤波器带宽为10时所得之图。 运行ex1017b.m，通过三阶巴特沃斯低通滤波器后所得之图。2022年年4月月20日星期三日星期三1021

44、0.18 二阶巴特沃斯低通数字滤波器的频率响应 【例10.18】二阶巴特沃斯低通滤波器的系统函数（传递函数）为： 求其频率响应并画图（）。2221( )(22)(22)zzH zz 2022年年4月月20日星期三日星期三103 解 （1）建模： 离散系统的频率响应函数为（这里，ts为抽样周期）： 其幅频特性为 ，相频特性为 ，按下面定义编程。信号处理工具箱中的freqz()函数与此程序有同样功能，不过在信号处理工具箱中，滤波器的系统函数是用z-1作为自变量来描述的，本题应为：e(e )( )jjzHH z()jH e()jH e121212()(22)(22)zzH zz2022年年4月月20

45、日星期三日星期三104 （2）MATLAB程序（ex1018.m）如下：%二阶巴特沃斯低通滤波器的离散系统函数为：%H(z)=(z2+2z+1)/(2+sqrt(2)z2+(2-sqrt(2);%求其频率响应可将z用exp(i*w)代入b=1,2,1;a=2+sqrt(2),0,2-sqrt(2); %给出滤波器分子分母系数N=input(取频率数组的点数N= );w=0:N-1*pi/N; %给出0到之间的频率数组H=polyval(b,exp(i*w)./polyval(a,exp(i*w); %求频率响应figure(1)subplot(2,1,1),plot(w,abs(H),grid

46、subplot(2,1,2),plot(w,unwrap(angle(H),gridfigure(2) %在线性坐标内画频率特性subplot(2,1,1),semilogx(w, 20*log10(abs(H),gridsubplot(2,1,2),semilogx(w, unwrap(angle(H),grid %在对数坐标内画频率特性2022年年4月月20日星期三日星期三105 （3）程序运行结果如图10-27所示。2022年年4月月20日星期三日星期三10610.19 时域采样频率与频谱混叠 【例10.19】函数如下： 计算 并图示及其幅频特性函数0( )esin( )atax tAt

47、u t444.128A 50 2a 050 2( )ax t()aXj2022年年4月月20日星期三日星期三107 解 （1）建模：对 进行等间隔采样，得到 ， 为采样周期。如果 ，则： 由以上关系式可见，采样信号的频谱函数是原模拟信号频谱函数的周期延拓，延拓周期为2/T。 ( )ax t( )()ax nxnT1/sTf()( )aaXjFT x t 12(e)( )( )e()njTjTankkXFT x nx nXjTT2022年年4月月20日星期三日星期三108 （2）MATLAB程序（ex1019.m）如下：%时域采样及其频谱clear;close all;fs=10000;fs1=

48、1000;fs2=400;fs3=200; %设置四种采样频率t=0:1/fs:0.1; %采集信号长度为0.1sA=444.128;a=50*sqrt(2)*pi;b=a;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); %近似模拟信号频谱T1=1/fs1;t1=0:T1:0.1; %采集信号长度为0.1s2022年年4月月20日星期三日星期三109x1=A*exp(-a.*t1).*sin(b*t1); %1kHz采样序列x1(n)X1=dtft(x1,2*pi*k/512); %x1(n)的512点dtftT2=1/fs2;t2=0:T2:0.1; %采集信号长度为0.1sx2=A*exp(-a.*t2).*sin(b.*t2); %400Hz采样序列x2(n)X2=dtft(x2,2*pi*k/512); %x2(n)的512点dtftT3=1/fs3;t3=0:T3:0.1; %采集信号长度为0.1sx3=A*exp(-a.*t3).*sin(b.*t3); %200Hz采样序列x3(