有限长序列线性卷积矢量乘法算法程序设计.docx
《有限长序列线性卷积矢量乘法算法程序设计.docx》由会员分享，可在线阅读，更多相关《有限长序列线性卷积矢量乘法算法程序设计.docx（22页珍藏版）》请在冰点文库上搜索。
有限长序列线性卷积矢量乘法算法程序设计
毕业设计（论文）
题目有限长序列线性卷积矢量乘法算法程序设计
学生姓名杨长兴学号1110064080
所在院（系）物理与电信工程学院
专业班级电子信息科学与技术1103班
指导教师龙姝明
完成地点博远楼A1109,C1207
有限长序列线性卷积矢量乘法算法程序设计
杨长兴
（陕西理工学院物电学院电子信息科学与技术1103班，陕西汉中723000）
指导老师：
龙姝明
[摘要]有限长序列线性卷积现有矩阵乘法、FFT方法、ListConvolve函数方法、卷积定义法等计算方法。
ListConvolve函数方法高效快捷，但我们无从知道它是如何完成计算的。
FFT方法对于两个长序列的卷积计算效率很高，但对短序列来说，反倒繁琐。
矩阵乘法算法对两个短序列的卷积计算简单明了，但长序列的卷积计算，在一般计算机上因内存的限制根本无法计算。
为此我们提出了卷积的矢量乘法算法，算法简洁，程序实现方便，计算速度很快，与其它方法相比，仅次于ListConvolve函数法和FFT方法。
[关键字]有限长序列;线性卷积;矢量乘法
Vectormultiplicationalgorithmprogrammingoffinitelengthsequencelinearconvolution
YangChangxing
（Grade11,Class3,MajorElectronicInformationScienceandTechnology
DepartmentofPhysics,ShannxiUniversityofTechnology,Hanzhong,723001,shannxi）Tutor:
LongShuming
AbstractFinitelengthsequencelinearconvolutionmatrixmultiplication,FFTmethod,existingListConvolvefunctionmethod,theconvolutioncalculationmethods,suchasdefinitionmethod.ListConvolvefunctionmethodishigh,butwedon'tknowifitishowtoperformthecalculation.FFTmethodfortwolongsequenceofconvolutioncomputationefficiencyisveryhigh,butforshortsequences,butcumbersome.Matrixmultiplicationalgorithmonthetwoshortsequenceofconvolutioncalculationsimpleandclear,butalongsequenceofconvolutioncomputation,becauseofthememorylimitongeneralcomputercannotcalculate.Forthis,weputforwardtheconvolutionvectormultiplicationalgorithm,algorithmisconcise,programimplementationisconvenient,fastcalculation,comparedwithothermethods,secondonlytoListConvolvefunctionmethodandtheFFTmethod.
Keywordsfinitelengthsequence,linearconvolution,vectormultiplication
目录
1序列卷积的意义1
1.1离散序列卷积的意义1
1.2卷积的用途1
1.3卷积的程序实现1
1.4课题研究工具1
2序列卷积常用算法及程序实现1
2.1定义法2
2.2矩阵乘法2
2.3ListConvolve算法3
2.4FFT算法4
3卷积矢量乘法算法的程序设计4
4各种时域卷积算法的优缺点以及适用对象5
5卷积算法程序应用实例6
5.1设计滤波器6
5.2滤波实例6
结语8
附录9
在很多人眼中，卷积这个概念是比较神秘且不容易理解的。
卷积在数学、物理学、电子工程、信号处理、计算机科学中极为重要。
离散信号卷积在电子通信领域的应用非常广泛，也是工程应用的基础。
所以快速有效的计算离散序列的卷积，是人们一直很关心的一个问题。
1序列卷积的意义
1.1离散序列卷积的意义
卷积是数学中的概念，对应在物理上的事实就是把输入信号加在系统上以后，系统对输入信号进行加工处理后有一个输出信号，输出信号和输入信号之间可以用卷积来连接它们。
也就是说系统对输入信号的加工过程等价于系统与输入信号做卷积。
卷积是系统功能的一种数学描述。
在物理实验中我们放大信号时，需要给定一个输入信号f（k），通过信号放大系统h（k），得到放大后的输出信号y（k）；他们之间的关系就是y（k）=f（k）*h（k）;在这里，我们就可以看出来卷积是信号放大系统功能的一种数学描述[1]。
序列卷积是求离散系统零状态响应yzs（k）的重要方法（“唯一”途径）。
1.2卷积的用途
卷积就是通过物理系统h（k）对输入信号f（k）进行加工处理（对其信号进行幅度的放大及相位的延时）后得到输出信号y（k），这是物理设备的运转；我们也可以用数学方法，对物理系统h（k）和输入信号f（k）做卷积运算得到y（k）即:
y（k）=h（k）*f（k）；这两种方法达到的是同样的效果。
卷积在分析系统的零状态响应过程中也有着重要作用[2]。
1.3卷积的程序实现
在计算机上进行卷积运算就是做窗口叠加运算，我们一般看的计算机信号处理方面的教材，例如《数字图像处理》，其中用卷积方法就实现了图像的积分、微分、锐化、平滑、去噪声等各种运算[3]。
1.4课题研究工具
随着计算机技术行业的发展,越来越多的计算问题都已经交给计算机来处理。
Mathematica作为优秀的科学计算的软件,在信号处理与通信、工程计算、图像处理技术等领域均得到了广泛的应用。
Matlab主要用于数据可视化、算法开发、数值计算及数据分析的交互式环境和高级语言。
Mathematica是一个教育软件，设计给Microsoft Windows，使用户能够解决数学的科学问题。
由微软开发并维护，它主要作为学生和老师的学习工具。
Mathematica拥有强大的符号及数值运算能力以及方便实用的绘图功能,所以应用Mathematica总会让人身心舒畅。
用Mathematica来计算积分、求和、画图都非常容易,所以可以利用这些特点,做好计算卷积的程序包或函数,使计算过程大大简化,并可以得到精确的数值解或者解析解[8]。
2序列卷积常用算法及程序实现
序列卷积的算法有很多种，有定义法、矩阵乘法、矢量乘法、ListConvolve算法以及FFT算法等。
本文着重讨论的是有限长序列线性卷积的矢量乘法算法，最后要进行各种算法间的优缺点进行比较。
研究之初，先用Mathematica软件构造两个有限长序列x（n）和y（n）。
用Mathematica编程，先定义两个连续信号x（t）与y（t），再对它们进行等间隔采样，产生两个离散信号。
程序实现如下：
Clear[x,y];
T=0.5;
x[t_]=5（2t/T）（UnitStep[t]-UnitStep[t-T/2]）+10（1-t/T）（UnitStep[t-T/2]-UnitStep[t-T]）；
Plot[x[Mod[t,T]],{t,0,2T}]
ts=Range[0,9999]*T/10;
xn=x[Mod[ts,T]];
y[t_]=1+Sin[πt]+2Sin[50πt];
Plot[y[t],{t,0,2T}]
yn=y[ts];
2.1定义法
离散线性卷积的定义为[5-7]：
设序列x（n）的长度为xl、序列y（n）长度为yl，对两序列做线性卷积得到长度为xl+yl-1的序列z（n）
举例如下:
x=（a,b,c）,y=（u,v,w,q）,则：
对于有限长序列x（n）与y（n）的定义法程序实现如下：
t1=Date[];
Clear[x,y];
T=0.5;
x[t_]=5（2t/T）（UnitStep[t]-UnitStep[t-T/2]）+10（1-t/T）（UnitStep[t-T/2]-UnitStep[t-T]）；
Plot[x[Mod[t,T]],{t,0,2T}]
ts=Range[0,9999]*T/10;
xn=x[Mod[ts,T]];
y[t_]=1+Sin[πt]+2Sin[50πt];
Plot[y[t],{t,0,2T}]
yn=y[ts];
z=Table[Sum[xn[[k+1]]If[0<=n-k<=yL-1,yn[[n-k+1]],0],{k,0,xL-1}],{n,0,xL+yL-2}];
t2=Date[];
dt=（t2-t1）.{0,0,0,3600,60,1}
ListPlot[xn,Filling->Axis]
ListPlot[yn,Filling->Axis]
ListPlot[z,Filling->Axis]
2.2矩阵乘法
矩阵乘法的规则如下：
（1）把第一个序列作为一个单行矩阵放在左边，作为第一个矩阵。
（2）把第二个序列中的元素作为第二个矩阵的第一行的前个yl元素，同时在后面补xl-1个零，使得第一行总长度等于xl+yl-1，也就是两序列长度之和减一。
第一行有了之后，从第二行开始始终是第一行依次右循环一个数据位置，第二个矩阵的行数与第一个矩阵的元素个数一致。
（3）两矩阵做矩阵乘法，给出卷积结果序列。
也就是第一个矩阵一行与第二个矩阵的一列进行对应元素相乘求和。
举例如下：
x=（a,b,c）,y=（u,v,w,q），则：
对于有限长序列x（n）与y（n）的矩阵乘法程序实现如下：
t1=Date[];
Clear[x,y];
T=0.5;
x[t_]=5（2t/T）（UnitStep[t]-UnitStep[t-T/2]）+10（1-t/T）（UnitStep[t-T/2]-UnitStep[t-T]）；
Plot[x[Mod[t,T]],{t,0,2T}]
ts=Range[0,9999]*T/10;
xn=x[Mod[ts,T]];
y[t_]=1+Sin[πt]+2Sin[50πt];
Plot[y[t],{t,0,2T}]
yn=y[ts];
y1=PadRight[yn,xL+yL-1];（*给y（n）补xL-1个零*）
ym=Table[0,{k,xL},{j,xL+yL-1}];（*造一个xL行xL+yL-1列的空矩阵*）
ym[[1]]=y1;
yt=y1;
For[j=2,j<=xL,j++,yt=RotateRight[yt];ym[[j]]=yt];（*循环右移生成矩阵*）
z2=xn.ym（*x（n）与y（n）对应的矩阵做矩阵乘法*）;
t2=Date[];
dt=（t2-t1）.{0,0,0,3600,60,1}
ListPlot[z2,Filling->Axis]
2.3ListConvolve算法
ListConvolve算法是Mathematica内部的一个函数，用于计算循环卷积[4]，用其计算线性卷积的调用格式是
ListConvolve[x,y,{1,-1},0]
下面是用Mathematica的ListConvolve函数编程计算线性卷积所需的程序语句组：
x={2,3,4};
y={1,2,9,0,3,5};
z=ListConvolve[x,y,{1,-1},0];
Print["z=",z];
输出结果：
z={2,7,28,35,42,19,27,20}
对于有限长序列x（n）与y（n）的Listconvolve算法程序实现如下：
Clear[x,y];
T=0.5;
x[t_]=5（2t/T）（UnitStep[t]-UnitStep[t-T/2]）+10（1-t/T）（UnitStep[t-T/2]-UnitStep[t-T]）；
Plot[x[Mod[t,T]],{t,0,2T}]
ts=Range[0,9999]*T/10;
xn=x[Mod[ts,T]];
y[t_]=1+Sin[πt]+2Sin[50πt];
Plot[y[t],{t,0,2T}]
yn=y[ts];
t1=Date[];
z5=ListConvolve[xn,yn,{1,-1},0];
t2=Date[];
dt=（t2-t1）.{0,0,0,3600,60,1}
ListPlot[z3,Filling->Axis]
2.4FFT算法
两序列线性卷积的FFT算法分四步：
（1）把两序列扩长，对两序列进行补零，使其长度都为两序列长度之和减一。
（2）对两序列求快速傅里叶变换，得到它们的变换结果数组。
（3）把两个结果数组做乘法。
（4）把两结果相乘后取逆变换。
对于有限长序列x（n）与y（n）的FFT算法程序实现如下：
Clear[x,y];
T=0.5;
x[t_]=5（2t/T）（UnitStep[t]-UnitStep[t-T/2]）+10（1-t/T）（UnitStep[t-T/2]-UnitStep[t-T]）；
Plot[x[Mod[t,T]],{t,0,2T}]
ts=Range[0,9999]*T/10;
xn=x[Mod[ts,T]];
y[t_]=1+Sin[πt]+2Sin[50πt];
Plot[y[t],{t,0,2T}]
yn=y[ts];
Clear[xc,yc,t1,t2,X,Y,z4];
t1=Date[];
xc=PadRight[xn,xL+yL-1];
yc=PadRight[yn,xL+yL-1];
X=Fourier[xc,FourierParameters->{1,-1}];
Y=Fourier[yc,FourierParameters->{1,-1}];
z4=Fourier[X*Y,FourierParameters->{-1,1}]//Chop;（*chop函数去掉结果中的虚部*）
t2=Date[];
dt=（t2-t1）.{0,0,0,3600,60,1}
ListPlot[z4,Filling->Axis]
3卷积矢量乘法算法的程序设计
矢量乘法算法是本课题着重研究的点，它的研究意义在于它易于程序实现，在空间时间等方面优秀于其它卷积算法。
具体比较会在后面给出。
对于给定两序列x={x1…xn}和y={y1…ym}，它的卷积矢量乘法算法步骤如下：
（1）yr={ymym-1…y1}（把第二个序列y先倒序）
（2）y2={0…0，ym…y1}（倒序后给它前面补若干个零，零的个数是第一个序列元素个数减一，即n-1个零）
（3）将y2右移一次成为yt={y1，0…0，ym…y2}，从头取n个数据即{y1，0…0}与序列x做矢量乘法，即对应元素相乘求和，结果是卷积的第一个数据。
（4）重复进行。
对于有限长序列x（n）与y（n）的矢量乘法程序实现如下：
t1=Date[];
Clear[x,y];
T=0.5;
x[t_]=5（2t/T）（UnitStep[t]-UnitStep[t-T/2]）+10（1-t/T）（UnitStep[t-T/2]-UnitStep[t-T]）；
Plot[x[Mod[t,T]],{t,0,2T}]
ts=Range[0,9999]*T/10;
xn=x[Mod[ts,T]];
y[t_]=1+Sin[πt]+2Sin[50πt];
Plot[y[t],{t,0,2T}]
yn=y[ts];
yr=yn[[-1;;1;;-1]];
yt=PadLeft[yr,xL+yL-1];Clear[yr];
z3=Table[0,{j,xL+yL-1}];
For[j=1,j<=xL+yL-1,j++,
yt=RotateRight[yt];
z3[[j]]=xn.yt[[1;;xL]]];
t2=Date[];
dt=（t2-t1）.{0,0,0,3600,60,1}
ListPlot[z3,Filling->Axis]
4各种时域卷积算法的优缺点以及适用对象
前面已经介绍完了序列线性卷积的各种算法以及程序实现，大家肯定会有疑惑介绍了这么多到底都有什么区别呢，不难发现每一个程序实现都有t=Date[]这样的表达式，这是测试时间的语句。
当然从程序中的复杂程度还可以看出空间占用的大小，下面来和大家一起分析分析。
给定两个各有10000个元素的序列，定义法需要一个一个元素的去找对应的元素进行卷积，判断、卷积、循环，每次都是乘法再作加法，这样非常的浪费时间，占用空间8*10^8字节；矩阵乘法计算的话，需要先把第二个序列补零，再把第二个序列造一个大大的空数组，已备保存数据，接下来把第二个序列补了零的放在矩阵的第一行，然后开始右循环移动填充矩阵的其余各行，最后进行矩阵乘法，这样时间倒是不多，就是占用空间太大，占用空间2*8*10^12个字节，对于一台4G电脑的电脑，近似有4*10^9个字节，这样一台电脑还算不过来，这只是只有10000个数据，如果是1000000个的话，那恐怕得成千上万的4G才能，是不可能实现的；对于矢量乘法，先把第二个序列倒序，前面补零，然后开始右循环一次，取和第一个序列一样的元素从头取，对应元素相乘求和，这样相当于还是两个10000数据的序列这样占用空间10^8字节，而且时间占用也不多，下面会给出具体比较。
对于这些分析可以看出，这些时域卷积算法里面更易于实现的是矢量乘法，矩阵乘法不适合用于长序列，否则会造成Mathematica内核自动停止运行，这时就需要关闭其它程序来再一次尝试。
在运行短序列时它的效率还是比较高的[9-10]。
正因为矩阵乘法只适合短序列所以教学过程中一般使用这个方法。
对于Listconvolve算法以及FFT算法，它们是普遍适用而且高效的。
所以这些算法的适用度从大到小排列应该是Listconvolve算法、FFT算法、矢量乘法、矩阵乘法、定义法。
下面对于给定x（n）与y（n）的有限长序列进行列表比较一下具体的时间。
表1各种序列卷积算法的时间比较
定义法
矩阵乘法
矢量乘法
Listconvolve
FFT算法
时间
129.34s
199.07s
11.75s
0.17s
0.76s
从表格可以看出上面的分析基本没有问题。
5卷积算法程序应用实例
在物理实验中，我们可以把滤波系统滤波这一过程看做卷积，滤波前的信号f（k）与滤波系统h（k）进行卷积可得到滤波后的信号y（k）。
5.1设计滤波器
滤波器设计主要是借助Matlab，设计程序如下：
w=fs/n;（*期望保留输出信号的最高频率*）
h=fir1（255,w）;
save'ycx'h
5.2滤波实例
为了更好的了解矢量乘法计算卷积的过程，我们设计一个离散滤波器，把数据存入h内。
我们用Matlab设计的滤波器单位冲击序列为：
h={0,-0.0000769399,-0.000144038,-0.00019133,-0.000211254,-0.000199738,-0.000156932,-0.0000874399,0,0.0000934227,0.000179066,0.0002432,0.000274144,0.000264206,0.000211247,0.000119582,0,-0.000131228,-0.000254311,-0.000348695,-0.000396259,-0.000384493,-0.000309138,-0.000175775,0,0.000194037,0.000376657,0.000516923,0.000587586,0.000569955,0.000457869,0.00026001,0,-0.000285975,-0.000553859,-0.000758204,-0.000859512,-0.000831324,-0.00066583,-0.000376929,0,0.000411904,0.00079513,0.00108489,0.00122577,0.00118165,0.000943302,0.000532265,0,-0.000577935,-0.00111214,-0.00151277,-0.00170407,-0.00163791,-0.00130378,-0.000733613,0,0.000792307,0.00152078,0.00206352,0.00231896,0.00222384,0.00176631,0.000991789,0,-0.00106697,-0.0020443,-0.00276918,-0.00310702,-0.00297515,-0.0023598,-0.00132336,0,0.00142058,0.00271932,0.00368064,0.00412695,0.00394968,0.00313152,0.00175569,0,-0.00188451,-0.00360809,-0.00488533,-0.00548059,-0.00524887,-0.00416528,-0.0023378,0,0.00251638,0.00482627,0.00654774,0.0073621,0.00706868,0.0056253,0.00316721,0,-0.00343434,-0.00661506,-0.00901694,-0.0101911,0.0586698,0.0373629,0.0173007,0,-0.0133965,-0.0222193,-0.0263148,-0.0260258,-0.0221179,-0.0156601,0.00176631,0.00222384,0.00231896,0.00206352,0.00152078,0.000792307,0,-0.000733613,-0.00130378,-0.00163791,-0.00170407,-0.00151277,-0.00111214,-0.000577935,0,0.000532265,0.000943302,0.00118165,0.00122577,0.00108489,0.00079513,0.000411904,0,-0.000376929,-0.00066583,-0.000831324,-0.000859512,-0.000758204,-0.000553859,-0.000285975,0,0.00026001,0.000457869,0.000569955,0.000587586,0.000516923,0.00037
升级会员 