secret
Created

discrete convolve vs. analytical

  • Download Gist
conv.py
Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
import numpy as np
import matplotlib.pyplot as plt
 
 
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.
'''
y1, y2 = map(np.asarray, [y1, y2])
P = len(y1) #Determine the length of the signal
z = np.zeros_like(y1) #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[k] = t
if dx is not None: #Is a step width specified?
z *= dx
return z
 
steps = 50 #Number of integration steps
maxtime = 5 #Maximum time
dt = float(maxtime) / steps #Obtain the width of a time step
 
t = np.linspace(0, maxtime, steps, endpoint=False) #Create an array of times
exp1 = np.exp(-t) #Create an array of function values
exp2 = 2*np.exp(-2*t)
 
# Calculate the analytical expression
analytical = exp2*(np.exp(t)-1)
# Calculate the trapezoidal convolution
trapezoidal = convolve(exp1, exp2, dt)
 
# Calculate the scipy convolution
sci = np.convolve(exp1, exp2)
# Slice the first half to obtain the causal convolution and multiply
# by dt to account for the step width
sci = sci[:steps]*dt
 
# shift right, multiply
sci = np.r_[0,sci[:steps-1]]*0.86
 
# Plot
plt.plot(t, analytical, label = 'analytical')
plt.plot(t, trapezoidal, 'o', label = 'trapezoidal')
plt.plot(t, sci, '.', label = 'numpy.convolve')
plt.legend()
plt.show()

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.