Re: Discrete Convolution

*To*: mathgroup at smc.vnet.net*Subject*: [mg19023] Re: [mg18943] Discrete Convolution*From*: "Wolf, Hartmut" <hwolf at debis.com>*Date*: Tue, 3 Aug 1999 13:44:55 -0400*Organization*: debis Systemhaus*References*: <199907300533.BAA18638@smc.vnet.net.>*Sender*: owner-wri-mathgroup at wolfram.com

Dear Alister Alister McAlister schrieb: > > I want a function that mimics Matlab's "conv" function for doing a discrete > convolution of two lists. > If you just want to emulate your Matlab's programs ... ok, you have to decide at which level to do that most economically. But perhaps you better switch to Mathematica style: she allows e.g. to multiply polynomials directly in a more explicit representation. So there is no reason to further use: > CONV Convolution and polynomial multiplication. > C = CONV(A, B) convolves vectors A and B. The resulting > vector is length LENGTH(A)+LENGTH(B)-1. > If A and B are vectors of polynomial coefficients, convolving > them is equivalent to multiplying the two polynomials. > > I wrote the following, but is there a way of either of > (1) speeding up the code by changing the algorithm ... > ignoring simple things like the multiple evaluations > of Length and so forth which I have left > in only for what I hope is clarity; or > (2) Using a built in function (possibly connected with polynomials) to do > the same thing? > > Mark R Diamond > No spam email: markd at psy dot uwa dot edu dot au > -------------------------------------------------------- > > convolve[a_List,b_List]:=Module[ > { > (* reverse one of the lists prior to the convolution *) > ra=Reverse[a], > > (* A variable that collects the indices of lists ra and b, > respectively *) > (* that will be Dot[ ]-ed together. *) > indices > }, > > (* Create the table of indices *) > indices=Table[ > { > { > Max[Length[a]+1-i,1], > Min[Length[a],Length[a]+Length[b]-i] > }, > { > Max[1,i-Length[a]+1],Min[Length[b],i] > } > }, > {i,Length[a]+Length[b]-1} > ]; > > (* Create a list of the appropriate pairs of dot products *) > Map[(Take[ra,#[[1]] ].Take[ b,#[[2]] ])&, indices] > ] /; (VectorQ[a,NumberQ]\[And]VectorQ[b,NumberQ]) Go ahead, that's ok, the program looks somewhat like: In[1]:= convolve[a_List, b_List] := With[{n = Length[a], m = Length[b]}, Table[Take[a, {Max[1, i - m + 1], Min[i, n]}].Reverse[ Take[b, {Max[1, i - n + 1], Min[i, m]}]], {i, m + n - 1}]] Now Mathematica can do things symbolically: In[8]:= coeffA = Array[a, {5}, {0}] Out[8]= {a[0], a[1], a[2], a[3], a[4]} In[9]:= coeffB = Array[b, {7}, {0}] Out[9]= {b[0], b[1], b[2], b[3], b[4], b[5], b[6]} In[10]:= convolve[coeffA, coeffB] Out[10]= {a[0] b[0], a[1] b[0] + a[0] b[1], a[2] b[0] + a[1] b[1] + a[0] b[2], a[3] b[0] + a[2] b[1] + a[1] b[2] + a[0] b[3], a[4] b[0] + a[3] b[1] + a[2] b[2] + a[1] b[3] + a[0] b[4], a[4] b[1] + a[3] b[2] + a[2] b[3] + a[1] b[4] + a[0] b[5], a[4] b[2] + a[3] b[3] + a[2] b[4] + a[1] b[5] + a[0] b[6], a[4] b[3] + a[3] b[4] + a[2] b[5] + a[1] b[6], a[4] b[4] + a[3] b[5] + a[2] b[6], a[4] b[5] + a[3] b[6], a[4] b[6]} In[11]:= convolve[coeffB, coeffA] == convolve[coeffA, coeffB] Out[11]= True However, to do that more Mathematica-styled, first load the Package "Streams" from Roman E. Maeder (you'll find it on his web pages, or in his book): In[2]:= << "c:/my Workplace/home/\ Mathematica/3.0/AddOns/Applications/TheMathematicaProgrammer/Streams.m" In[3]:= ?Accumulate "Accumulate[f, s] returns the stream {s1, f[s1, s2], f[s1, f[s2, s3]],...}." (* You have to make an obvious slight adoption for version 4 in the notebook, because the obsolete Accumulate no longer is defined there. *) In[4]:= X := MakeStream[1, x X] This is a "list" of all (infinite many) exponents of x, starting with x^0 (==1). Now you can easily generate the polynomials from lists of coefficients: In[12]:= coeffA.Take[X, Length[coeffA]] Out[12]= a[0] + x a[1] + x^2 a[2] + x^3 a[3] + x^4 a[4] (If you take a minute, you may incorporate Dotting a List with a Stream into the package -- R.E.M. left that out to be done by us -- then this simply might look somewhat like coeffA.X ) Or (for our demonstration here) do In[5]:= poly[A] = With[{l = 5, name = a}, Array[name, {l}, {0}].Take[X, l]] Out[5]= a[0] + x a[1] + x^2 a[2] + x^3 a[3] + x^4 a[4] In[6]:= poly[B] = With[{l = 7, name = b}, Array[name, {l}, {0}].Take[X, l]] Out[6]= b[0] + x b[1] + x^2 b[2] + x^3 b[3] + x^4 b[4] + x^5 b[5] + x^6 b[6] Now the Mathematica function you searched for simply is In[7]:= CoefficientList[poly[A] poly[B], x] Out[7]= {a[0] b[0], a[1] b[0] + a[0] b[1], a[2] b[0] + a[1] b[1] + a[0] b[2], a[3] b[0] + a[2] b[1] + a[1] b[2] + a[0] b[3], a[4] b[0] + a[3] b[1] + a[2] b[2] + a[1] b[3] + a[0] b[4], a[4] b[1] + a[3] b[2] + a[2] b[3] + a[1] b[4] + a[0] b[5], a[4] b[2] + a[3] b[3] + a[2] b[4] + a[1] b[5] + a[0] b[6], a[4] b[3] + a[3] b[4] + a[2] b[5] + a[1] b[6], a[4] b[4] + a[3] b[5] + a[2] b[6], a[4] b[5] + a[3] b[6], a[4] b[6]} To refabricate a "true" polynomial you may use In[12] above, or do that directly: In[16]:= Collect[poly[A] poly[B], x] Out[16]= a[0] b[0] + x (a[1]\ b[0] + a[0] b[1]) + x^2 (a[2] b[0] + a[1] b[1] + a[0] b[2]) + x^3 (a[3] b[0] + a[2] b[1] + a[1] b[2] + a[0] b[3]) + x^4 (a[4] b[0] + a[3] b[1] + a[2] b[2] + a[1] b[3] + a[0] b[4]) + x^5 (a[4] b[1] + a[3] b[2] + a[2] b[3] + a[1] b[4] + a[0] b[5]) + x^10 a[4] b[6] + x^6 (a[4] b[2] + a[3] b[3] + a[2] b[4] + a[1] b[5] + a[0] b[6]) + x^7 (a[4] b[3] + a[3] b[4] + a[2] b[5] + a[1] b[6]) + x^8 (a[4] b[4] + a[3] b[5] + a[2] b[6]) + x^9 (a[4] b[5] + a[3] b[6]) (Perhaps you won't like the ordering of the x^10 term,... but addition is commutative ;-) Kind regards, hw