Writing a program to convolve two inputs [closed]

L = Length[h];

M = Length[x];

y = Table[0,M - L + 1];

convolve[inputx_, inputh_] :=Module[{k,n},

Do[

x = Table[inputx];

h = Table[inputh];

L = Length[h];

M = Length[x];

For[k = 1, k < L + 1, k++, y[[n]] += h[[k]]*x[[L - k + n]]], {n, M - L + 1}]];

convolve[{1, 2,  2, -1}, {1, -1}]
y


It returns 1,0,-3 but thats wrong. It should be 1,1,0,-3,1 why is it leaving out the 1's on the outsides?

• (1) You cannot define the array y until you know its length. Which is not known until the x and h inputs are known. (2) A construct such as Table[{1,2,2,-1}] basically makes no sense in Mathematica. Check documentation for further details. (3) If this code works at all, then y will have length 4-2+1=3. So a result with 5 elements is impossible. Mar 16 '20 at 23:01

There are many ways to convolve (filter) data x with a kernel h: convolution, correlation, using the frequency domain method, and directly in the time domain. The only difference (other than numerical factors) is in the way edge conditions are handled with padding. Here is the setup:

h = {1, -1, 2, -2, 3, -3};
x = {1, 2, 3, 4, 5, 6, -5, -4, -3, -2, -1};
n = Length[x] + Length[h] - 1;


In the convolution method, the kernel h is thought of as the impulse response of a linear time-invariant system and the x is thought of as the input to that system. The convolution yConv is then the output of the system.

yConv = ListConvolve[h, x, {1, 1}, 0];


In the correlation method, the kernel h is thought of as a marker or mask and x is thought of as the data that is to be examined. The correlation yCorr is then how much like x the kernel is at each place in the sequence.

yCorr = ListCorrelate[Reverse[h], x, {-1, -1}, 0];


The Fourier method exploits the fact from Fourier Transforms that the product of the transfoms is equal to the convolution of the time domain signals. The following calculate the Fourier transform of h (ffth) and the Fourier transform of x (fftx), after padding to the same length. The element-by-element product is then inverse transformed, giving yFourier. which is numerically the same as the above methods.

ffth = Fourier[PadRight[h, n], FourierParameters -> {1, -1}];
fftx =  Fourier[PadRight[x, n], FourierParameters -> {1, -1}];
yFourier = InverseFourier[ffth fftx, FourierParameters -> {1, -1}]


In the time-domain method, the output of the system with impulse response h is calculated once for each time k, as the input takes on all values in x.

z = PadLeft[x, n];
yTim = ConstantArray[0, Length[x]];
Do[yTim[[k]] = Total[Reverse[h] z[[k ;; k + Length[h] - 1]]];, {k, 1, Length[x]}]


Here's a time domain version that's like one might program it in Java or C. Normally one would truncate the initial string of zeros.

z = PadLeft[x, n];
yJav = ConstantArray[0, n + 1];
Do[
Do[yJav[[k]] = yJav[[k]] + h[[j]] z[[k - j]];, {j, 1, Length[h]}];,
{k, Length[h] + 1, Length[x] + Length[h]}];


Without ListConvolve:

myListConvolve[kernel_, list_][n_] :=
Dot[Reverse@kernel, Take[RotateLeft[list, n - 1], Length[kernel]]]


example:

myListConvolve[{x, y}, {a, b, c, d, e, f}] /@ Range[5]


myListConvolve will be faster than your procedural program (maybe 20 or 30 x) but slower than ListConvolve.

When working out your solution, do not hesitate to use symbols instead of numbers, so that you can see more clearly what is going on.

Have a look at the documentation: ListConvolve does this (and may be easily 100 x faster than your code).

• I have used the function in the past but I am required to use this method
– Noah
Mar 16 '20 at 21:39
• Sorry I answered too fast. Are you required to use a procedural program for your exercise? Mar 16 '20 at 22:10