Implement conv in python #408

jasmaa · 2018-10-01T16:51:28Z

Note: Discrete convolution implemented differently because I couldn't get the previous implementations to get the same answer as the fft ones.

jiegillet · 2018-10-02T04:14:48Z

That's normal, the FFT method does not work the same way.

jiegillet

A few things to change. I think if you want to use numpy that's fine, but then we should use more of its tools to make your code cleaner. But maybe for the purpose of the AAA, not using is is better? I'm not sure.

jiegillet · 2018-10-02T04:17:09Z

contents/convolutions/code/python/conv.py

@@ -0,0 +1,46 @@
+import numpy as np


You're only using numpy to make examples, I don't think it's worth it...

jiegillet · 2018-10-02T04:20:36Z

contents/convolutions/code/python/conv.py

+    """Discrete convolution by definition"""
+
+    out = []
+    n = min(len(signal1), len(signal2))


You should end up with len(signal1) + len(signal2) - 1 points, not min(len(signal1), len(signal2))

jiegillet · 2018-10-02T04:23:33Z

contents/convolutions/code/python/conv.py

+    for i in range(n):
+        s = 0
+        for j in range(n):
+            s += signal1[j] * signal2[i - j]


You're not overlapping all the points. Look at the Julia code for reference.
Don't try to mimic the FFT method, they have different assumptions.

Thanks for the feedback!
I'm not that familiar with convolutions so please bear with me:

The discrete convolution from the definition is as listed:

f*g[n] = sum[-inf,inf][ f[m]*g[n-m] ]

In the code, I am going for each element of f*g using the outer for loop.

Taking the [-inf,inf] to effectively be the full length of the signal arrays,
for each element in f*g, I sum the values f[m]*g[n-m] where n and m are represented
by i in the outer for loop and j in the inner for loop respectively.
This j should go the length of the signal as m is in the range, [-inf,inf].

The code follows the discrete convolution definition to my knowledge and,
additionally, also appears to give the same answer as the fft version,
so I am unsure as to why it is not valid.

I will experiment with the Julia code when I can get a chance but a cursory
look at the code shows the length of the outputted signal to be the sum of
the two inputted signals. Since the outer for loop runs through this combined
length, isn't there a chance of an array-out-of-bounds when indexing the smaller
signal2 with i-j when i is big and j is small?

FFT also appears to output a signal of equal length to the input signal, so
I am unclear if there is ever a case where the fft version of the convolution
could give back the same output signal described by the Julia implementation.

I will try to find more information since there is likely a misunderstanding on
my part. Meanwhile, it would be helpful if some of these questions can be cleared up.

The trick is to simulate a [-inf,inf] range by adding an infinite number of zeroes to the left and right of both signals, convolute them and clip the zeros at the end. For example, I want to convolve a=[1,2,3,4] and b=[8,9]. So if I flip b and slide it across a, I have:

0 0 1 2 3 4 0 0 9 8 0 0 0 0 0 0 => sum of product = 0 0 0 1 2 3 4 0 0 0 9 8 0 0 0 0 0 => sum of product = 0*9 + 1*8 = 8 0 0 1 2 3 4 0 0 0 0 9 8 0 0 0 0 => sum of product = 1*9 + 2*8 = 25 0 0 1 2 3 4 0 0 0 0 0 9 8 0 0 0 => sum of product = 2*9 + 3*8 = 42 0 0 1 2 3 4 0 0 0 0 0 0 9 8 0 0 => sum of product = 3*9 + 4*8 = 59 0 0 1 2 3 4 0 0 0 0 0 0 0 9 8 0 => sum of product = 4*9 + 0*8 = 36 0 0 1 2 3 4 0 0 0 0 0 0 0 0 9 8 => sum of product = 0

so the output should be [8, 25, 42, 59, 36], which is of length 5 = len(a)+len(b)-1.
This is the best way to mimic the continuous case.

The FFT method reproduces the same result only if the signals you want to convolve are infinite and periodic. If you use the FFT method on finite signals it will assume that both signals repeat left and right forever, so if I did the same process above, we would have something like:

. . 1 2 3 4 . . 9 8 9 8 9 8 9 8 => sum of product = 1*9 + 2*8 + 3*9 + 4*8 = A . . 1 2 3 4 . . 8 9 8 9 8 9 8 9 => sum of product = 1*8 + 2*9 + 3*8 + 4*9 = B . . 1 2 3 4 . . 9 8 9 8 9 8 9 8 => sum of product = 1*9 + 2*8 + 3*9 + 4*8 = A . . 1 2 3 4 . . 8 9 8 9 8 9 8 9 => sum of product = 1*8 + 2*9 + 3*8 + 4*9 = B

