Created
January 17, 2012 10:56
-
-
Save tillahoffmann/1626219 to your computer and use it in GitHub Desktop.
Problem with numerical convolution
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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