Понимание scipy deconvolve

Question

Я пытаюсь понять scipy.signal.deconvolve.

С математической точки зрения свертка просто умножение в Фурье пространстве, так что я бы ожидать , что для двух функций f и g:
Deconvolve(Convolve(f,g) , g) == f

В NumPy/SciPy это или не так, или Мне не хватает важного момента. Хотя есть еще вопросы, связанные с deconvolve на SO уже (например, here и here), они не затрагивают этот момент, другие остаются неясными (this) или без ответа (here). Также есть два вопроса по SignalProcessing SE (this и this), ответы на которые не помогают понять, как работает функция deconvolve от scipy.

вопрос будет:

• Как вы восстановить исходный сигнал f из свернутого сигнала, если вы знаете функцию свертки г.?
• Или другими словами: Как этот псевдокод Deconvolve(Convolve(f,g) , g) == f перевести на numpy/scipy?

Edit: Обратите внимание, что этот вопрос не направлен на предотвращении численных неточностей (хотя это тоже open question), но при понимании того, как скручивать/деконволюцию работать вместе в SciPy.

Следующий код пытается сделать это с помощью функции Heaviside и gaussian filter. Как видно из изображения, результат деконволюции свертки не равен всей исходной функции Хевисайда. Я был бы рад, если бы кто-то мог пролить свет на этот вопрос.

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

Редактировать: Обратите внимание, что есть matlab example, показывающий, как свертка/деконволюция прямоугольного сигнала с использованием

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

В духе этого вопроса также помогло бы, если бы кто-то смог перевести этот пример в python.

Answer 1

После некоторых проб и ошибок я выяснил, как интерпретировать результаты scipy.signal.deconvolve(), и я отправляю свои результаты в качестве ответа.

Начнем с рабочего кода примера

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

Этот код дает следующее изображение, показывая именно то, что мы хотим (Deconvolve(Convolve(signal,gauss) , gauss) == signal)

Некоторые важные выводы:

• Фильтр должен быть shor чем сигнал
• Фильтр должен быть намного больше нуля везде (здесь> 0,013 достаточно)
• Использование аргумента ключевого слова mode = 'same' для свертки гарантирует, что он живет в той же форме массива, что и сигнал.
• Deconvolution имеет n = len(signal) - len(gauss) + 1 баллов. Итак, чтобы он также находился на той же исходной форме массива, нам нужно развернуть его на s = (len(signal)-n)/2 с обеих сторон.

Конечно, дальнейшие выводы, комментарии и предложения по этому вопросу по-прежнему приветствуются.

Answer 2

Как написано в комментариях, я не могу помочь с примером, который вы опубликовали первоначально. Как отметил @Stelios, деконволюция может не сработать из-за числовых проблем.

Я могу, однако, воспроизвести пример вы размещены в вашем Edit:

То есть код, который является прямым переводом из исходного 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()