0

I'm looking to convolve a Gaussian with an asymmetric Gaussian. My attempt at this is below.

import numpy as np
import matplotlib.pyplot as plt

x=np.linspace(0,5,500)
dx=x[1]-x[0]

# Gaussian
sigma1=0.1
peak1=1
gauss1=np.exp(-(x-peak1)**2/(2*sigma1**2))

# Asymmetric Gaussian
gauss2=0.5*np.exp(-(x-1.5)**2/(-0.2+x*0.4)**2)

# convolution
conv=np.convolve(gauss1,gauss2,mode='same')*dx

plt.plot(x,gauss1,label='Gaussian')
plt.plot(x,gauss2,label='Asymmetric Gaussian')
plt.plot(x,conv,label='Convolution')
plt.xlim(0,5)
plt.legend()
plt.show()

Convolution plot

I don't understand why the resultant curve is positioned where it is. I would have thought it would have a peak at some x location between that of the two original curves?

freja
  • 47
  • 1
  • 6

1 Answers1

1

The problem is that in the convolution function, you used the parameter mode='same'. This mode only returns the middle values of the convolution as described in the official documentation. To get the actual convolution output, you need to use mode='full' to get the entire convolution result. In your case, given two inputs of 500 values, the output will have a size of 500 + 500 - 1 = 999.

Here's the edited code.

import numpy as np
import matplotlib.pyplot as plt

x=np.linspace(0,5,500)
dx=x[1]-x[0]

# Gaussian
sigma1=0.1
peak1=1
gauss1=np.exp(-(x-peak1)**2/(2*sigma1**2))

# Asymmetric Gaussian
gauss2=0.5*np.exp(-(x-1.5)**2/(-0.2+x*0.4)**2)

# convolution
x1 = np.linspace(0,9.99,999)
conv=np.convolve(gauss1,gauss2,mode='full')*dx

plt.plot(x,gauss1,label='Gaussian')
plt.plot(x,gauss2,label='Asymmetric Gaussian')
plt.plot(x1, conv,label='Convolution')
plt.legend()
plt.show()

The resultant peak need not lie in between that of the original curves. Say your peak was at 1 and 1.5. The resultant peak will be at 1 + 1.5 = 2.5.

enter image description here

NITHIN C BABU
  • 101
  • 4
  • Cheers, that helps a lot. – freja Dec 17 '22 at 21:02
  • Whilst the x1 array needs to have N*M+1 elements when using mode='full', should it also run to 2*max(x)? If I set the x1 array to have 999 elements in this example, but set it have the same limits as x (0 to 5), the convolved curve peak is of course in a different location. – freja Dec 18 '22 at 13:35
  • Theoretically, if you have two signals of size 5, the output will be of size 10. As we are dealing with discrete data, it is N+M-1. Squeezing the output to the range [0, 5] is not a standard procedure. It does not represent the actual convolution. You can shift the inputs to control the starting position of the output. – NITHIN C BABU Dec 21 '22 at 10:22