毕业论文-因果序列卷积的矢量乘法算法程序设计.doc

陕西理工学院毕业设计 毕业设计 题 目 因果序列卷积的矢量乘法算法程序设计 学生姓 2014 年 5 月 14 日 陕西理工学院毕业设计 因果序列卷积的矢量乘法算法程序设计 摘要 卷积是信号后期处理中的重要概念。离散信号的卷积计算在文献中有多种解法。我们提出了用矢量乘 法计算离散信号卷积和的方法。与传统计算方法相比, 该方法避开了繁琐的做图求和过程, 使计算更简便、快速。 同时，也对几种方法进行了对比和讨论。 关键字 因果序列;卷积和;矢量乘法程序设计 The vector multiplication of the causal sequence of convolution algorithm programming Zhou Mixue (Grade10,Class1,Major Electronic Information Science and Technology Institute of Physics and Electronics,Shannxi University of Technology,Hanzhong,723001) Tutor:Long Shuming Abstract: Convolution is an important concept in “Signals post-processing”.Methods for calculating causal sequence of convolution are given.The method to compute the convolution sum of discrete signals byusing Mathematica is introduced. In comparison with the traditional computation, this method avoids the complicated procedures such as drawing.Finally,comparisons and discussions of the these methods are provided。 Keywords: causal sequence; convolution sum; vector multiplication;program 陕西理工学院毕业设计 第 1 页 共 11 页 目录 1 引言 1 1.1 离散序列卷积的意义 1 1.2 卷积的用途 2 2 课题研究方法 2 2.1 卷积的计算机实现 2 2.2 课题实现工具 2 3 序列卷积常用算法及其实现思路 2 3.1 定义法 2 3.2 矩阵乘法 2 3.3 LISTCONVOLVE算法3 4 卷积矢量乘法算法的程序设计 3 5 各种卷积算法优劣比较 4 5.1 矩阵乘法 4 5.2 LISTCONVOLVE算法4 5.3 矢量乘法算法 4 5.4 FFT 算法5 5.5 列表比较 5 6 卷积算法程序应用实例 6 6.1 设计滤波器 6 6.2 物理实例 6 7 结语 10 致谢 10 参考文献 11 陕西理工学院毕业设计 第 2 页 共 11 页 1 引言 在很多人眼中，卷积这个概念是比较神秘且不容易理解的。卷积在数学、物理学、电子工程、 信号处理、计算机科学中极为重要。离散信号卷积在电子通信领域的应用非常广泛，也是工程应用 的基础。所以快速有效的计算离散序列的卷积，是人们一直很关心的一个问题。 1.1 离散序列卷积的意义 在科学技术飞速发展的今天，现代通信正在蓬勃发展，如何快速准确、方便易行地进行远距离 信息传递，这是人类追求的目标。要想实现这一目标就必须掌握离散卷积的运算。卷积是在信号与 线性系统的基础上和背景中出的，卷积关系很重要的一种情况就是，在信号与线性系统中或数字信 号处理中的卷积定理。因为利用该定理，可以将进行时间域或空间域中的卷积运算等价为频率域的 相乘运算，从而利用 FFT 等快速算法，实现更有效的计算，节省运算时间，提升效率。 信号通过信号放大系统就是卷积，因为系统是对信号做了线性加工。进一步推广照相系统是卷 积，视觉系统也是卷积。我们看到的物体是物体本身与我们的视觉系统的卷积，因为同一个物体， 视力不同的人看到的形象不尽相同。同样我们听到的声音是声源发出的声音与听觉系统的卷积。世 界上所有的系统都是卷积，例如脉搏是心跳的卷积，水压的波动时泵水电机运转动力的卷积，只要 系统是线性的或近似线性的。所以当我们正确认识了卷积的意义之后就会发现卷积是无处不在。但 是卷积的意义建立在这两个前提之下：一是任意的数字信号都可以表示成单位脉冲之线性组合，二 是该系统也是线性的。实际上，我们是根据我们需要待处理信号的形式，来设计所谓的系统传递函 数，那么这个系统的传递函数和输入信号，在数学上的形式就是所谓的卷积关系。 1.2 卷积的用途 卷积的用途就是通过物理系统对信号进行滤波运算，然后对其信号进行幅度的放大及延时。卷 积在信号中主要是时域、频域的相互转换。卷积在其它领域也有很多运用,例如在编码中，有卷积 码，就是运用原码中前后序列的码字确定当前的编码输出。卷积在分析系统的零状态响应过程中也 有着重要作用, 2 课题研究方法 2.1 卷积的计算机实现 在计算机上进行卷积运算就是做一个窗口叠加运算，我们一般看的计算机信号处理方面的教材， 例如数字图像处理 ，其中用卷积方的法就实现了图像的积分、微分、锐化、平滑、去噪声等各 种运算。 2.2 课题实现工具 随着计算机技术行业的发展, 越来越多的计算问题都已经交给计算机来处理。Matlab 和 Mathematica 作为优秀的科学计算的软件, 在信号处理与通信、工程计算、图像处理技术等领域均 得到了广泛的应用。MATLAB 主要用于数据可视化、算法开发、数值计算及数据分析的交互式环 境和高级计算语言。Mathematics 是一个教育软件，设计给 Microsoft Windows，使用户能够解决数 学的科学问题。由微软开发并维护，它主要作为学生和老师的学习工具。Mathematica 拥有强大的 符号及数值运算能力以及方便实用的绘图功能, 所以应用 Mathematica 总会让人身心舒畅。用 Mathematica 来计算积分、求和、画图都非常容易, 所以可以利用这些特点, 做好计算卷积的程序 包或函数, 使计算过程大大简化, 并可以得到精确的数值解或者解析解。 3 序列卷积常用算法及其实现思路 3.1 定义法 离散信号的卷积和定义为1： (式 3.1)( )( )* ( )( ) () i y kf kg kf i g ki 陕西理工学院毕业设计 第 3 页 共 11 页 式中：f(k)和 g(k)为离散时间序列；i 为虚设求和变量，卷积的计算结果为另一新的序列。 用如下例子说明： 已知, ,求卷积。( )( )(01) n x na u na( )( )h nu n( )( )* ( )y nx nh n 解答： 1 0 1 ( )( )* ( )( ) ()( )( ) 1 mm n n mm a y nx nh na u m u nma u nu n a 当时n 1 ( ) 1 y n a 3.2 矩阵乘法 用如下例子说明： 求序列 x=2,3,4和序列 y=1,2,9,0,3,5的卷积，可以等效为矩阵和矩阵432 相乘。 53092100 05309210 00530921 运用 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,5 Print z=,z; 输出结果： z=2,7,28,35,42,19,27,20 3.3 ListConvolve 算法 ListConvolve 算法是 mathematics 内部一个算法，主要是进行循环卷积，基本用法:ListConvolv eker, list 构造核 ker 和 list 的卷积。出现的新概念是核。核即指核心，便是值输入的信号。list 自然是单位冲击信号对应的响应了。下面是在 mathematics 上进行的一个简单编程： x=2,3,4; y=1,2,9,0,3,5; z=ListConvolvex,y,1,-1,0; Print z=, z; 输出结果： z= 2,7,28,35,42,19,27,20 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+, 陕西理工学院毕业设计 第 4 页 共 11 页 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,p2 Print矢量乘法用时, (t2-t1).tk , secnListConvolvey 用时, (t3-t2).tk,sec; 5 各种卷积算法优劣比较 前文已主要介绍了几种求解卷积的算法，但是各种算法从运算效率和占用内存上来说肯定存在 一定的差异，在工作和学习中为了选择更合适的，很有必要对各种算法进行讨论。在选择算法时我 们必须本着简洁，实用，便捷的原则。 5.1 矩阵乘法 tk=0,0,0,3600,60,1; x=Table Cos*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;*) 陕西理工学院毕业设计 第 5 页 共 11 页 t3=Date; ListConvolvex,y,1,-1,0; dt2=(Date-t3).0,0,0,3600,60,1; Printdt2=,dt2; 运算结果： dt2=_0.0010001 5.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; 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,p2 Print矢量乘法用时,(t2-t1).tk, secnListConvolve 用时,(t3-t2).tk,sec ; 运算结果： 矢量乘法用时 10.6316081sec ListConvolve 用时 0.0040002sec 5.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; 陕西理工学院毕业设计 第 6 页 共 11 页 Printdtlc=,(Date-t2).0,0,0,3600,60,1; er=MaxAbszfft-zlc 5.5 列表比较 表 5.1 算法复杂度比较 序列一长度 序列二长度 计算耗时 占用存储空间 算法明晰度 矩阵乘法算法 999 999 0.1 159976 明晰 ListConvolve 算法 9999 9999 0.005 28747 明晰 矢量乘法算法 9999 9999 0.004 159976 明晰 FFT 算法 9999 9999 0.01 159976 明晰 结论：矩阵乘法不适合用于长序列，否则会造成 mathematics 内核自动停止运行，这时就需要 关闭其它程序来再一次尝试。在运行短序列时它的效率还是比较高的。正因为矩阵乘法只适合短序 列所以教学过程中一般使用这个方法。ListConvolve 算法不仅适用于长序列也适用于短序列，而且 效率最高，非常节省时间。矢量乘法算法不仅适用于短序列，也适用于长序列，但是长序列时相当 费时。FFT 算法长序列和短序列都适用，较长序列一般考虑 ListConvolve 算法和 FFT 算法。相较之 下 ListConvolve算法最优，接下来是FFT算法，矢量乘法算法，最后是矩阵乘法。 6 卷积算法程序应用实例 6.1 设计滤波器 滤波器设计主要是借助 MATLAB，设计程序如下： (*确定阶数 255；确定通带上限频率； fp=fs/m; b=fir1(255,w); save zmx b cd *) 6.2 物理实例 为了更好的了解矢量乘法算法，我找了一份正常人的心电图作为实例来说明。首先从 MATLAB 调用滤波器,程序如下： b=ImportD:Program FilesMATLABR2012abinzxmfilter.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/3 bL=Lengthb; n=FloorbL/2; Ts=0.04;fs=1/Ts; tsb=Range0,bL-1*Ts; Print滤波器单位冲激响应时域波形; ListPlotTransposetsb,b,JoinedTrue,AxesLabelt,h(t),PlotRangeAll, InterpolationOrder2 H=Fourierb,FourierParameters1,-1; 陕西理工学院毕业设计 第 7 页 共 11 页 HA=AbsH1;n; F=Range0,bL/2-1*fs/bL; Print滤波器幅频特性; ListPlotTransposeF,HA,JoinedTrue,AxesLabelf/Hz,H(f),PlotRangeAll ListPlotTransposeF,ArgH1;n,JoinedTrue,AxesLabelf/Hz,(f),PlotRangeAll ynL2=FloorynL/2; ts=Range0,ynL-1*0.04; Print待处理信号滤波前波形; ListPlotTransposets,yn,JoinedTrue,AxesLabelt,y(t),PlotRangeAll,InterpolationO rder2 Yn=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 Print待处理信号滤波前相频特性; ListPlotTransposeF2,ArgYn1;ynL2,FillingAxis,AxesLabelf/Hz,yn: (f),PlotRangeAll L=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), PlotRangeAll ListPlotTransposeFa,ArgYa1;ynL2,FillingAxis,AxesLabelf/Hz,滤波后: (f),PlotRangeAll dt=(Date-t1).0,0,0,3600,60,1; Printdt=,dt; 输出结果： 滤波器单位冲激响应时域波形如图 6.1.1： 陕西理工学院毕业设计 第 8 页 共 11 页 0.5 1.0 1.5 2.0 2.5 t 0.1 0.0 0.1 0.2 0.3 0.4 0.5 ht 图 6.1.1 滤波器单位冲激响应时域波形 滤波器幅频特性如图 6.1.2： 02468 10 12 fHz 0.2 0.4 0.6 0.8 1.0 Hf 图 6.1.2 滤波器幅频特性波形 滤波器相频特性如图 6.1.3： 2468 10 12 fHz 3 2 1 1 f 图 6.1.3 滤波器相频特性波形 待处理信号滤波器波形如图 6.1.4： 陕西理工学院毕业设计 第 9 页 共 11 页 0.51.0 1.52.0 2.53.0 t 0.5 0.0 0.5 1.0 yt 图 6.1.4 待处理信号滤波器波形 待处理信号滤波前幅频特性如图 6.1.5： 图 6.1.5 待处理信号滤波前幅频特性 待处理信号滤波前相频特性如图 6.1.6： 图 6.1.6 待处理信号滤波前相频特性 待处理信号滤波后波形如图 6.1.7： 陕西理工学院毕业设计 第 10 页 共 11 页 20406080 t 0.2 0.0 0.2 0.4 0.6 滤波后:yt 图 6.1.7 待处理信号滤波后波形 待处理信号滤波后幅频特性如图 6.1.8： 图 6.1.8 待处理信号滤波后幅频特性波形 待处理信号滤波后相频特性如图 5.1.9： 图 6.1.9 待处理信号滤波后相频特性波形 陕西理工学院毕业设计 第 11 页 共 11 页 7 结语 毕业设计是培养我们综合运用所学知识，并发现，提