Понимание scipy-деконволюции ⇐ Python

[Anonymous]

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

С математической точки зрения свертка — это просто умножение в пространстве Фурье, поэтому я ожидаю, что
что для двух функций f и g< /код>:

Код: Выделить всё

Deconvolve(Convolve(f,g) , g) == f

В numpy/scipy это либо не так, либо я упускаю важный момент.
Хотя уже есть некоторые вопросы, связанные с деконволюцией на SO (например, здесь и здесь), они не затрагивают этот вопрос, другие остаются неясными (это) или без ответа (здесь). Есть также два вопроса о SignalProcessing SE (этот и этот), ответы на которые не помогают понять, как работает функция деконволюции scipy.

Вопрос будет такой:

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

Изменить: Обратите внимание, что этот вопрос направлен не на предотвращение числовых неточностей (хотя это также открытый вопрос), а на понимание того, как свертывать/обратить свертку работать вместе в scipy.

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

Код: Выделить всё

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, показывающий, как свернуть/обратить свертку прямоугольного сигнала с помощью

Код: Выделить всё

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

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

Подробнее здесь: https://stackoverflow.com/questions/406 ... deconvolve