Noções básicas sobre deconvolve scipy
estou tentando entenderscipy.signal.deconvolve
.
Do ponto de vista matemático, uma convolução é apenas a multiplicação no espaço fourier, então eu esperaria que para duas funçõesf
eg
:Deconvolve(Convolve(f,g) , g) == f
Em numpy / scipy, esse não é o caso ou estou perdendo um ponto importante. Embora já existam algumas questões relacionadas à deconvolução no SO (comoaqui eaqui) eles não abordam esse ponto, outros permanecem obscuros (esta) ou sem resposta (aqui) Há também duas perguntas sobre o SignalProcessing SE (esta eesta) cujas respostas não são úteis para entender como a função deconvolve do scipy funciona.
A questão seria:
Como você reconstrói o sinal originalf
de um sinal complicado, supondo que você conheça a função convolutiva g.?Ou, em outras palavras: como esse pseudocódigoDeconvolve(Convolve(f,g) , g) == f
traduzir em numpy / scipy?Editar: Observe que esta pergunta não tem como objetivo evitar imprecisões numéricas (embora essa também seja umaquestão aberta), mas ao entender como convolve / deconvolve trabalham juntos em scipy.
O código a seguir tenta fazer isso com uma função Heaviside e um filtro gaussiano. Como pode ser visto na imagem, o resultado da desconvolução da convolução não é de modo algum a função Heaviside original. Ficaria feliz se alguém pudesse esclarecer esse assunto.
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: Note que existe umexemplo matlab, mostrando como convolver / desconvolver um sinal retangular usando
yc=conv(y,c,'full')./sum(c);
ydc=deconv(yc,c).*sum(c);
No espírito desta pergunta, também ajudaria se alguém pudesse traduzir este exemplo em python.