Zrozumienie scipy deconvolve

Question

Próbuję zrozumieć scipy.signal.deconvolve.

Z matematycznego punktu widzenia splot jest tylko mnożenie w przestrzeni Fouriera więc spodziewałbym że dla dwóch funkcji f i g:
Deconvolve(Convolve(f,g) , g) == f

W numpy/scipy to jest albo nie przypadek lub Brakuje mi ważnej sprawy. Mimo że istnieją pewne pytania związane z dekonvolve na SO już (jak here i here), nie odnoszą się do tego punktu, inne pozostają niejasne (this) lub nieodebrane (here). Są również dwa pytania dotyczące SignalProcessing SE (this i this), na które odpowiedzi nie są pomocne w zrozumieniu działania funkcji dekonwolucji scipy.

pytanie byłoby:

• Jak odtworzyć oryginalny sygnał f z zawiłe sygnału zakładając wiesz funkcję splatanie g.?
• Innymi słowy: w jaki sposób ten pseudokod Deconvolve(Convolve(f,g) , g) == f przekształca się w numpy/scipy?

Edit: Należy zauważyć, że kwestia ta nie jest ukierunkowane na zapobieganie nieścisłości numeryczne (chociaż jest to także open question), ale na zrozumieniu, jak prace convolve/deconvolve razem w scipy.

Następujący kod próbuje to zrobić za pomocą funkcji Heaviside i filtru gaussowskiego. Jak widać na obrazku, wynikiem dekonwolucji splotu nie jest żadna oryginalna funkcja Heaviside. Byłbym szczęśliwy, gdyby ktoś mógł rzucić trochę światła na tę kwestię.

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() 

Edit: Zauważ, że nie jest matlab example, pokazujący jak convolve/deconvolve prostokątny sygnał używając

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

W duchu tym pytaniu również pomogłoby, gdyby ktoś mógł przetłumaczyć ten przykład na pytona.

Answer 1

Po kilku próbach i błędach dowiedziałem się, jak interpretować wyniki scipy.signal.deconvolve() i zamieszczam moje wyniki jako odpowiedź.

Zacznijmy działającego kodu przykład

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

# let the signal be box-like 
signal = np.repeat([0., 1., 0.], 100) 
# and use a gaussian filter 
# the filter should be shorter than the signal 
# the filter should be such that it's much bigger then zero everywhere 
gauss = np.exp(-((np.linspace(0,50)-25.)/float(12))**2) 
print gauss.min() # = 0.013 >> 0 

# calculate the convolution (np.convolve and scipy.signal.convolve identical) 
# the keywordargument mode="same" ensures that the convolution spans the same 
# shape as the input array. 
#filtered = scipy.signal.convolve(signal, gauss, mode='same') 
filtered = np.convolve(signal, gauss, mode='same') 

deconv, _ = scipy.signal.deconvolve(filtered, gauss) 
#the deconvolution has n = len(signal) - len(gauss) + 1 points 
n = len(signal)-len(gauss)+1 
# so we need to expand it by 
s = (len(signal)-n)/2 
#on both sides. 
deconv_res = np.zeros(len(signal)) 
deconv_res[s:len(signal)-s-1] = deconv 
deconv = deconv_res 
# now deconv contains the deconvolution 
# expanded to the original shape (filled with zeros) 


#### Plot #### 
fig , ax = plt.subplots(nrows=4, figsize=(6,7)) 

ax[0].plot(signal,   color="#907700", label="original",  lw=3) 
ax[1].plot(gauss,   color="#68934e", label="gauss filter", lw=3) 
# we need to divide by the sum of the filter window to get the convolution normalized to 1 
ax[2].plot(filtered/np.sum(gauss), color="#325cab", label="convoluted" , lw=3) 
ax[3].plot(deconv,   color="#ab4232", label="deconvoluted", lw=3) 

for i in range(len(ax)): 
    ax[i].set_xlim([0, len(signal)]) 
    ax[i].set_ylim([-0.07, 1.2]) 
    ax[i].legend(loc=1, fontsize=11) 
    if i != len(ax)-1 : 
     ax[i].set_xticklabels([]) 

plt.savefig(__file__ + ".png") 
plt.show()  

Ten kod daje następujący obraz, pokazujący dokładnie to, co chcemy (Deconvolve(Convolve(signal,gauss) , gauss) == signal)

Niektóre ważne ustalenia:

• Filtr powinien być shor ter niż sygnał
• Filtr powinien być wszędzie dużo większy niż zero (tutaj> 0,013 jest wystarczająco dobry)
• Użycie słowa kluczowego argument mode = 'same' do splotu zapewnia, że ​​żyje on w tym samym kształcie tablicy co sygnał.
• Dekonwolucja ma n = len(signal) - len(gauss) + 1 punktów. Aby więc mógł on również znajdować się na tym samym oryginalnym kształcie tablicy, musimy go rozwinąć po obu stronach, s = (len(signal)-n)/2.

Oczywiście, dalsze ustalenia, komentarze i sugestie na to pytanie są nadal mile widziane.

Answer 2

Jak napisano w komentarzach, nie mogę pomóc z przykładami, które pierwotnie wysłałeś. Jak zauważył @Stelios, dekonwolucja może nie zadziałać z powodu problemów numerycznych.

mogę jednak odtworzyć przykład pan pisał w swojej EDIT:

Jest to kod, który jest bezpośrednie tłumaczenie z kodu źródłowego Matlab:

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

x = np.arange(0., 20.01, 0.01) 
y = np.zeros(len(x)) 
y[900:1100] = 1. 
y += 0.01 * np.random.randn(len(y)) 
c = np.exp(-(np.arange(len(y)))/30.) 

yc = scipy.signal.convolve(y, c, mode='full')/c.sum() 
ydc, remainder = scipy.signal.deconvolve(yc, c) 
ydc *= c.sum() 

fig, ax = plt.subplots(nrows=2, ncols=2, figsize=(4, 4)) 
ax[0][0].plot(x, y, label="original y", lw=3) 
ax[0][1].plot(x, c, label="c", lw=3) 
ax[1][0].plot(x[0:2000], yc[0:2000], label="yc", lw=3) 
ax[1][1].plot(x, ydc, label="recovered y", lw=3) 

plt.show()