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