了解 scipy 反卷积

Question

我正在努力理解 scipy.signal.deconvolve。

从数学的角度来看，卷积只是傅里叶的乘法 space 所以我希望对于两个函数 f 和 g:
Deconvolve(Convolve(f,g) , g) == f

在 numpy/scipy 中，情况并非如此，或者我遗漏了重要的一点。尽管已经有一些与 SO 上的反卷积相关的问题（如 here and here) they do not address this point, others remain unclear (this) or unanswered (here). There are also two questions on SignalProcessing SE (this and this），但这些问题的答案对理解 scipy 的反卷积函数的工作原理没有帮助。

问题是：

如何从卷积信号重建原始信号 f，假设你知道卷积函数 g.?
或者换句话说：这个伪代码 Deconvolve(Convolve(f,g) , g) == f 是如何转化为 numpy / scipy?

编辑：请注意，这个问题的目的不是为了防止数值不准确（尽管这也是一个 open question），而是为了理解如何 convolve/deconvolve 在 scipy.

一起工作

以下代码尝试使用 Heaviside 函数和高斯滤波器来实现。从图中可以看出，卷积反卷积的结果并不在所有原始的 Heaviside 函数。如果有人能阐明这个问题，我会很高兴。

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，展示了如何使用convolve/deconvolve矩形信号

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

本着这个问题的精神，如果有人能够将这个例子翻译成 python。

也会有所帮助

Answer 1

经过反复试验，我发现了如何解释 scipy.signal.deconvolve() 的结果，并且我 post 我的发现作为答案。

让我们从一个工作示例代码开始

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)

一些重要的发现是：

过滤器应该比信号短
过滤器应该比零大得多（这里 > 0.013 就足够了）
对卷积使用关键字参数 mode = 'same' 可确保它与信号具有相同的数组形状。
反卷积有 n = len(signal) - len(gauss) + 1 个点。所以为了让它也驻留在同一个原始数组形状上，我们需要在两边将它扩展s = (len(signal)-n)/2。

当然，仍然欢迎对这个问题进一步的研究、评论和建议。

Answer 2

正如评论中所写，我无法帮助您最初发布的示例。正如@Stelios 指出的那样，由于数值问题，反卷积可能无法解决。

但是，我可以重现您在编辑中发布的示例：

这是直接翻译自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()

了解 scipy 反卷积

Understanding scipy deconvolve

python

numpy

scipy

signals

deconvolution