因果序列卷积的矢量乘法算法程序设计定稿

因果序列卷积的矢量乘法算法程序设计摘要卷积是信号后期处理中的重要概念。离散信号的卷积计算在文献中有多种解法。我们提出了用矢量乘法计算离散信号卷积和的方法。与传统计算方法相比,该方法避开了繁琐的做图求和过程,使计算更简便、快速。同时，也对几种方法进行了对比和讨论。关键字因果序列;卷积和;矢量乘法程序设计ThevectormultiplicationofthecausalsequenceofconvolutionalgorithmprogrammingAbstract:Convolutionisanimportantconceptin“Signalspost-processing”.Methodsforcalculatingcausalsequenceofconvolutionaregiven.ThemethodtocomputetheconvolutionsumofdiscretesignalsbyusingMathematicaisintroduced.Incomparisonwiththetraditionalcomputation,thismethodavoidsthecomplicatedproceduressuchasdrawing.Finally,comparisonsanddiscussionsofthethesemethodsareprovided。Keywords:causalsequence;convolutionsum;vectormultiplication;program第0页共11页目录1引言11.1离散序列卷积的意义11.2卷积的用途22课题研究方法22.1卷积的计算机实现22.2课题实现工具23序列卷积常用算法及其实现思路23.1定义法23.2矩阵乘法23.3LISTCONVOLVE算法34卷积矢量乘法算法的程序设计35各种卷积算法优劣比较45.1矩阵乘法45.2LISTCONVOLVE算法45.3矢量乘法算法45.4FFT算法55.5列表比较56卷积算法程序应用实例66.1设计滤波器66.2物理实例67结语10致谢10参考文献11第1页共11页1引言在很多人眼中，卷积这个概念是比较神秘且不容易理解的。卷积在数学、物理学、电子工程、信号处理、计算机科学中极为重要。离散信号卷积在电子通信领域的应用非常广泛，也是工程应用的基础。所以快速有效的计算离散序列的卷积，是人们一直很关心的一个问题。1.1离散序列卷积的意义在科学技术飞速发展的今天，现代通信正在蓬勃发展，如何快速准确、方便易行地进行远距离信息传递，这是人类追求的目标。要想实现这一目标就必须掌握离散卷积的运算。卷积是在信号与线性系统的基础上和背景中出的，卷积关系很重要的一种情况就是，在信号与线性系统中或数字信号处理中的卷积定理。因为利用该定理，可以将进行时间域或空间域中的卷积运算等价为频率域的相乘运算，从而利用FFT等快速算法，实现更有效的计算，节省运算时间，提升效率。信号通过信号放大系统就是卷积，因为系统是对信号做了线性加工。进一步推广照相系统是卷积，视觉系统也是卷积。我们看到的物体是物体本身与我们的视觉系统的卷积，因为同一个物体，视力不同的人看到的形象不尽相同。同样我们听到的声音是声源发出的声音与听觉系统的卷积。世界上所有的系统都是卷积，例如脉搏是心跳的卷积，水压的波动时泵水电机运转动力的卷积，只要系统是线性的或近似线性的。所以当我们正确认识了卷积的意义之后就会发现卷积是无处不在。但是卷积的意义建立在这两个前提之下：一是任意的数字信号都可以表示成单位脉冲之线性组合，二是该系统也是线性的。实际上，我们是根据我们需要待处理信号的形式，来设计所谓的系统传递函数，那么这个系统的传递函数和输入信号，在数学上的形式就是所谓的卷积关系。1.2卷积的用途卷积的用途就是通过物理系统对信号进行滤波运算，然后对其信号进行幅度的放大及延时。卷积在信号中主要是时域、频域的相互转换。卷积在其它领域也有很多运用,例如在编码中，有卷积码，就是运用原码中前后序列的码字确定当前的编码输出。卷积在分析系统的零状态响应过程中也有着重要作用,第2页共11页2课题研究方法2.1卷积的计算机实现在计算机上进行卷积运算就是做一个窗口叠加运算，我们一般看的计算机信号处理方面的教材，例如数字图像处理，其中用卷积方的法就实现了图像的积分、微分、锐化、平滑、去噪声等各种运算。2.2课题实现工具随着计算机技术行业的发展,越来越多的计算问题都已经交给计算机来处理。Matlab和Mathematica作为优秀的科学计算的软件,在信号处理与通信、工程计算、图像处理技术等领域均得到了广泛的应用。MATLAB主要用于数据可视化、算法开发、数值计算及数据分析的交互式环境和高级计算语言。Mathematics是一个教育软件，设计给MicrosoftWindows，使用户能够解决数学的科学问题。由微软开发并维护，它主要作为学生和老师的学习工具。Mathematica拥有强大的符号及数值运算能力以及方便实用的绘图功能,所以应用Mathematica总会让人身心舒畅。用Mathematica来计算积分、求和、画图都非常容易,所以可以利用这些特点,做好计算卷积的程序包或函数,使计算过程大大简化,并可以得到精确的数值解或者解析解。第3页共11页3序列卷积常用算法及其实现思路3.1定义法离散信号的卷积和定义为1：(式3.1)()*()()iykfgkfgki式中：f(k)和g(k)为离散时间序列；i为虚设求和变量，卷积的计算结果为另一新的序列。用如下例子说明：已知,求卷积。()()01)nxau()hnu()*()ynxh解答：10*mmnayauu当时n1()ya3.2矩阵乘法用如下例子说明：求序列x=2,3,4和序列y=1,2,9,0,3,5的卷积，可以等效为矩阵和矩阵432相乘。5309210运用mathematics运算，编程如下：x=2,3,4；y=1,2,9,0,3,5；L=8;z=x.1,2,9,0,3,5,0,0,0,1,2,9,0,3,5,0,0,0,1,2,9,0,3,5Printz=,z;输出结果：z=2,7,28,35,42,19,27,203.3ListConvolve算法ListConvolve算法是mathematics内部一个算法，主要是进行循环卷积，基本用法:ListConvolveker,list构造核ker和list的卷积。出现的新概念是核。核即指核心，便是值输入的信号。list自然是单位冲击信号对应的响应了。下面是在mathematics上进行的一个简单编程：x=2,3,4;y=1,2,9,0,3,5;z=ListConvolvex,y,1,-1,0;Printz=,z;输出结果：z=2,7,28,35,42,19,27,20第4页共11页4卷积矢量乘法算法的程序设计运用mathematics设计矢量乘法程序，编写程序如下：x=TableCos*k/171Exp-0.0002k,k,0,9999;y=Table(0.999)k,k,0,99999;t1=Date;m=Lengthx;n=Lengthy;L=m+n-1;z=ConstantArray0,L;s=JoinConstantArray0,m-1,yn;1;-1;Forj=1,jL,j+,s=RotateRights;zj=x.s1;m;t2=Date;z2=ListConvolvex,y,1,-1,0;t3=Date;(*Printx*y=,z*)Printerror=,Maxz-z2;p1=ListPlotz1;15000,PlotRangeAll,PlotStylePointSize0.006,Red;p2=ListPlotz21;15000,PlotRangeAll,PlotStylePointSize0.003,Blue;Showp1,p2Print矢量乘法用时,(t2-t1).tk,secnListConvolvey用时,(t3-t2).tk,sec;第5页共11页5各种卷积算法优劣比较前文已主要介绍了几种求解卷积的算法，但是各种算法从运算效率和占用内存上来说肯定存在一定的差异，在工作和学习中为了选择更合适的，很有必要对各种算法进行讨论。在选择算法时我们必须本着简洁，实用，便捷的原则。5.1矩阵乘法tk=0,0,0,3600,60,1;x=TableCos*k/171Exp-0.0002k,k,0,999;y=Table(0.999)k,k,0,999;t1=Date;n=Lengthx;m=Lengthy;L=n+m-1;y0=PadRighty,L;tem=y0;m=tem;Forj=1,j=km,yk-m+1,0,m,0,Lengthx-1;zn=Tablezk,k,0,Lengthx+Lengthy-2;dt=(Date-t1).0,0,0,3600,60,1;Printdt=,dt;*)t3=Date;ListConvolvex,y,1,-1,0;dt2=(Date-t3).0,0,0,3600,60,1;Printdt2=,dt2;运算结果：dt2=_0.00100015.3矢量乘法算法tk=0,0,0,3600,60,1;x=TableCos*k/171Exp-0.0002k,k,0,9999;y=Table(0.999)k,k,0,9999;t1=Date;第6页共11页m=Lengthx;n=Lengthy;L=m+n-1;z=ConstantArray0,L;s=JoinConstantArray0,m-1,yn;1;-1;L,j+,Forj=1,jL,j+,s=RotateRights;zj=x.s1;m;t2=Date;z2=ListConvolvex,y,1,-1,0;t3=Date;Printerror=,Maxz-z2;p1=ListPlotz1;15000,PlotRangeAll,PlotStylePointSize0.006,Red;p2=ListPlotz21;15000,PlotRangeAll,PlotStylePointSize0.003,Blue;Showp1,p2Print矢量乘法用时,(t2-t1).tk,secnListConvolve用时,(t3-t2).tk,sec;运算结果：矢量乘法用时10.6316081secListConvolve用时0.0040002sec5.4FFT算法x=RandomInteger1,5,10000;xL=Lengthx;y=RandomInteger1,5,20000;yL=Lengthy;L=xL+yL-1;t1=Date;xn=PadRightx,L;yn=PadRighty,L;X=Fourierxn,FourierParameters1,-1;Y=Fourieryn,FourierParameters1,-1;Z=X*Y;zfft=FourierZ,FourierParameters-1,1;Printdtfft=,(Date-t1).0,0,0,3600,60,1;t2=Date;zlc=ListConvolvex,y,1,-1,0;Printdtlc=,(Date-t2).0,0,0,3600,60,1;er=MaxAbszfft-zlc5.5列表比较表5.1算法复杂度比较序列一长度序列二长度计算耗时占用存储空间算法明晰度矩阵乘法算法9999990.1159976明晰ListConvolve算法999999990.00528747明晰矢量乘法算法999999990.004159976明晰FFT算法999999990.01159976明晰结论：矩阵乘法不适合用于长序列，否则会造成mathematics内核自动停止运行，这时就需要关闭其它程序来再一次尝试。在运行短序列时它的效率还是比较高的。正因为矩阵乘法只适合短序第7页共11页列所以教学过程中一般使用这个方法。ListConvolve算法不仅适用于长序列也适用于短序列，而且效率最高，非常节省时间。矢量乘法算法不仅适用于短序列，也适用于长序列，但是长序列时相当费时。FFT算法长序列和短序列都适用，较长序列一般考虑ListConvolve算法和FFT算法。相较之下ListConvolve算法最优，接下来是FFT算法，矢量乘法算法，最后是矩阵乘法。第8页共11页6卷积算法程序应用实例6.1设计滤波器滤波器设计主要是借助MATLAB，设计程序如下：(*确定阶数255；确定通带上限频率；fp=fs/m;b=fir1(255,w);savezmxbcd*)6.2物理实例为了更好的了解矢量乘法算法，我找了一份正常人的心电图作为实例来说明。首先从MATLAB调用滤波器,程序如下：b=ImportD:ProgramFilesMATLABR2012abinzxmfilter.mat/Flatten;记录心电图上的数据，然后调用到mathematics，程序如下：s=ImportC:UsersAdministratorDesktop111.txt,Table/Flatten;yn=-0.2,0.05,0.09,0.1,0.1,0.1,0.2,0.3,0.3,0.2,0.2,0.2,0.29,0.29,0.29,0.21,0.21,0.3,0.2,1.3,-0.2,0.1,0.12,0.12,0.1,0.18,0.2,0.21,0.12,0.11,0.1,0.09,0,-0.01,-0.1,-0.12,-0.15,-0.1,-0.1,-0.2,0.9,-0.5,-0.2,-0.15,-0.2,-0.1,0,0.15,0,-0.02,-0.1,-0.1,-0.08,-0.05,0,0,0,0.1,0,-0.01,1,-0.3,0,0.02,0,0,0.1,0.2,0.12,0.1,0,-0.1,-0.09,-0.09,0,0,0,0.1,0.1,0,1.1;在mathematics上进行矢量乘法进行卷积，程序如下：ynL=Lengthyn;ListPlotyn,JoinedTrue,PlotRangeAll,AspectRatio1/3bL=Lengthb;n=FloorbL/2;Ts=0.04;fs=1/Ts;tsb=Range0,bL-1*Ts;Print滤波器单位冲激响应时域波形;ListPlotTransposetsb,b,JoinedTrue,AxesLabelt,h(t),PlotRangeAll,InterpolationOrder2H=Fourierb,FourierParameters1,-1;HA=AbsH1;n;F=Range0,bL/2-1*fs/bL;Print滤波器幅频特性;ListPlotTransposeF,HA,JoinedTrue,AxesLabelf/Hz,H(f),PlotRangeAllListPlotTransposeF,ArgH1;n,JoinedTrue,AxesLabelf/Hz,(f),PlotRangeAllynL2=FloorynL/2;ts=Range0,ynL-1*0.04;Print待处理信号滤波前波形;ListPlotTransposets,yn,JoinedTrue,AxesLabelt,y(t),PlotRangeAll,InterpolationOrder2Yn=2Fourieryn,FourierParameters-1,-1;Yn1/=2;YA=AbsYn1;ynL2;F2=Range0,ynL/2-1*fs/ynL;Print待处理信号滤波前幅频特性;ListPlotTransposeF2,YA,FillingAxis,AxesLabelf/Hz,yn:A(f),PlotRangeAll第9页共11页Print待处理信号滤波前相频特性;ListPlotTransposeF2,ArgYn1;ynL2,FillingAxis,AxesLabelf/Hz,yn:(f),PlotRangeAllL=bL+ynL-1;t1=Date;z=ConstantArray0,L;s=JoinConstantArray0,bL-1,s2ynL;1;-1;Forj=1,jL,j+,s=RotateRights,1;zj=b.s1;bL;result=zn;L-n(*n滤波器单位冲激响应采样长度bL的一半*);Print待处理信号滤波后波形;ListPlotresult,JoinedTrue,InterpolationOrder2,AxesLabelt,滤波后:y(t)Ya=2Fourierresult,FourierParameters-1,-1;Ya1/=2;YaA=AbsYa1;ynL2;Fa=Range0,ynL2-1*fs/ynL;Print待处理信号滤波后幅频特性;ListPlotTransposeFa,YaA,FillingAxis,AxesLabelf/Hz,滤波后:A(f),PlotRangeAllListPlotTransposeFa,ArgYa1;ynL2,FillingAxis,AxesLabelf/Hz,滤波后:(f),PlotRangeAlldt=(Date-t1).0,0,0,3600,60,1;Printdt=,dt;输出结果：滤波器单位冲激响应时域波形如图6.1.1：0.51.01.52.02.5t0.10.00.10.20.30.40.5ht图6.1.1滤波器单位冲激响应时域波形滤波器幅频特性如图6.1.2：第10页共11页024681012fHz0.20.40.60.81.0Hf图6.1.2滤波器幅频特性波形滤波器相频特性如图6.1.3：24681012fHz3211f图6.1.3滤波器相频特性波形待处理信号滤波器波形如图6.1.4：0.51.01.52.02.53.0t0.50.00.51.0yt图6.1.4待处理信号滤波器波形待处理信号滤波前幅频特性如图6.1.5：第11页共11页图6.1.5待处理信号滤波前幅频特性待处理信号滤波前相频特性如图6.1.6：图6.1.6待处理信号滤波前相频特性待处理信号滤波后波形如图6.1.7：20406080t0.20.00.20.40.6滤波后:yt图6.1.7待处理信号滤波后波形第12页共11页待处理信号滤波后幅频特性如图6.1.8：图6.1.8待处理信号滤波后幅频特性波形待处理信号滤波后相频特性如图5.1.9：图6.1.9待处理信号滤波后相频特性波形第13页共11页7结语毕业设计是培养我们综合运用所学知识，并发现，提出，分析和解决实际问题，锻炼实践能力的重要环节，也是对学生实际工作能力的具体训练和考察过程。通过两个多月的全心投入和辛苦努力，终于完成了因果序列卷积的矢量乘法算法程序设计，设计的过程就是综合运用所学知识和学习新知识的过程。由于经验不足，刚开始遇到了不少的困难，但是最后在指导老师的帮助下都一一解决。本课题对因果序列卷积的矢量乘法算法程序进行了设计，还进行了物理实例的讨论，更贴近生活，也对其他几种算法进行全面、详细地阐述与对比。本次设计已经达到了设计题目的要求，但在程序上还留有进一步改善的余地。当然由于我在理论和实践方面存在一定的不足，所以在算法设计和实现功能上难免有不足和没有想到的地方，还请各位老师给予指正。第14页共11页致谢时光荏苒，大学四年的生活也转瞬即逝，在此，向曾经关心、支持、鼓励我的老师、同学致以最真诚的谢意。感谢指导我论文的龙老师，在整个过程中，从选题、资料搜集到论文撰写、修改定稿对我的悉心指导，老师传授给我的不仅是专业知识及科研技能，更重要的是老师深厚的学术造诣、严谨的治学态度和对事业的忘我精神深深地影响着我，使我终身受益。同时也要感谢其他帮助我的老师和同学,他们是我人生中最大的一笔财富。毕业在即，我将努力秉承“明德、励志、博学、笃行”的校训,认真学习，努力工作，以自己的实际行动努力践行科学发展观，为社会做贡献，为母校添光彩！第15页共11页参考文献1吴大正.信号与线性系统M.北京:北京高等教育出版社.2005.84-110.2李昌利,有限长序列卷积和求解法J.电气电子教学学报.2008,3D(01)45-473李培芳.离散信号与系统分析M.浙江:浙江大学出版社.2001.1.14高西全,丁玉美.数字信号处理（第三版M.西安:西安电子科技大学出版社.2008.45吴晓艳,李莹等.C#语言程序设计M.北京:清华大学出版社.2011.56冯博琴,信号、系统与信号处理M.北京:机械工业出版社.20017姜立东.VHDL语言程序设计及应用M.北京:北京邮电大学出版社.20048孔凡年.时域卷积反演的新方法N.电子学报,1985.79吴大正,杨林耀,张永瑞.信号与线性系统分析J.北京邮电大学出版社.199810卢稳超,石成业.控制系统中的信号J.黑龙江科技信息.2009第16页共11页附录：1利用矩阵乘法计算卷积：tk=0,0,0,3600,60,1;x=TableCos*k/171Exp-0.0002k,k,0,999;y=Table(0.999)k,k,0,999;t1=Date;n=Lengthx;m=Lengthy;L=n+m-1;y0=PadRighty,L;tem=y0;m=tem;Forj=1,jn,j+,tem=RotateRighttem,1;AppendTom,temz2=x.m;Clearm;ListPlotz2,AspectRatio1/4,PlotRangeAll,PlotStylePointSize0.0001Print按矩阵乘法计算耗时:,(Date-t1).tk,sec;(*x=x1,x2,x3,.;y=y1,y2,y3,.;2ListConvolve算法x=2,3,4;y=1,2,9,0,3,5;t1=Date;z=ListConvolvex,y,1,-1,0;dt=(Date-t1).0,0,0,3600,60,1;Printdt=,dt;Printz=,z;3矢量乘法Z1=ListConvolvex1,x2,x3,u1,u2,u3,u4,u5,1,-1,0x=x1,x2,x3;xL=Lengthx;y=u1,u2,u3,u4,u5;yL=Lengthy;L=xL+yL-1;t1=Date;z=ConstantArray0,Ls=JoinConstantArray0,xL-1,yyL;1;-1Forj=1,jL,j+,s=RotateRights,1;zj=x.s1;xL;z(*Printdt=,(Date-t1).0,0,0,3600,60,1;t2=Date;ListConvolvex,y,1,-1,0;第17页共11页Printdt=,(Date-t2).0,0,0,3600,60,1;*)4序列卷积的FFT算法x=RandomInteger1,5,10000;xL=Lengthx;y=RandomInteger1,5,20000;yL=Lengthy;L=xL+yL-1;t1=Date;xn=PadRightx,L;yn=PadRighty,L;X=Fourierxn,FourierParameters1,-1;Y=Fourieryn,FourierParameters1,-1;Z=X*Y;zfft=FourierZ,FourierParameters-1,1;Printdtfft=,(Date-t1).0,0,0,3600,60,1;t2=Date;zlc=ListConvolvex,y,1,-1,0;Printdtlc=,(Date-t2).0,0,0,3600,60,1;er=MaxAbszfft-zlc6设计滤波器(*确定阶数255；确定通带上限频率；fp=fs/m;b=fir1(255,w);savezmxbcd*)b=ImportD:ProgramFilesMATLABR2012abinzxmfilter.mat/Flatten;s=ImportC:UsersAdministratorDesktop111.txt,Table/Flatten;yn=-0.2,0.05,0.09,0.1,0.1,0.1,0.2,0.3,0.3,0.2,0.2,0.2,0.29,0.29,0.29,0.21,0.21,0.3,0.2,1.3,-0.2,0.1,0.12,0.12,0.1,0.18,0.2,0.21,0.12,0.11,0.1,0.09,0,-0.01,-0.1,-0.12,-0.15,-0.1,-0.1,-0.2,0.9,-0.5,-0.2,-0.15,-0.2,-0.1,0,0.15,0,-0.02,-0.1,-0.1,-0.08,-0.05,0,0,0,0.1,0,-0.01,1,-0.3,0,0.02,0,0,0.1,0.2,0.12,0.1,0,-0.1,-0.09,-0.09,0,0,0,0.1,0.1,0,1.1;ynL=Lengthyn;ListPlotyn,JoinedTrue,PlotRangeAll,AspectRatio1/3bL=Lengthb;n=FloorbL/2;Ts=0.04;fs=1/Ts;tsb=Range0,bL-1*Ts;Print滤波器单位冲激响应时域波形;ListPlotTransposetsb,b,JoinedTrue,AxesLabelt,h(t),PlotRangeAll,InterpolationOr