CODEWITHSUNDEEP
pythonnumpyconvolutiongaussian
NorrinRadd
NorrinRadd

Reputation: 565

Python: size of the resulting function of the convolution of two Gaussians with np.convolve

I am interested to optimize a function which is the convolution of two functions. The main problem is that my resulting function is completly of scale and i do not understand what np.convolve actually does.

I wrote a small script that should convolve two Gaussian, but the resulting Gaussian is much larger in size than the input functions:

from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import numpy as np


# https://stackoverflow.com/questions/18088918/combining-two-gaussians-into-another-guassian
def gauss(x, p): # p[0]==mean, p[1]==stdev, p[2]==heightg, p[3]==baseline                   
    a = p[2]
    mu = p[0]
    sig = p[1]
    #base = p[3]
    return a * np.exp(-1.0 * ((x - mu)**2.0) / (2.0 * sig**2.0)) #+ base

p0 = [0, 0.3, 1] # Inital guess is a normal distribution
p02 = [0, 0.2, 0.5]


xp = np.linspace(-4, 4, 2000)
convolved = np.convolve(gauss(xp, p0),gauss(xp, p02), mode="same")


fig = plt.figure()



plt.subplot(2, 1, 1)

plt.plot(xp, gauss(xp, p0), lw=3, alpha=2.5)
plt.plot(xp, gauss(xp, p02), lw=3, alpha=2.5)
plt.xlim([-2, 2])

plt.subplot(2, 1, 2)
plt.plot(xp, gauss(xp, p0), lw=3, alpha=2.5)
plt.plot(xp, gauss(xp, p02), lw=3, alpha=2.5)
plt.plot(xp, convolved, lw=3, alpha=2.5,label="too damn high?")
plt.legend()
plt.xlim([-2, 2])

plt.tight_layout()
plt.show()

The resulting gaussian after the convolution is much higher

Result

than what i expected (wikipedia):

convolution convolution1

Upvotes: 3

Views: 4197

Answers (1)

Crivella
Crivella

Reputation: 1007

You gotta renormalize for the dx between two x ticks.

Numpy is substituting an integration for a summation, but since the functions takes only the Y values it doesn't care about the volume element on the integration axis which you need to include manually.
I've had to deal with this problem as well and it is a pain when you start doing stuff with dx=1 and than all of a sudden you get a wrong result because of different x axis distribution.

xp = np.linspace(-4, 4, 2000)
dx = xp[1] - xp[0]
convolved = np.convolve(gauss(xp, p0),gauss(xp, p02), mode="same") * dx

!!NOTE: don't put the renormalization inside the function definition. dx should be counted only once because of the integral going into a summation. If you put it inside the function it will actually be counted twice because bot gaussian are generated using it.

enter image description here

PS: To try and understand this better you can try to generate the x axis data with different spacing and without renormalizing you will see that the height of your convolution will differ (the smaller the spacing the greater the height)

fig = plt.figure()
ax = fig.add_subplot(111)
for spacing in (100,500,1000,2000):
    spacing += 1
    xp = np.linspace(-4, 4, spacing)
    dx = xp[1] - xp[0]
    convolved = np.convolve(gauss(xp, p01),gauss(xp, p02), mode="same") * dx
    ax.plot(xp, convolved, lw=3, alpha=2.5,label="spacing = {:g}".format(8/spacing))

ax.set_title("Convolution with different x spacing. With renormalization")
fig.legend()
plt.show()

enter code here enter image description here

Upvotes: 8

Related Questions