diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md
index 29d9e8e09..d594ad218 100644
--- a/CONTRIBUTORS.md
+++ b/CONTRIBUTORS.md
@@ -56,3 +56,4 @@ This file lists everyone, who contributed to this repo and wanted to show up her
 - Ishaan Verma
 - Delphi1024
 - ntindle
+- Ridham177
diff --git a/contents/convolutions/2d/2d.md b/contents/convolutions/2d/2d.md
index d026c8a52..7fc21aeb6 100644
--- a/contents/convolutions/2d/2d.md
+++ b/contents/convolutions/2d/2d.md
@@ -21,6 +21,8 @@ In code, a two-dimensional convolution might look like this:
 {% method %}
 {% sample lang="jl" %}
 [import:4-28, lang:"julia"](../code/julia/2d_convolution.jl)
+{% sample lang="py" %}
+[import:5-19, lang:"python"](../code/python/2d_convolution.py)
 {% endmethod %}
 
 This is very similar to what we have shown in previous sections; however, it essentially requires four iterable dimensions because we need to iterate through each axis of the output domain *and* the filter.
@@ -48,6 +50,8 @@ At this stage, it is important to write some code, so we will generate a simple
 {% method %}
 {% sample lang="jl" %}
 [import:30-47, lang:"julia"](../code/julia/2d_convolution.jl)
+{% sample lang="py" %}
+[import:21-33, lang:"python"](../code/python/2d_convolution.py)
 {% endmethod %}
 
 Though it is entirely possible to create a Gaussian kernel whose standard deviation is independent on the kernel size, we have decided to enforce a relation between the two in this chapter.
@@ -135,6 +139,8 @@ In code, the Sobel operator involves first finding the operators in $$x$$ and $$
 {% method %}
 {% sample lang="jl" %}
 [import:49-63, lang:"julia"](../code/julia/2d_convolution.jl)
+{% sample lang="py" %}
+[import:36-52, lang:"python"](../code/python/2d_convolution.py)
 {% endmethod %}
 
 With that, I believe we are at a good place to stop discussions on two-dimensional convolutions.
@@ -148,6 +154,8 @@ We have also added code to create the Gaussian kernel and Sobel operator and app
 {% method %}
 {% sample lang="jl" %}
 [import, lang:"julia"](../code/julia/2d_convolution.jl)
+{% sample lang="py" %}
+[import, lang:"python"](../code/python/2d_convolution.py)
 {% endmethod %}
 
 <script>
diff --git a/contents/convolutions/code/python/2d_convolution.py b/contents/convolutions/code/python/2d_convolution.py
new file mode 100644
index 000000000..48b79b611
--- /dev/null
+++ b/contents/convolutions/code/python/2d_convolution.py
@@ -0,0 +1,122 @@
+import numpy as np
+from contextlib import suppress
+
+
+def convolve_linear(signal, filter, output_size):
+    out = np.zeros(output_size)
+    sum = 0
+
+    for i in range(output_size[0]):
+        for j in range(output_size[1]):
+            for k in range(max(0, i-filter.shape[0]), i+1):
+                for l in range(max(0, j-filter.shape[1]), j+1):
+                    with suppress(IndexError):
+                        sum += signal[k, l] * filter[i-k, j-l]
+            out[i, j] = sum
+            sum = 0
+
+    return out
+
+
+def create_gaussian_kernel(kernel_size):
+    kernel = np.zeros((kernel_size, kernel_size))
+
+    # The center must be offset by 0.5 to find the correct index
+    center = kernel_size*0.5 + 0.5
+
+    sigma = np.sqrt(0.1*kernel_size)
+
+    def kernel_function(x, y):
+        return np.exp(-((x-center+1)**2 + (y-center+1)**2)/(2*sigma**2))
+
+    kernel = np.fromfunction(kernel_function, (kernel_size, kernel_size))
+    return kernel / np.linalg.norm(kernel)
+
+
+def create_sobel_operators():
+    Sx = np.dot([[1.0], [2.0], [1.0]], [[-1.0, 0.0, 1.0]]) / 9
+    Sy = np.dot([[-1.0], [0.0], [1.0]], [[1.0, 2.0, 1.0]]) / 9
+
+    return Sx, Sy
+
+def sum_matrix_dimensions(mat1, mat2):
+    return (mat1.shape[0] + mat2.shape[0], 
+            mat1.shape[1] + mat2.shape[1])
+
+def compute_sobel(signal):
+    Sx, Sy = create_sobel_operators()
+    
+    Gx = convolve_linear(signal, Sx, sum_matrix_dimensions(signal, Sx))
+    Gy = convolve_linear(signal, Sy, sum_matrix_dimensions(signal, Sy))
+
+    return np.sqrt(np.power(Gx, 2) + np.power(Gy, 2))
+
+
+def create_circle(image_resolution, grid_extents, radius):
+    out = np.zeros((image_resolution, image_resolution))
+
+    for i in range(image_resolution):
+        x_position = ((i * grid_extents / image_resolution)
+                      - 0.5 * grid_extents)
+        for j in range(image_resolution):
+            y_position = ((j * grid_extents / image_resolution)
+                          - 0.5 * grid_extents)
+            if x_position ** 2 + y_position ** 2 <= radius ** 2:
+                out[i, j] = 1.0
+
+    return out
+
+
+def main():
+
+    # Random distribution in x
+    x = np.random.rand(100, 100)
+
+    # Gaussian signals
+    def create_gaussian_signals(i, j):
+        return np.exp(-(((i-50)/100) ** 2 +
+                        ((j-50)/100) ** 2) / .01)
+    y = np.fromfunction(create_gaussian_signals, (100, 100))
+
+    # 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
+    full_linear_output = convolve_linear(x, y, sum_matrix_dimensions(x, y))
+
+    # simple boundaries
+    simple_linear_output = convolve_linear(x, y, x.shape)
+
+    np.savetxt("full_linear.dat", full_linear_output)
+    np.savetxt("simple_linear.dat", simple_linear_output)
+
+    # creating simple circle and 2 different Gaussian kernels
+    circle = create_circle(50, 2, 0.5)
+
+    circle = circle / np.linalg.norm(circle)
+
+    small_kernel = create_gaussian_kernel(3)
+    large_kernel = create_gaussian_kernel(25)
+
+    small_kernel_output = convolve_linear(circle, small_kernel,
+                                          sum_matrix_dimensions(circle,
+                                                                small_kernel))
+
+    large_kernel_output = convolve_linear(circle, large_kernel,
+                                          sum_matrix_dimensions(circle,
+                                                                large_kernel))
+
+    np.savetxt("small_kernel.dat", small_kernel_output)
+    np.savetxt("large_kernel.dat", large_kernel_output)
+
+    circle = create_circle(50, 2, 0.5)
+
+    # Normalization
+    circle = circle / np.linalg.norm(circle)
+
+    # using the circle for sobel operations as well
+    sobel_output = compute_sobel(circle)
+
+    np.savetxt("sobel_output.dat", sobel_output)
+