Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Problem with numerical convolution
import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import math
import datetime
from matplotlib.font_manager import FontProperties
def convolveoriginal(x, y):
'''
The original algorithm from http://www.physics.rutgers.edu/~masud/computing/WPark_recipes_in_python.html.
'''
P, Q, N = len(x), len(y), len(x) + len(y) - 1
z = []
for k in range(N):
t, lower, upper = 0, max(0, k - (Q - 1)), min(P - 1, k)
for i in range(lower, upper + 1):
t = t + x[i] * y[k - i]
z.append(t)
return np.array(z) #Modified to include conversion to numpy array
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.
'''
P = len(y1) #Determine the length of the signal
z = [] #Create a list of convolution values
for k in range(P):
t = 0
lower = max(0, k - (P - 1))
upper = min(P - 1, k)
for i in range(lower, upper):
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 = [math.exp(-t) for t in time] #Create an array of function values
#exp2 = [2 * math.exp(-2 * t) for t in time]
exp2 = [t ** 2 * math.exp(-t) for t in time]
#Calculate the analytical expression
#analytical = [2 * math.exp(-2 * t) * (-1 + math.exp(t)) for t in time]
analytical = [1. / 3 * math.exp(-t) * t ** 3 for t in time]
#Calculate the trapezoidal convolution
trapezoidal = convolve(exp1, exp2, dt)
#trapezoidal2 = convolve(exp1, sqexp, dt)
#Calculate the scipy convolution
sci = signal.convolve(exp1, exp2, mode = 'full')
#Slice the first half to obtain the causal convolution and multiply by dt
#to account for the step width
sci = sci[0:steps] * dt
sci = np.r_[0, sci[:steps - 1]]
ratio = [a / b for a, b in zip(sci, analytical)]
#Calculate the convolution using the original Riemann sum algorithm
riemann = convolveoriginal(exp1, exp2)
riemann = riemann[0:steps] * dt
#Plot
plt.plot(time, analytical, label = 'analytical')
plt.plot(time, trapezoidal, 'o', label = 'trapezoidal')
plt.plot(time, sci, '.', label = 'signal.convolve')
plt.plot(time, ratio, '.', label = 'ratio: signal.convolve/analytical')
fontP = FontProperties()
fontP.set_size('small')
plt.legend(prop = fontP)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment