diff --git a/CMakeLists.txt b/CMakeLists.txt
index ec59879..5c25fb7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -11,6 +11,7 @@ include_directories(${CUDAToolkit_INCLUDE_DIRS})
 set(LIBRARY_SOURCES
     src/utils/cuda_helper.cu
     src/kernels/activations.cu
+    src/kernels/convolution.cu
     src/kernels/padding.cu
     src/kernels/matrix_math.cu
     src/layers/dense.cu
diff --git a/src/layers/conv2d.cu b/src/layers/conv2d.cu
index 8277c20..f05b069 100644
--- a/src/layers/conv2d.cu
+++ b/src/layers/conv2d.cu
@@ -3,6 +3,7 @@
 
 #include "activations.cuh"
 #include "conv2d.cuh"
+#include "convolution.cuh"
 #include "cuda_helper.cuh"
 #include "padding.cuh"
 
@@ -31,24 +32,22 @@ Layers::Conv2d::Conv2d(
         outputSize  = (inputSize - kernelSize) / stride + 1;
     }
 
-    kernels.resize(kernelSize * kernelSize);
+    kernels.resize(kernelSize * kernelSize * numFilters);
     initializeKernels();
 
     d_kernels = nullptr;
 
     CUDA_CHECK(
-        cudaMalloc((void**)&d_kernels, sizeof(float) * kernelSize * kernelSize)
+        cudaMalloc((void**)&d_kernels, sizeof(float) * kernelSize * kernelSize * numFilters)
     );
     toCuda();
 
     d_padded = nullptr;
 
-    if (paddingSize > 0) {
-        CUDA_CHECK(cudaMalloc(
-            (void**)&d_padded, sizeof(float) * (inputSize + 2 * paddingSize) *
-                                   (inputSize + 2 * paddingSize) * inputChannels
-        ));
-    }
+    CUDA_CHECK(cudaMalloc(
+        (void**)&d_padded, sizeof(float) * (inputSize + 2 * paddingSize) *
+                                (inputSize + 2 * paddingSize) * inputChannels
+    ));
 }
 
 Layers::Conv2d::~Conv2d() {
@@ -67,22 +66,24 @@ void Layers::Conv2d::setKernels(const std::vector<float>& kernels_input) {
 
 void Layers::Conv2d::toCuda() {
     CUDA_CHECK(cudaMemcpy(
-        d_kernels, kernels.data(), sizeof(float) * kernelSize * kernelSize,
+        d_kernels, kernels.data(), sizeof(float) * kernelSize * kernelSize * numFilters,
         cudaMemcpyHostToDevice
     ));
 }
 
 void Layers::Conv2d::forward(const float* d_input, float* d_output) {
-    // Padd input
-    int THREADS_PER_BLOCK = 256;
-    int BLOCKS =
-        (outputSize * outputSize * inputChannels) / THREADS_PER_BLOCK + 1;
+    // Pad input
+    int THREADS_PER_BLOCK =  (inputSize + 2 * paddingSize) * (inputSize + 2 * paddingSize) * inputChannels;
 
-    pad_matrix_kernel<<<BLOCKS, THREADS_PER_BLOCK>>>(
+    pad_matrix_kernel<<<1, THREADS_PER_BLOCK>>>(
         d_input, d_padded, inputSize, inputSize, inputChannels, paddingSize
     );
 
-    // TODO: Implement 2D convolution
+    // Convolve
+    THREADS_PER_BLOCK = outputSize * outputSize * numFilters;
+    convolution_kernel<<<1, THREADS_PER_BLOCK>>>(
+        d_padded, d_kernels, d_output, inputSize + (2 * paddingSize), inputChannels, kernelSize, stride, numFilters, outputSize
+    );
 }
 
 /*
@@ -101,8 +102,6 @@ void Layers::Conv2d::host_conv(const float* input, float* output) {
                 
                 float sum = 0.0f;
 
-                // std::cout << "f: " << f << ", i: " << i << ", j: " << j << std::endl;
-
                 // Iterate over kernel and input matrix
                 for (int k = 0; k < kernelSize; k++) {
                     for (int l = 0; l < kernelSize; l++) {
@@ -111,15 +110,11 @@ void Layers::Conv2d::host_conv(const float* input, float* output) {
                             int kernelIndex = k * (kernelSize * inputChannels * numFilters) + l * (inputChannels * numFilters) + c * (numFilters) + f;
                             int inputIndex  = (i * stride + k) * (inputSize * inputChannels) + (j * stride + l) * (inputChannels) + c;
 
-                            // std::cout << "kernelIndex: " << kernelIndex << ", kernel value: " << kernels[kernelIndex] << ", inputIndex: " << inputIndex << ", input value: " << input[inputIndex] << std::endl;
-
                             sum += kernels[kernelIndex] * input[inputIndex];
                         }                      
                     }
                 }
 
-                // std::cout << "sum: " << sum << std::endl;
-
                 output[i * (outputSize * numFilters) + j * (numFilters) + f] = sum;
             }
         }
diff --git a/test/layers/test_conv2d.cu b/test/layers/test_conv2d.cu
index e4b675a..32313a0 100644
--- a/test/layers/test_conv2d.cu
+++ b/test/layers/test_conv2d.cu
@@ -5,7 +5,13 @@
 
 #include "conv2d.cuh"
 
-TEST(Conv2dTest, SimpleExample) {
+class Conv2dTest : public::testing::Test {
+    protected:
+        cudaError_t    cudaStatus;
+};
+
+
+TEST_F(Conv2dTest, SimpleExample) {
 
     int inputSize = 4;
     int inputChannels = 1;
@@ -38,18 +44,38 @@ TEST(Conv2dTest, SimpleExample) {
         1.0f, 2.0f, 3.0f, 4.0f,
     };
 
+    float* d_input;
+    float* d_output;
+
     conv2d.setKernels(kernels);
 
-    
-    std::vector<float> output(outputSize * outputSize * numFilters);
+            // Allocate device memory
+    cudaStatus = cudaMalloc((void**)&d_input, sizeof(float) * inputSize * inputSize * inputChannels);
+    EXPECT_EQ(cudaStatus, cudaSuccess);
 
-    conv2d.host_conv(input.data(), output.data());
+    cudaStatus = cudaMalloc((void**)&d_output, sizeof(float) * outputSize * outputSize * numFilters);
+    EXPECT_EQ(cudaStatus, cudaSuccess);
+
+    // // Copy input to device
+    cudaStatus = cudaMemcpy(
+        d_input, input.data(), sizeof(float) * input.size(), cudaMemcpyHostToDevice
+    );
+    EXPECT_EQ(cudaStatus, cudaSuccess);
+
+    conv2d.forward(d_input, d_output);
 
     std::vector<float> expected = {
         44.0f, 54.0f, 64.0f,
         84.0f, 94.0f, 104.0f,
         124.0f, 134.0f, 144.0f
-    };  
+    };
+    std::vector<float> output(outputSize * outputSize * numFilters);
+
+    cudaStatus = cudaMemcpy(
+        output.data(), d_output, sizeof(float) * output.size(),
+        cudaMemcpyDeviceToHost
+    );
+    EXPECT_EQ(cudaStatus, cudaSuccess);
 
     for (int i = 0; i < output.size(); ++i) {
         EXPECT_FLOAT_EQ(expected[i], output[i]);