0

So I'm doing some research which involves analyzing the convolution between two functions. I've been using scipy's built in convolve function, but I'm getting results that differ from when I just manually calculate the convolution integral.

Specifically, I've found that increasing the resolution of the input vectors also increases the output of the convolve function by that same factor. In other words, I have to scale the result back down by a factor of dx.

For example, I wish to convolve these two functions, f(x) and g(x):

After both calculating manually and using scipy's built-in function, I get these curves:

Multiplying the scipy curve by a factor of dx makes the two lie on each other:

Why is this the case?

Below is the code used to make the plots:

import numpy as np
import scipy as sp
from matplotlib import pyplot as plt

def func(x):
    return a/(np.pi*w*(1+((x)/w)**2))

def gauss(x,mx,o,I):
    return I*np.exp(-(((x-mx)/(o))**2))

def weightSignal(x,mx,o,I):
    return gauss(x,mx,o,I)*func(x)

I = 1
o = 0.5

a = 1
w = 1

n = 1001

xi = -1
xf = 1
x = np.linspace(2*xi,2*xf,num=n)

dx = (2*xf-2*xi)/(n-1)

f = func(x)
g = gauss(x,0,o,I)

plt.figure(0)
plt.plot(x,f,label='f(x)')
plt.plot(x,g,label='g(x)')
plt.legend()

fg = sp.signal.convolve(f,g,'same')

fgm = np.zeros((len(x),1))
for i in range(len(fgm)):
    mx = x[i]
    avg = sp.integrate.quad(weightSignal,-np.inf,np.inf,args=(mx,o,I))
    fgm[i] = avg[0]

plt.figure(1)
plt.plot(x,fg,label='Convolution Function')
plt.plot(x,fgm,label='Manual')
plt.xlim([xi, xf])
plt.legend()

plt.figure(2)
plt.plot(x,fg*dx,label='Convolution Function')
plt.plot(x,fgm,label='Manual')
plt.xlim([xi, xf])
plt.ylim([0, 1])
plt.legend()
plt.show()
3
  • What are you asking and how does it relate to the title of your post (How to extract a square array from the center of a larger array?)?
    – Woodford
    Commented Dec 11, 2023 at 17:43
  • Sorry, title was an old question I was going to ask but ended up figuring out. Must have saved it on accident. It's now updated. What I'm asking is why do I need to multiply scipy's convolution function by the step size in order to get the right output? I asked this in the Math stack exchange, and it got deleted for not being related to math for some reason
    – 3edw
    Commented Dec 11, 2023 at 17:44
  • is this relevant - stackoverflow.com/questions/47373050/…? Something about doing a sum rather than integral.
    – hpaulj
    Commented Dec 11, 2023 at 17:52

1 Answer 1

3

The reason why fg = sp.signal.convolve(f, g, 'same') doesn't automatically correct for the scaling by dx is because of the difference between the continuous definition of convolution and the discrete implementation used in computational libraries like SciPy.

Continuous Convolution: The mathematical definition of convolution is a continuous operation that involves integrating the product of two functions over all possible values.

Discrete Convolution: Scipy - no actually you by defining a linespace - discretize(s) the function(s). This means we only have values at specific points.

In the case of sp.signal.convolve(f, g, 'same'), the 'same' mode means that the output will have the same size as the input. However, the result will be a sum of products at each point, and it doesn't automatically account for the spacing between these points.

Because the convolution in SciPy is a discrete operation, the spacing between data points (dx in your case) needs to be explicitly considered. That's why you need to multiply the result by dx to match the scaling in your manual calculation.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.