因果序列卷积的 ListConvolve算法程序设计定稿

因果序列卷积的ListConvolve算法程序设计摘要采用调用Mathematica的内部函数ListConvolve进行因果序列卷积算法的程序设计，并编写因果序列卷积的ListConvolve算法和其它两种算法的程序，以及与其它两种算法的比较。并从实例展示其程序用法和高效性。关键字Mathematica程序;卷积;ListConvolve算法;因果序列CausalsequenceconvolutionalgorithmdesignListConvolveAbstractInvoketheMathematicaListConvolveinternalfunctionisadoptedtoimprovetheprogramLcausalsequenceconvolutionalgorithmdesign,andwriteacausalsequenceofconvolutionListConvolvealgorithmprogramandothertwokindsofalgorithm,andcomparedwiththeothertwoalgorithms.Fromexamplesshowitsprogramusageandefficiency.KeywordsMathematicaprogram;Convolution;ListConvolvealgorithm;Causalsequence；目录引言.11因果序列卷积的意义.22常见有限长因果序列卷积的计算方法.22.1ListConvolve算法.22.2矩阵乘法算法.22.3矢量乘法算法.23各种卷积计算方法程序实现的编程思路.23.1各种算法编程前的准备.23.1.1构造滤波器.23.1.2计算滤波器时域系统函数的采样数据.33.1.3构造信号.33.2因果序列卷积的各种算法编程实现.53.2.1ListConvolve算法程序设计.53.2.2矩阵乘法算法程序设计.63.2.3矢量乘法算法程序设计.74各种卷积算法的比较.84.1ListConvolve算法本身的比较.84.2卷积的三种算法比较.85应用实例.96结语.11致谢.12参考文献.12第0页共9页引言Mathematica是由英国科学家斯蒂芬沃尔夫勒姆领导的著名公司沃尔夫勒姆研究公司开发的一款广泛使用的科学计算软件。它拥有强大的数值计算和符号运算能力，是使用最广泛的数学软件之一。Mathematica的发布标志着现代科技计算的开始1。当然Mathematica是世界上通用计算系统中最强大的系统。自从1988发布以来，它已经对如何在科技和其它领域运用计算机产生了深刻的影响。就像人们常说，Mathematica的发布标志着现代科技计算的开始。自从上世纪六十年代以来，在数值、代数、图形、以及其它方面一直有个别的软件包存在。但是，Mathematica的基本概念是用一个连贯的和统一的方法创造一个能适用于科技计算各个方面的软件系统。至于如何实现这一点的关键之处是发明了一种新的计算机符号语言。这种语言可以仅仅用很少量的基本元素制造出广泛的物体，而满足科技计算的广泛性。这在人类历史上还是第一次出现。而且长期以来，Mathematica核心设计的普遍性使得它涉及的领域不断增长，扩大。从刚开始的是一个主要用于数学和科技计算的系统，到而今发展成为许多计算领域的主要力量，Mathematica已经成为世界上最强大的通用计算系统之一。Mathematica应用广泛，Mathematica的用户群中最主要的是科技工作者和其它专业研究人士。不仅如此，Mathematica还被广泛地用于教学中。数学中的许多计算是非常繁琐的，特别是函数的作图费时又费力，而且所画的图形很不规范，所以现在流行用Mathematica符号计算系统进行学习，从高中到研究生院的数以百计的课程都是在使用它，很多问题便可以迎刃而解。此外，随着学生版的出现，Mathematica已经在全世界的学生中流行起来，成为了一个著名的学习，研究工具。在运算领域里，Mathematica作为当前国际控制界最流行的面向工程与科学计算的高级语言，近年来得到了业界的一致认可2-7。而且在控制各种系统的分析、设计算法方面得到了非常广泛的应用，其自身也因此得到了迅速的发展，功能也在不断扩充。便于对通信系统的物理概念和运行过程的直观理解，所以近年来在通信工程专业中得到了广大师生的特别重视和广泛的应用，在理论教学、课程实践环节，以及理论和技术前沿的研究中发挥了非常重要作用。而传统的因果序列卷积运算方法主要可以归纳为图解法,解析法,位移法和乘式法，矩阵乘法和矢量乘法等等。至于图解法直观明了但烦琐,不便于编程计算;而解析法适合于短序列编程,当序列较长时,表达式冗长,不便于编程计算;乘式法与解析法有相同特点,不适合较长的两个序列的卷积和运算;移位法有效的解决了较长序列的卷积和运算,是一种较好的计算方法,特别对于长序列更显示出其计算能力,但是需要考虑数据存储结构问题,而且也不便于理解。而本次毕业设计中主要提出了另一种计算方法,即ListConvolve算法，易于编程，算法高速8。ListConvolve算法适用于因果序列卷积的运算,而且算法高速，占据内存小，因此是一种有效的计算方法.在工程实际中具有较大的应用价值。在此调用Mathematica软件中的软件ListConvolve算法进行程序设计，使其能够实现因果序列卷积的ListConvolve算法的运算，使卷积运算更加简便。许多文献讨论了卷积的计算方法以及对各种算法的总结，离散信号的卷积的多种求解方法，包括解析法、列表法、利用卷积和性质求解，应用单位样值信号求解，ListConvolve算法，矩阵算法及矢量空间算法。一般利用定义式、图解法和卷积性质这些通用方法进行求解连续信号的卷积或离散信号的卷积和。其中，利用定义式及卷积性质的求解，需要根据信号形式做繁琐的积分或求和；利用图解法进行的反转、平移、相乘、求和四个步骤，绘图过程麻烦9。而利用ListConvolve算法则简单高速，但算法不够透明，而矩阵算法则属于算法过程透明，但却占内存过大，矢量空间算法则运算时间过长。总之各种算法都有利有弊。本文将利用矩阵乘法和矢量乘法与ListConvolve算法相比较来突显出ListConvolve算法的优越性。对于本文来说，主要研究了ListConvolve算法的如何调用，矩阵乘法算法和矢量乘法算法的计算方法是如何实现的。以及如何各种算法的程序编写。因为卷积在信号处理过程中是很常见的，所以研究如何实现卷积滤波是非常只有用的。相信卷积滤波会得到广泛应用。第1页共12页1因果序列卷积的意义任意两个序列的卷积并无意义的，只有一个代表离散信号的序列f(k)与描述离散LTI系统的序列g(k)（g(k)是离散LTI系统在单位脉冲序列k在它激励下的零状态响应。）卷积y(k)=f(k)*g(k)才有真实的物理意义，它就是离散LTI系统在输入离散信号f(k)激励下的零状态响应，即系统的输出，能够代表真是的离散信号的序列，只能是因果序列，所以，本文只讨论一个描述输入信号的有限长序列与一个描述系统滤波特性的有限长序列之间的卷积，其结果是输入信号经过滤波器滤波后的结果。因为在信号处理过程中会常遇到卷积计算问题，所以对于研究因果序列卷积的算法是非常有用的。第2页共12页2常见有限长因果序列卷积的计算方法2.1ListConvolve算法ListConvolve算法，对于学者而言，只知道经过编程实践发现，它的算法很优越，省时省空间，目前来说它的时间复杂度和空间复杂度都是最低的。它是由世界著名软件公司（Mathematica软件）的精英们设计的，所以也只有他们清楚其内部算法。这种算法的透明度不高，至于为什么还要对其进行研究，是因为要学习如何熟练的使用它，如何熟练地调用其内部函数，以充分的利用其高速性与占据内存小的特点。2.2矩阵乘法算法矩阵乘法是一种高效的算法可以把一些一维递推优化到log（n)，还可以求路径方案等，所以更是一种应用性极强的算法。矩阵，是线性代数中的基本概念之一。一个mn的矩阵就是mn个数排成m行n列的一个数阵。由于它把许多数据紧凑的集中到了一起，所以有时候可以简便地表示一些复杂的模型。例如：序列x=2,3,4和序列y=4,5,6,7的卷积，可以等效为矩阵234和矩阵作卷积结果为16,22,43,52,45,284567002.3矢量乘法算法将两个长度分别为m与n的序列x，y进行卷积求和，则其卷积结果序列长度为L=m+n-1,其计算过程为，第一个序列保持不变，第二个序列补m-1个零，放在前面，后面加上倒置的第二个序列，构成新的第二个序列，在相乘的过程中，新生成的第二个序列每乘一次就右循环滚动一次（第一次乘就要右循环滚动一次），其求和方式为每一次对应位置相乘求和，即为矢量乘法算法。例如：序列x=2,3,4和序列y=4,5,6,7的卷积，y序列改写为y1=4,0,0,7,6,5y6=0,0,7,6,5,4卷积结果如下：8，22，43，52，45，28第3页共12页3各种卷积计算方法程序实现的编程思路因为有物理意义的卷积是，一个序列代表输入信号，另一个描述滤波器的滤波特性，它们的卷积结果就是信号加到滤波器上以后，滤波的输出信号，所以，要讨论序列卷积，都必须先准备描述滤波器的序列。因此，有如下编程。3.1各种算法编程前的准备3.1.1构造滤波器滤波器设计主要是借助MATLAB，设计程序如下：(*确定阶数255；确定通带上限频率；n=255;w=0.25;b=fir1(n,w);savezjfilter,b;3.1.2计算滤波器时域系统函数的采样数据Ts=0.01;fs=1/Ts;L=Lengthb;Lh=FloorL/2;ts=Range0,L-1*Ts;ListPlotTransposets/Ts,b,Filling-Axis,AxesLabel-t/Ts,h(n),PlotRange-All,PlotStyle-PointSize0.005B=2Fourierb,FourierParameters-1,-1;B1/=2;BA=AbsB1;Lh;CurlyPhi=ArgB1;Lh;F=Range0,Lh-1*(fs/L);ListPlotTransposeF,BA,Filling-Axis,PlotStyle-PointSize0.005,PlotRange-All,AxesLabel-f/Hz,A(f)ListPlotTransposeF,CurlyPhi,Filling-Axis,PlotStyle-PointSize0.005,PlotRange-All,AxesLabel-f/Hz,CurlyPhi(f)3.1.3构造信号A=10,5,4,1,0.5,7,2,4;f=2,6,8,10,13,15,18,20;CurlyPhi0=0,0,-Pi/2,Pi/4,-Pi/8,Pi/12,0,Pi/20;LA=LengthA;xt_=SumAkCos2Pifkt+CurlyPhi0k,k,1,LA;(*p1=Plotxt,t,0,2.56*)sL=2000;m=100;ts2=Range0,sL-1*Ts;xn=xts2;p2=ListPlotTransposets2,xn1;m,(*FillingRuleAxis*)Joined-True,PlotStyle-PointSize0.005,AxesLabel-t/Ts,x(n)程序运行结果如下：第4页共12页滤波器单位冲激响应时域波形如图3.1：图3.1滤波器单位冲激响应时域波形图滤波器幅频特性如图3.2：第5页共12页图3.2滤波器幅频特性波形图滤波器相频特性如图3.3：图3.3滤波器相频特性波形图第6页共12页待处理信号滤波器波形如图3.4：0.20.40.60.81.0tTs2010102030xn图3.4待处理信号滤波器波形3.2因果序列卷积的各种算法编程实现3.2.1ListConvolve算法程序设计Clearx,y,a,b,c,d,u,v,w;tk=0,0,0,3600,60,1;x=TableRandom,k,5000;y=TableRandom,k,8000;t3=Date;z3=ListConvolvex,y,1,-1,0;t4=Date;ListPlotz3,Filling-Axis,AspectRatio-1/4,PlotRange-All,PlotStyle-PointSize0.0001,AxesOrigin-0,0Print引用ListConvolve函数方法计算耗时:,(t4-t3).tk,sec;运行结果：引用ListConvolve函数方法计算耗时:0.*10-8sec运行结果：图3.5ListConvolve算法运行图第7页共12页In3：=Clearu,v,w,x,y,z,p;ListConvolveu,v,w,x,y,z,p,1,-1,0Out4=ux,vx+uy,wx+vy+uz,pu+wy+vz,pv+wz,pw3.2.2矩阵乘法算法程序设计Clearx,y,z,s;tk=0,0,0,3600,60,1;x=TableCosPi*k/171Exp-0.0002k,k,0,79;y=Table(0.999)k,k,0,59;t1=Date;z2=ListConvolvex,y,1,-1,0;t3=Date;n=Lengthx;m=Lengthy;L=n+m-1;y0=PadRighty,L;tem=y0;m=tem;Forj=1,j1/4,PlotRange-All,PlotStyle-PointSize0.0001Print按矩阵乘法计耗时:,(Date-t1).tk,secnListConvolve用时,(t3-t1).tk,sec;204060801001201401020304050图3.6矩阵乘法算法运行图按矩阵乘法计算耗时:0.1080062secListConvolve用时0.0010000sec3.2.3矢量乘法算法程序设计Clearx*,y*,z*,s;tk=0,0,0,3600,60,1;x=TableCosPi*k/171Exp-0.0002k,k,0,9999;第8页共12页y=Table(0.999)k,k,0,9;t1=Date;m=Lengthx;n=Lengthy;L=m+n-1;z=ConstantArray0,L;s=JoinConstantArray0,m-1,yn;1;-1;Forj=1,jAll,PlotStyle-PointSize0.006,Red;p2=ListPlotz21;15,PlotRange-All,PlotStyle-PointSize0.003,Blue;Showp1,p2Print矢量乘法用时,(t2-t1).tk,secnListConvolve用时,(t3-t2).tk,sec;2468101214246810图3.7矢量乘法算法运行结果矢量乘法用时12.3147043secListConvolve用时0.0840048sec第9页共12页4各种卷积算法的比较4.1ListConvolve算法本身的比较表4.1ListConvolve算法不同长序列耗时比较算法名称序列（长度）耗时1250008000极短，可以忽略不计ListConvolve算法10000200000.0156000sec50000800000.0468001secListConvolve算法本身是一种非常高速的算法，在序列数值不大的情况下算法耗时几乎看不到，但在其对序列长度差异很大的序列进行卷积运算时还是会有耗时上的差异，即长序列耗时大于短序列耗时。4.2卷积的三种算法比较表4.2卷积的三种算法比较算法名称序列（长度）耗时时间，空间复杂度12ListConvolve算法极短都很小矩阵乘法算法50008000无法计算空间复杂度过大矢量乘法算法17.1444299sec都一般根据本次毕业设计对以上三种算法的对比，可以明显体现出ListConvolve算法的高速性，当然，这三种算法也各有利弊，对于ListConvolve算法当然就用来计算长序列数值大的这类序列卷积，而矩阵算法则用于学者做研究，至于矢量空间算法一般不太用。所以，在选择运用那种算法时，要结合它们的特点来选择。第10页共12页5应用实例n=63;wp=0.3;fp=wp*fs/2;b=fir1(n,wp);savezj009b;cdh=-0.000798867,-0.000390074,0.000436796,0.00110852,0.000950907,-0.000255716,-0.0017797,-0.00217304,-0.000463884,0.00253415,0.00424909,0.00232847,-0.00275781,-0.00706324,-0.005921,0.00152669,0.0100959,0.0116871,0.00237057,-0.0123678,-0.0199452,-0.0106029,0.0123052,0.0312472,0.026339,-0.00694942,-0.0481115,-0.0601859,-0.0138395,0.088797,0.208625,0.289004,0.289004,0.208625,0.088797,-0.0138395,-0.0601859,-0.0481115,-0.00694942,0.026339,0.0312472,0.0123052,-0.0106029,-0.0199452,-0.0123678,0.00237057,0.0116871,0.0100959,0.00152669,-0.005921,-0.00706324,-0.00275781,0.00232847,0.00424909,0.00253415,-0.000463884,-0.00217304,-0.0017797,-0.000255716,0.000950907,0.00110852,0.000436796,-0.000390074,-0.000798867;fs=100;Ts=1/fs;ta=0;tb=5;n=(tb-ta)/Ts;ts=Range0,n-1*Ts;L=Lengthh;ts1=Range0,L-1*Ts;ListPlotTransposets1,h,Filling-Axis,PlotRange-All,PlotStyle-PointSize0.008,AxesLabel-Subscriptt,k/s,h(Subscriptt,k)F=Range0,L/2-1*fs/L;H=Fourierh,FourierParameters-1,-1;ListPlotTransposeF,AbsH1;L/2,Joined-True,PlotStyle-PointSize0.008,AxesLabel-f/Hz,|H(f)|st_=10Cos10Pit+5Cos40Pit+2Cos80Pit+0.5Cos100Pit;sn=sts;ListPlotTransposets,sn1;100,Filling-Axis,PlotRange-All,PlotStyle-PointSize0.008,AxesLabel-Subscriptt,k/s,s(Subscriptt,k)zn=ListConvolveh,sn,1,-1,0;yn=znL/2;n+L/2-1;ListPlotTransposets,yn1;100,Filling-Axis,PlotRange-All,PlotStyle-PointSize0.008,AxesLabel-Subscriptt,k/s,y(Subscriptt,k)程序运行结果：第11页共12页图5.1滤波器的单位脉冲序列杆状图010203040fHz0.20.40.60.81.0Hf图5.2幅频特性曲线折线图第12页共12页图5.3输入信号杆状图图5.4输出信号杆状图第13