了解 scipy 反卷积

Question

我正在尝试理解

scipy.signal.deconvolve

。

从数学的角度来看，卷积只是傅里叶空间中的乘法，所以我期望 对于两个函数

f

和

g

:

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

在 numpy/scipy 中，情况要么不是这样，要么我遗漏了重要的一点。 尽管已经存在一些与 SO 上的反卷积相关的问题（例如here和here），但它们没有解决这一点，但其他问题仍然不清楚（this）或未得到解答（here）。还有两个关于 SignalProcessing SE 的问题（this 和 this），这些问题的答案对于理解 scipy 的反卷积函数如何工作没有帮助。

问题是：

• 如何从卷积信号中重建原始信号

  f

  ， 假设你知道卷积函数 g.?
• 或者换句话说：这个伪代码

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

  如何翻译成numpy / scipy？

编辑：请注意，这个问题的目的不是为了防止数字不准确（尽管这也是一个开放问题），而是为了理解卷积/反卷积如何在 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示例，展示了如何使用

对矩形信号进行卷积/解卷积

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

)

一些重要的发现是：

• 滤波器应比信号短
• 过滤器在任何地方都应该比零大得多（这里 > 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()

Answer 3

如果您有二进制模式的脉冲响应系统：“开启”模式或“关闭”模式 - 在过滤器中仅定义这 2 个模式就足够了。如果你想让它们变大（正如你在OP中提到的）简单做：

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter

original = np.repeat([0., 1., 0., 0., 0., 1.], 100)
impulse_response_filterKernel = x =  [20, 1]
x= np.arange(len(original))

filtered = signal.convolve(impulse_response_filterKernel, original)
recovered, remainder = signal.deconvolve(filtered, impulse_response_filterKernel)
print(recovered)

fig, ax = plt.subplots(1,1)
ax.plot(x, original, label="original", lw=7, alpha=0.2)
ax.plot(x, filtered[:len(filtered)-1], label="filtered", lw=3)
ax.plot(x, recovered, label="recovered", color='black', lw=1, )
plt.legend()
plt.show()

如果您需要高斯滤波器平滑所有噪声点，可以参见here手动进行纠正