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..e77e68d09
--- /dev/null
+++ b/contents/convolutions/1d/code/python/1d_convolution.py
@@ -0,0 +1,53 @@
+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))
+    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
+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 /= 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)
+
+# 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)