and the output would be [B, A] (or maybe [A,B]?) of length 2 = min(len(a), len(b)).

But it gets a bit more messy when you actually use the FFT when you have signals of different lengths, you need to either cut the longest signal (which you have done, and is actually a valid choice) or pad the shorter signal with zeros (that would be my preference, but it's not more valid).

Depending on your data you might chose different methods. The FFT is a faster algorithm but comes with different assumptions on the data.

Understood. Thanks for the explanation!

jiegillet · 2018-10-02T04:25:35Z

contents/convolutions/code/python/conv.py

+    fft_s2 = fft(signal2)
+    out = []
+
+    for i in range(min(len(signal1), len(signal2))):


I'm not sure if there is an official method in the case the signals have different lengths. I think instead of cutting the longer signal, I would add zeros to the shorter one until they have the same length. That's the usual method when using FFTs.

This was based off of what I could gather from the C implementation on AAA, but I'll try to change it.

According to Wikipedia, both ways of doing it are valid, the results will a bit look different. As long as you're making an informed decision it's all good.

jiegillet · 2018-10-02T04:38:30Z

contents/convolutions/code/python/conv.py

+    # Example convolution with sin and cos
+    x = np.linspace(0, 1, 5)
+    s1 = np.sin(x)
+    s2 = np.cos(x)


I think you can get away with something like s1 = [ math.sin(x) for x in range(5)] if we want to skip numpy.

jiegillet

The output is correct but you can simplify your code by a lot!

jiegillet · 2018-10-03T02:42:24Z

contents/convolutions/code/python/conv.py

+    """Discrete convolution by definition"""
+
+    n = len(signal1) + len(signal2)
+    signal1 = signal1.copy()


You don't need to modify the signals if your loops are done right. Lines 8 to 17 can be deleted without any effect on the algorithm.

jiegillet · 2018-10-03T02:44:10Z

contents/convolutions/code/python/conv.py

+    out = []
+
+    for i in range(n):
+        s = complex(0)


Why complex? just 0 is fine, if the data is complex it will get updated anyway.

jiegillet · 2018-10-03T02:45:22Z

contents/convolutions/code/python/conv.py

+
+        out.append(s)
+
+    return out[1:]


return out.

jiegillet · 2018-10-03T02:48:33Z

contents/convolutions/code/python/conv.py

+def conv(signal1, signal2):
+    """Discrete convolution by definition"""
+
+    n = len(signal1) + len(signal2)


Define n = len(signal1) + len(signal2) -1 that's the most natural value

jiegillet · 2018-10-03T02:49:38Z

contents/convolutions/code/python/conv.py

+    for i in range(n):
+        s = complex(0)
+
+        for j in range(i):


change to range(i+1) this way all the i-1 become i and it's clearer.

jiegillet · 2018-10-03T02:49:56Z

contents/convolutions/code/python/conv.py

+        s = complex(0)
+
+        for j in range(i):
+            if j < len(signal1) and i - j - 1  < len(signal2):


i - j - 1 => i - j

jiegillet · 2018-10-03T02:50:26Z

contents/convolutions/code/python/conv.py

+
+        for j in range(i):
+            if j < len(signal1) and i - j - 1  < len(signal2):
+                s += signal1[j] * signal2[i - j - 1]


i - j - 1 => i - j

jiegillet

Very nice! Thank you for the quick updates and your great contribution.

* Implement conv in python * Fix conv func * fix range * remove numpy usage * Fix discrete convolution * Revert julia file * Add padding for fft * Make requested changes

jasmaa added 3 commits October 1, 2018 12:45

Implement conv in python

82ed04f

Fix conv func

260c750

fix range

d928569

june128 added Implementation This provides an implementation for an algorithm. (Code and maybe md files are edited.) Hacktoberfest The label for all Hacktoberfest related things! labels Oct 1, 2018

jiegillet requested changes Oct 2, 2018

View reviewed changes

jiegillet self-assigned this Oct 2, 2018

jasmaa added 4 commits October 2, 2018 07:47

remove numpy usage

bbba44c

Fix discrete convolution

a12dbb2

Revert julia file

18c144c

Add padding for fft

5a39c7e

jiegillet requested changes Oct 3, 2018

View reviewed changes

Make requested changes

28549a9

jiegillet approved these changes Oct 3, 2018

View reviewed changes

jiegillet merged commit 0fb0325 into algorithm-archivists:master Oct 3, 2018

Uh oh!

Implement conv in python #408

Implement conv in python #408

Uh oh!

Conversation

jasmaa commented Oct 1, 2018

Uh oh!

jiegillet commented Oct 2, 2018

Uh oh!

jiegillet left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jiegillet left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jiegillet left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!