A quick implementation of a trapezoidal convolution in C.

performancemodule.c

#include <Python.h>
//Normally #include "arrayobject.h" should be sufficient. I need to fix my includes.
#include "/Library/Frameworks/Python.framework/Versions/7.2/lib/python2.7/site-packages/numpy/core/include/numpy/arrayobject.h"

/*
 * Convolution method.
 */
static PyObject* convolve(PyObject* self, PyObject* args)
{
	PyArrayObject *vec1, *vec2, *conv; //The python object proxies
	double *c1, *c2, *cconv; //Pointers to the data of the numpy arrays
	int length; //The length of the vectors (assumed equal)
	int k, i; //Iteration variables

	//Parse the argument tuple and extract two vectors
	if (!PyArg_ParseTuple(args, "O!O!", &PyArray_Type, &vec1, &PyArray_Type, &vec2) )
		return NULL;
	//Check whether the vectors are ok
	if (NULL == vec1 || NULL == vec2)
		return NULL;
	//Obtain pointers for reference
	c1=(double*)vec1->data;
	c2=(double*)vec2->data;
	//Create a numpy array as an output variable and obtain a pointer
	conv = (PyArrayObject *) PyArray_FromDims(1, vec1->dimensions, NPY_DOUBLE);
	cconv = (double*)conv->data;
	//Obtain the vector length
	length=vec1->dimensions[0];
	//Calculate the convolution
	for(k = 0; k < length; k++)
		for(i = 0; i < k; i++)
			cconv[k] += (c1[i] * c2[k - i] + c1[i + 1] * c2[k - (i + 1)]) / 2.0;

	//Return the new array
	return PyArray_Return(conv);
}

/*
 * Boiler plate code.
 */

static PyMethodDef performancemoduleMethods[] =
{
	 {"convolve", convolve, METH_VARARGS, "Convolve two arrays of equal length."},
	 {NULL, NULL, 0, NULL}
};

PyMODINIT_FUNC

initperformancemodule(void)
{
	 (void) Py_InitModule("performancemodule", performancemoduleMethods);
	 import_array(); //Required call for numpy to work
}

performancemodulesetup.py

from distutils.core import setup, Extension

module1 = Extension('performancemodule', sources = ['performancemodule.c'])

setup(name = 'Performance module',
        version = '1.0',
        description = 'This is a package to speed up computationally intensive tasks.',
        ext_modules = [module1])

performancetest.py

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import math
import datetime
import performancemodule

def convolve(y1, y2, dx = None):
    '''
    Compute the finite convolution of two signals of equal length.
    @param y1: First signal.
    @param y2: Second signal.
    @param dx: [optional] Integration step width.
    @note: Based on the algorithm at http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
    '''
    z = [] #Create a list of convolution values
    for k in range(len(y1)):
        t = 0
        for i in range(k):
            t += (y1[i] * y2[k - i] + y1[i + 1] * y2[k - (i + 1)]) / 2
        z.append(t)
    z = np.array(z) #Convert to a numpy array
    if dx != None: #Is a step width specified?
        z *= dx
    return z

steps = 50 #Number of integration steps
maxtime = 7 #Maximum time
dt = float(maxtime) / steps #Obtain the width of a time step
time = [dt * i for i in range (steps)] #Create an array of times
exp1 = np.array([math.exp(-t) for t in time]) #Create an array of function values
exp2 = np.array([t ** 2 * math.exp(-t) for t in time])
#Calculate the analytical expression
analytical = [1. / 3 * math.exp(-t) * t ** 3 for t in time]
#Calculate the trapezoidal convolution
trapezoidal = convolve(exp1, exp2, dt)
trapezoidalC = performancemodule.convolve(exp1, exp2) * dt

#Plot
plt.plot(time, analytical, label = 'analytical')
plt.plot(time, trapezoidal, 'o', label = 'trapezoidal')
plt.plot(time, trapezoidalC, '.', label = 'trapezoidal in C')

plt.legend()
plt.show()

#Primitive runtime measurement
runs = 10000
last = datetime.datetime.now()
for i in range(runs):
    signal.convolve(exp1, exp2, mode = 'full')
delta = datetime.datetime.now() - last
print "scipy.signal.convolve (seconds, microseconds)", delta.seconds, delta.microseconds

last = datetime.datetime.now()
for i in range(runs):
    performancemodule.convolve(exp1, exp2)
delta = datetime.datetime.now() - last
print "convolve in C (seconds, microseconds)", delta.seconds, delta.microseconds

last = datetime.datetime.now()
for i in range(runs):
    convolve(exp1, exp2)
delta = datetime.datetime.now() - last
print "convolve (seconds, microseconds)", delta.seconds, delta.microseconds