MathGroup Archive 1999

[Date Index] [Thread Index] [Author Index]

Search the Archive

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



  • Prev by Date: Re: Removing Graphics3D polygon edges
  • Next by Date: Re: Wheel mouse under X.
  • Previous by thread: Re: Discrete Convolution
  • Next by thread: ListPlot: corner problem ??