基于MATLAB的卷积积分 计算卷积的通用函数

中国西部科技2008年6月(中旬)第07卷第17期总第142期基于mrⅡ。AB的卷积积分金波(长江大学电子信息学院，湖北荆州43402；3)摘要：本文介绍了卷积积分的数值计算、符号计算和GuI程序设计的方法和技，亍。给出了卷积积分的数值计算通用函数CSCONV()、符号计算的通用函数CSCONVS()，通过实例计算，验证了这些函数的有效性。展示了MATLAB在卷积积分计算方面的优点。关键词：卷积积分；MATLAB；程序设计l引言在线性时不变连续系统中，利用系统的冲激响应和叠加原理来求系统对任意激励信号作用时的零状态响应，这就是卷积方法的原理。因此，在时域内，卷积积分是求解线性非时变系统零状态响应的重要方法。特别是激励信号为时限信号时尤其如此。卷积积分的计算比较复杂，是信号与系统分析中的重点和难点，特别适合用计算机来计算。以往的卷积积分多用fortran、C、VB等语言编程，不仅编程烦琐，而且可视性差。用MATLAB来计算卷积积分问题要比用C、FORTRAN等语言完成相同的事情简捷得多。ICATLAB软件是由美国Mathworks公司推出的用于数值计算和图形处理的科学计算系统环境。MATLAB是英文以TrixLABoratory(矩阵实验室)的缩写。它包括高效的数值计算、矩阵运算、信号处理和图形生成等功能。在blATLAB环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件管理等各项操作。MATLAB的出现给信号与系统的分析提供了很大的方便。本文结合卷积积分给出了如何应用MATLABI程计算语言，用数值计算和符号计算来解决卷积积分的计算问题。2卷积积分的数值计算卷积积分计算实际上可用信号的分段求和来实现，即：．一=，(f)2爪，)．f2(t)2JiA(x)12(t—t胁=l。iⅢra∑I,(kh)12(t一尬)‘△㈣如果我们只求当r=nA(n为整数)时f(t)的值，(以)，则由上式可得：：，(H△)=A艺州尬)，2【∽一★)△】t=—z当时间间隔足够小时，，(”△)就是f(t)的数值近似。Matlab的conv(X，h)函数可以用来计算卷积积分的数值解。为此编写了计算卷积积分的通用函数如下。functionx=CSCONV(F1，tl_s，tl_e，F2，t2一S，t2一e)；％计算有限长连续信号的卷积：％F1，F2为两信号的字符串，tl—S，t2一S为两信号的起始点，tl—e，t2一e为两信号的终止点；％同时显示f1(t)。f2(t)和f(t)=fl*f2一---个波形，并且三个波形自动显示在图形中间。％例：fl=’2*rectpuls(t，4)’；f2=’3*rectpuls(t-0．5，3)’；figure(1)％CSCONV(fl，一2，2，f2，一1，2)％t—s=tl—s+t2一S：t_e=tl—e+t2一e；t0=min([tl—S，t2一S，t—S])一1；tl=max“tl—e，t2一e，t_e])+l；dt=O．005：t=tO：dt：t1：L=length(t)：tp=[2*t(1)：dt：2*t(L)]；fl=eval(F1)：f2=eval(F2)：y=dt*conv(fl，f2)；y_max=max(Y)：y_min=min(Y)：dy=(y_max—y_min)／10：fl—max=max(f1)：fl—min=min(f1)：dfl=(fl—max—fl_min)／10．f2一max=max(f2)：f2一min=mifl(f2)：df2=(f2一max—f2一min)／10：subplot(3，1，1)，plot(t，fl，’1inewidth’，2)：title(’信号fl(t)的波形’，’color’，’b’，’FontSize’，8)axis([t(1)t(L)flmin-dflfl_max+dfl])：grid，set(gca，’FontSize’，8)subplot(3，1，2)，plot(t，f2，’linewidth’，2)：title(’信号f2(t)的波形’，’color’，’b’，’FontSize’，8)AxiS([t(1)t(L)f2min-dflf2一max+dfl])：grid，set(gca，’FontSize’，8)subplot(3，1，3)，plot(tp，Y，’linewidth’，2)：title(’卷积f(t)=fl*f2的波形’，’color’，’b’，’FontSize’，8)ify_min>=OyO=(y_max—y_min)／2：elseyO=O；end收稿日期：2008-05-10修回日期：2008—06—06作者简介：金波l1955一)．男．副教授．现从事电路、信号与系统、信号处理等方面的教学和研究工作。40万方数据axis([t(1)，t(L)，y_min—dy，y_max+dy])；grid：set(gca，’Ytick’，[y_min，yO，y_max]，’FontSize’，8)调用以上函数可以方便地计算画出两个任意有限长波形和卷积波形。用下面的例子说明。例1用Matlab画出图中的信号的卷积波形。l，l‘f)厂]2。r]～-，l【"厂]～一lOl(a)¨’2卅～-2∥O2_-2(b)图1Matlab程序如下：fl=’2*rectpuls(t+1．5)+2*rectpuls(t-1．5)’：f2=’2*rectpuls(t-2。4)’；figure(1)CSCONV(fl，一2，2，f2，0，4)fl=’2*rectpuls(t，2)’：f2=’t．*rectpuls(t。4)’；figure(2)CSCONV(fl，一1，l，f2，一2，2)运行程序后显示的图形如图2所示。(a)(B)图23卷积积分的符号计算Matlab不仅有数值计算功能，还有强大的符号计算功能。即有推理功能町得解析表达式。卷积积分用Matlab的积分int0函数，可以计算不定积分和定积分。调用格式为Int(f)，int(f，a，b)其中，f为被积函数，a，b为积分的上下限。为此编写了计算卷积积分解析解的通用函数如下：厂]墅壅皇里functiony=CSCONVS(f，h，t_s，t_e，a，b)；％用计算卷积积分的解析解％f激励信号，含有阶跃函数Heaviside(t)；％h冲激响应，卷积时要反折的信号，不可含阶跃函数，默认起始点为0；％t_s，t_e为系统零状态响应y=f*h波形的起始点和终止点，％a，b为卷积积分的上，下限。例：计算exp(一2t)u(t+1)与u(t-3)卷积。％symst％h=sym(’exp(一2．t+2)’)：％f=sym(’HeaviSide(t-2)’)：％CSCONVS(f，h，1，5，2，t)％symsttaoftao=subs(f。t。tao)；ht—tao=subs(h，t，t-tao)；y=simple(int(ftao*ht—tao，tao，a，b))；t=t—S：0．02：t—e；yt=subs(Y)：y_max--max(yt)：y_min=min(yt)：dy=(y_max—y_min)／i0：plot(t，yt，’linewidth’，2)，axis([t_st_ey_min—dyy_max+dy])；％axis([一23011])；xlabel(’t(sea)’)；grid；ify_min>=OyO=(y_max—Y—min)／2；elseyO=O：endset(gca，’Ytick’，[y_min，yO，y_max]，’FontSize’，8)titie(’卷积的波形’)；Disp(’零状态响应’)，Y现用实例说明其计算方法。例2线性非时变系统的输入信号f(t)和冲激响应h(t)由下列各式给出，试求系统的零状态响应％(，)。(a)，(f)=e“”I；o)一￡(r一2)】，^(，)=e一’￡“)；(B)，(f)=e-2,￡O+3)，^(r)=e-2t￡(1一1)。解(a)系统的零状态响应为：比(，)=，(f)·^(f)=J-ie-“f”【￡o)一60-2)]e-(t-．,t)E(卜T)出=j二e4“e1。l￡o碡(f—t)dr—j二e““e⋯1’e(t一2)e(t—t)dr=【e“【e““出】e(f)-∥p“出)E(f一2)=2e4“一e-t)s(，)一2e“，一e-O-I))e(t-2)(b)系统的零状态响应为：)，。(f)=，(f)·向(f)=12e一“eo+3)p一2‘’1】￡(f—T—1)dr=【P一2‘It，-I(打k(f+2)：(r+2)口一2’￡(，+2)为了用CSCONVS()函数计算，利用卷积的时移性质将被积函数要变换如下：41万方数据中国西部科技2008年6月(中旬)第07卷第17期总第142期用MATLAB并调用CSCONVS()函数的计算程序为：y口(f)=e．2’eO+3)·e一2’￡O—1)=e-2(t-I)￡O+2)’e一2‘⋯’e(f)symsth--sym(’exp(一t)’)；f=sym(’exp(一t／2)木(HeaviSide(t)一Heaviside(t一2))’)：figure(I)yI=CSCONVS(f，h，一0．5，5，0，t)H=sym(’exp(一2*t一2)’)；f=sym(’exp(-2．t+2)：Ic(Heaviside(t+2))’)；figure(2)y2：CSCONVS(f，h，一2．5，3，一2，t)disp(’零状态响应’)，yl，y2，运行程序后显示的图形如图3所示。在命令窗口显示其卷积的解析式如下：零状态响应。yl=2*Heaviside(t)*exp(-1／2*t)一2*HeaviSide(t)*exp(一t)一2*HeaviSide(t-2)*exp(一I／2*t)+2*HeaviSide(t一2)*exp(卜t)y2=Heaviside(t+2)木(t+2)*exp(一2*t)其中，Heaviside(t)表示阶跃函数g(t)，所以，y．(，)=2(e1“一e-i)e(f)一2(e“”一e-O-i))60一2)y2(，)=(f+2)P4’￡(，+2)与理论计算结果一致。(a)y。(f)的波形(b)Y2(t)的波形田34结束语由于MATLAB的强大数值计算和符号计算功能，使卷积积分的计算和可视化变为可能。通过以上实例分析可以看到，应用MATLAB不仅可以对卷积积分进行数值计算，画出波形。而且可以进行符号计算得到解析表达式，画出波形。本文提出的用MATLAB来解决问题的一些方法和技巧出可推广到信号与系统分析和其它领域。参考文献：[1]金波．信号与系统基础[M]．武汉：华中科技大学出版社，2006，8．[2]王宏．^IATLAB6．5及其在信号处理中的应用[M]．北京：清华大学出版社，2004，10．[3]张昱，周绮敏．信号与系统实验教程[M]．北京：人民邮电出版社，2005．面iiiiiiijijab舀西iiiiij：j=j石i6西幻tiiiii—i一⋯⋯⋯一云⋯⋯wI—gI—I⋯III一1～1⋯：1⋯1～1—1Il—l⋯lll～l～。l⋯lIll⋯l⋯⋯～一⋯wlg—ll：l：iK—：：lO—D(上接第57页)(3)反应介质及酸度。①选择盐酸、硫酸、硝酸进行测定，加入硝酸后二硫腙被氧化故不能在此介质中测定。在盐酸、硫酸介质中均能测定。②使用309／L的KBH。及在最佳工作条件下，研究了盐酸、硫酸酸度对反应的影响：在0．22～0．32mol／L盐酸范围内，0．18～0．26mol／L硫酸范围内效果均好，但使用硫酸灵敏度要比盐酸高，因此，我们选择了0．25mol／L硫酸酸度下测定。(d)还原剂的浓度。KBH。浓度对测定有显著影响，为提高KBH．水溶液的稳定性，加入KOH，当KOH浓度为59／L时，试验表明，KBH。浓度在309／L以上出现平台，因此我们选择]'309／L浓度的KBH。。(5)共存离子的干扰和消除。本方法在Cd1ug／L浓度下，对生物样品中常见离子进行干扰实验，结果表明：105倍的Ca、Mg、K、Na、，104倍的Fe、Al，103倍的Zn、Mn、Cr、Mo、Ni、Co，100倍的As、Sb、Bi、Hg、Sn、Se，80倍的Pb不干扰测定，40倍的Cu对测定有干扰。因为生物样品中Cu是常存在的离子，且浓度也超过干扰倍数，故去除Cu的干扰是必要的。我们研究了许多种去除干扰的方法，最后采用0．59／L二硫腙一四氯化碳5．OmL，在0．25mol／L硫酸酸度下可去除500倍Cu的干扰，而萃取后不用分离即可进行镉的测定。(6)线性范围、检出线。按实验确定的最佳条件进行测定绘制工作曲线，镉浓度在o～50ug／L范围内荧光强度极42限If与其浓度c呈线性关系，直结回归方程：If=71．787C-4．249；相关系数r=O．9997；检出限为0．12IJg／L。(7)准确度、精密度、回收率。同一样品同时处理两份，测定5次，其RSD为4．096。两种物质的分析结果见表2。对不同的实际样品进行标准加入回收试验，结果见表3。表2标准物质分析结果I标准物质标准僵(“g／L)l测定僵(4l，L)l始刃GBl-085714．5±0．54．6l甘兰GBW08504o．o∞土0．006o．026表3回收率测定‘吉果l样品名称含量(rig)加入量(nlO测得僵(吨)I回收率(％)l词科l624005821050l大米0433口0335I卯334结论实验表明，原子荧光光度法用于测定生物样品中的镉，是一种可行的方法。该法简便、灵敏度高、准确性好。把研究结果投入到批量品测试生产中，并用标准物质进行监控，提高了测试效率，保证了测试质量，已经取得了很好的经济效益。参考文献：[1]北京科创海光仪器有限公司编写组．原子荧光光度计分析方法[R]．北京：北京科创海光仪器有限公司，2004．[2]岩石矿物分析编写组．岩石矿物分析(第三版)[M]．北京：地质出版社，1991．万方数据基于MATLAB的卷积积分作者：金波作者单位：长江大学电子信息学院,湖北,荆州,434023刊名：中国西部科技英文刊名：SCIENCEANDTECHNOLOGYOFWESTCHINA年，卷(期)：2008,7(17)参考文献(3条)1.张昱;周绮敏信号与系统实验教程20052.王宏MATLAB6.5及其在信号处理中的应用20043.金波信号与系统基础2006本文链接：http://d.g.wanfangdata.com.cn/Periodical_zgxbkj200817019.aspx