Entender scipy deconvolve

estoy tratando de entenderscipy.signal.deconvolve.

Desde el punto de vista matemático, una convolución es solo la multiplicación en el espacio de Fourier, por lo que esperaría que para dos funcionesf yg:
Deconvolve(Convolve(f,g) , g) == f

En numpy / scipy, este no es el caso o me falta un punto importante. Aunque ya hay algunas preguntas relacionadas con la deconvolución en SO (comoaquí yaquí) no abordan este punto, otros siguen sin estar claros (esta) o sin respuesta (aquí) También hay dos preguntas sobre SignalProcessing SE (esta yesta) las respuestas a las cuales no son útiles para comprender cómo funciona la función deconvolución de scipy.

La pregunta sería:

¿Cómo reconstruyes la señal original?f de una señal enrevesada, suponiendo que conozca la función enrevesada g.O en otras palabras: ¿cómo funciona este pseudocódigoDeconvolve(Convolve(f,g) , g) == f traducir a numpy / scipy?

Editar: Tenga en cuenta que esta pregunta no está dirigida a prevenir imprecisiones numéricas (aunque esto también es unpregunta abierta) pero al comprender cómo convolve / deconvolve trabajan juntos en scipy.

El siguiente código intenta hacer eso con una función Heaviside y un filtro gaussiano. Como se puede ver en la imagen, el resultado de la deconvolución de la convolución no es en absoluto la función Heaviside original. Me alegraría si alguien pudiera arrojar algo de luz sobre este tema.

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt

# Define heaviside function
H = lambda x: 0.5 * (np.sign(x) + 1.)
#define gaussian
gauss = lambda x, sig: np.exp(-( x/float(sig))**2 )

X = np.linspace(-5, 30, num=3501)
X2 = np.linspace(-5,5, num=1001)

# convolute a heaviside with a gaussian
H_c = np.convolve( H(X),  gauss(X2, 1),  mode="same"  )
# deconvolute a the result
H_dc, er = scipy.signal.deconvolve(H_c, gauss(X2, 1) )


#### Plot #### 
fig , ax = plt.subplots(nrows=4, figsize=(6,7))
ax[0].plot( H(X),          color="#907700", label="Heaviside",    lw=3 ) 
ax[1].plot( gauss(X2, 1),  color="#907700", label="Gauss filter", lw=3 )
ax[2].plot( H_c/H_c.max(), color="#325cab", label="convoluted" ,  lw=3 ) 
ax[3].plot( H_dc,          color="#ab4232", label="deconvoluted", lw=3 ) 
for i in range(len(ax)):
    ax[i].set_xlim([0, len(X)])
    ax[i].set_ylim([-0.07, 1.2])
    ax[i].legend(loc=4)
plt.show()

Editar: Tenga en cuenta que hay unejemplo de matlab, que muestra cómo convolucionar / desconvolucionar una señal rectangular usando

yc=conv(y,c,'full')./sum(c); 
ydc=deconv(yc,c).*sum(c); 

En el espíritu de esta pregunta, también ayudaría si alguien pudiera traducir este ejemplo a Python.