From 5152de72b2f499ea4e5800776704198aa9cc3b02 Mon Sep 17 00:00:00 2001 From: Ayman Lafaz Date: Sun, 10 Oct 2021 19:31:58 +0100 Subject: [PATCH 1/4] added 1d convolution implementation in python --- contents/convolutions/1d/1d.md | 14 ++++++ .../1d/code/python/1d_convolution.py | 49 +++++++++++++++++++ 2 files changed, 63 insertions(+) create mode 100644 contents/convolutions/1d/code/python/1d_convolution.py diff --git a/contents/convolutions/1d/1d.md b/contents/convolutions/1d/1d.md index d5d77652b..e3f9a3770 100644 --- a/contents/convolutions/1d/1d.md +++ b/contents/convolutions/1d/1d.md @@ -56,6 +56,8 @@ With this in mind, we can almost directly transcribe the discrete equation into [import:27-46, lang:"julia"](code/julia/1d_convolution.jl) {% sample lang="cs" %} [import:63-84, lang:"csharp"](code/csharp/1DConvolution.cs) +{% sample lang="py" %} +[import:18-27, lang:"python"](code/python/1d_convolution.py) {% endmethod %} The easiest way to reason about this code is to read it as you might read a textbook. @@ -189,6 +191,8 @@ Here it is again for clarity: [import:27-46, lang:"julia"](code/julia/1d_convolution.jl) {% sample lang="cs" %} [import:63-84, lang:"csharp"](code/csharp/1DConvolution.cs) +{% sample lang="py" %} +[import:18-27, lang:"python"](code/python/1d_convolution.py) {% endmethod %} Here, the main difference between the bounded and unbounded versions is that the output array size is smaller in the bounded case. @@ -199,6 +203,8 @@ For an unbounded convolution, the function would be called with a the output arr [import:58-59, lang:"julia"](code/julia/1d_convolution.jl) {% sample lang="cs" %} [import:96-97, lang:"csharp"](code/csharp/1DConvolution.cs) +{% sample lang="py" %} +[import:37-38, lang:"python"](code/python/1d_convolution.py) {% endmethod %} On the other hand, the bounded call would set the output array size to simply be the length of the signal @@ -208,6 +214,8 @@ On the other hand, the bounded call would set the output array size to simply be [import:61-62, lang:"julia"](code/julia/1d_convolution.jl) {% sample lang="cs" %} [import:98-99, lang:"csharp"](code/csharp/1DConvolution.cs) +{% sample lang="py" %} +[import:40-41, lang:"python"](code/python/1d_convolution.py) {% endmethod %} Finally, as we mentioned before, it is possible to center bounded convolutions by changing the location where we calculate the each point along the filter. @@ -218,6 +226,8 @@ This can be done by modifying the following line: [import:35-35, lang:"julia"](code/julia/1d_convolution.jl) {% sample lang="cs" %} [import:71-71, lang:"csharp"](code/csharp/1DConvolution.cs) +{% sample lang="py" %} +[import:22-22, lang:"python"](code/python/1d_convolution.py) {% endmethod %} Here, `j` counts from `i-length(filter)` to `i`. @@ -252,6 +262,8 @@ In code, this typically amounts to using some form of modulus operation, as show [import:4-25, lang:"julia"](code/julia/1d_convolution.jl) {% sample lang="cs" %} [import:38-61, lang:"csharp"](code/csharp/1DConvolution.cs) +{% sample lang="py" %} +[import:5-15, lang:"python"](code/python/1d_convolution.py) {% endmethod %} This is essentially the same as before, except for the modulus operations, which allow us to work on a periodic domain. @@ -269,6 +281,8 @@ For the code associated with this chapter, we have used the convolution to gener [import, lang:"julia"](code/julia/1d_convolution.jl) {% sample lang="cs" %} [import, lang:"csharp"](code/csharp/1DConvolution.cs) +{% sample lang="py" %} +[import, lang:"python"](code/python/1d_convolution.py) {% endmethod %} At a test case, we have chosen to use two sawtooth functions, which should produce the following images: diff --git a/contents/convolutions/1d/code/python/1d_convolution.py b/contents/convolutions/1d/code/python/1d_convolution.py new file mode 100644 index 000000000..f48e5f0af --- /dev/null +++ b/contents/convolutions/1d/code/python/1d_convolution.py @@ -0,0 +1,49 @@ +from scipy import fft +import numpy as np + + +def convolve_cyclic(signal, filter_array): + output_size = max(len(signal), len(filter_array)) + out = np.zeros(output_size) + s = 0 + for i in range(output_size): + for j in range(output_size): + if (i - j)% output_size < len(filter_array): + s += signal[(j - 1) % output_size] * filter_array[(i - j) % output_size] + out[i] = s + s = 0 + return out + + +def convolve_linear(signal, filter_array, output_size): + out = np.zeros(output_size) + s = 0 + for i in range(output_size): + for j in range(max(0, i - len(filter_array)), i + 1): + if j < len(signal) and (i - j) < len(filter_array): + s += signal[j] * filter_array[i - j] + out[i] = s + s = 0 + return out + +# sawtooth functions for x and y +x = [float(i)/200 for i in range(1, 201)] +y = [float(i)/200 for i in range(1, 201)] + +# Normalization is not strictly necessary, but good practice +x = x/np.linalg.norm(x) +y = y/np.linalg.norm(x) + +# full convolution, output will be the size of x + y - 1 +full_linear_output = convolve_linear(x, y, len(x) + len(y) - 1) + +# simple boundaries +simple_linear_output = convolve_linear(x, y, len(x)) + +# cyclic convolution +cyclic_output = convolve_cyclic(x, y) + +# outputting convolutions to different files for plotting in external code +np.savetxt('full_linear.dat', full_linear_output) +np.savetxt('simple_linear.dat', simple_linear_output) +np.savetxt('cyclic.dat', cyclic_output) From 2a00d87577b5c101a9bb0d3bd870bea5e66a953b Mon Sep 17 00:00:00 2001 From: Ayman Lafaz Date: Mon, 11 Oct 2021 17:16:10 +0100 Subject: [PATCH 2/4] fixed some mistakes in the code so it outputs correct results --- .../convolutions/1d/code/python/1d_convolution.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/contents/convolutions/1d/code/python/1d_convolution.py b/contents/convolutions/1d/code/python/1d_convolution.py index f48e5f0af..53a2ffadf 100644 --- a/contents/convolutions/1d/code/python/1d_convolution.py +++ b/contents/convolutions/1d/code/python/1d_convolution.py @@ -1,6 +1,6 @@ -from scipy import fft import numpy as np +def mod1(x, y): return ((x%y) + y) % y def convolve_cyclic(signal, filter_array): output_size = max(len(signal), len(filter_array)) @@ -8,8 +8,8 @@ def convolve_cyclic(signal, filter_array): s = 0 for i in range(output_size): for j in range(output_size): - if (i - j)% output_size < len(filter_array): - s += signal[(j - 1) % output_size] * filter_array[(i - j) % output_size] + if(mod1(i - j, output_size) < len(filter_array)): + s += signal[mod1(j - 1, output_size)] * filter_array[mod1(i - j, output_size)] out[i] = s s = 0 return out @@ -27,15 +27,15 @@ def convolve_linear(signal, filter_array, output_size): return out # sawtooth functions for x and y -x = [float(i)/200 for i in range(1, 201)] -y = [float(i)/200 for i in range(1, 201)] +x = [float(i + 1)/200 for i in range(200)] +y = [float(i + 1)/200 for i in range(200)] # Normalization is not strictly necessary, but good practice x = x/np.linalg.norm(x) -y = y/np.linalg.norm(x) +y = y/np.linalg.norm(y) # full convolution, output will be the size of x + y - 1 -full_linear_output = convolve_linear(x, y, len(x) + len(y) - 1) +full_linear_output = convolve_linear(x, y, len(x) + len(y) -1) # simple boundaries simple_linear_output = convolve_linear(x, y, len(x)) From c83d6449097af5f337f111dbdfacc46bebabadbd Mon Sep 17 00:00:00 2001 From: Ayman Lafaz Date: Mon, 11 Oct 2021 17:23:48 +0100 Subject: [PATCH 3/4] making the code look better --- contents/convolutions/1d/code/python/1d_convolution.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/contents/convolutions/1d/code/python/1d_convolution.py b/contents/convolutions/1d/code/python/1d_convolution.py index 53a2ffadf..88ef0b579 100644 --- a/contents/convolutions/1d/code/python/1d_convolution.py +++ b/contents/convolutions/1d/code/python/1d_convolution.py @@ -1,6 +1,6 @@ import numpy as np -def mod1(x, y): return ((x%y) + y) % y +def mod1(x, y): return ((x % y) + y) % y def convolve_cyclic(signal, filter_array): output_size = max(len(signal), len(filter_array)) @@ -31,11 +31,11 @@ def convolve_linear(signal, filter_array, output_size): y = [float(i + 1)/200 for i in range(200)] # Normalization is not strictly necessary, but good practice -x = x/np.linalg.norm(x) -y = y/np.linalg.norm(y) +x /= np.linalg.norm(x) +y /= np.linalg.norm(y) # full convolution, output will be the size of x + y - 1 -full_linear_output = convolve_linear(x, y, len(x) + len(y) -1) +full_linear_output = convolve_linear(x, y, len(x) + len(y) - 1) # simple boundaries simple_linear_output = convolve_linear(x, y, len(x)) From cba2d7026138ec320c18e6a4128bd684daa140d4 Mon Sep 17 00:00:00 2001 From: Ayman Lafaz Date: Tue, 12 Oct 2021 01:17:42 +0100 Subject: [PATCH 4/4] spacing code properly for readability --- contents/convolutions/1d/code/python/1d_convolution.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/contents/convolutions/1d/code/python/1d_convolution.py b/contents/convolutions/1d/code/python/1d_convolution.py index 88ef0b579..e77e68d09 100644 --- a/contents/convolutions/1d/code/python/1d_convolution.py +++ b/contents/convolutions/1d/code/python/1d_convolution.py @@ -6,24 +6,28 @@ def convolve_cyclic(signal, filter_array): output_size = max(len(signal), len(filter_array)) out = np.zeros(output_size) s = 0 + for i in range(output_size): for j in range(output_size): if(mod1(i - j, output_size) < len(filter_array)): s += signal[mod1(j - 1, output_size)] * filter_array[mod1(i - j, output_size)] out[i] = s s = 0 + return out def convolve_linear(signal, filter_array, output_size): out = np.zeros(output_size) s = 0 + for i in range(output_size): for j in range(max(0, i - len(filter_array)), i + 1): if j < len(signal) and (i - j) < len(filter_array): s += signal[j] * filter_array[i - j] out[i] = s s = 0 + return out # sawtooth functions for x and y