了解科学解卷积 [英] Understanding scipy deconvolve

问题描述

我正在尝试了解 scipy.signal.deconvolve .

I'm trying to understand scipy.signal.deconvolve.

从数学的角度来看，卷积只是傅立叶空间中的乘法，所以我希望 两个功能f和g的功能:
Deconvolve(Convolve(f,g) , g) == f

From the mathematical point of view a convolution is just the multiplication in fourier space so I would expect that for two functions f and g:
Deconvolve(Convolve(f,g) , g) == f

在numpy/scipy中不是这种情况，或者我错过了重要的一点. 尽管已经有一些与去卷积相关的问题(例如，此处和此处)，它们没有解决这一点，其他方面仍不清楚(). SignalProcessing SE上还有两个问题( this 和

In numpy/scipy this is either not the case or I'm missing an important point. Although there are some questions related to deconvolve on SO already (like 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) the answers to which are not helpful in understanding how scipy's deconvolve function works.

问题可能是:

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

• How do you reconstruct the original signal f from a convoluted signal, assuming you know the convolving function g.?
• Or in other words: How does this pseudocode Deconvolve(Convolve(f,g) , g) == f translate into numpy / scipy?

编辑:请注意，此问题并非旨在防止数值不正确(尽管这也是

Edit: Note that this question is not targeted at preventing numerical inaccuracies (although this is also an open question) but at understanding how convolve/deconvolve work together in scipy.

以下代码尝试使用Heaviside函数和高斯滤波器进行此操作. 从图像中可以看出，卷积反卷积的结果不是 所有原始的Heaviside功能.如果有人可以阐明这个问题，我将感到非常高兴.

The following code tries to do that with a Heaviside function and a gaussian filter. As can be seen in the image, the result of the deconvolution of the convolution is not at all the original Heaviside function. I would be glad if someone could shed some light into this issue.

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示例，展示了如何使用

Edit: Note that there is a matlab example, showing how to convolve/deconvolve a rectangular signal using

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

本着这个问题的精神，如果有人能够将此示例转换为python，这也将有所帮助.

推荐答案

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

After some trial and error I found out how to interprete the results of scipy.signal.deconvolve() and I post my findings as an answer.

让我们从一个有效的示例代码开始

Let's start with a working example code

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)

This code produces the following image, showing exactly what we want (Deconvolve(Convolve(signal,gauss) , gauss) == signal)

一些重要的发现是:

• 滤波器应比信号短
• 过滤器在任何地方都应该大于零(这里> 0.013足够好)
• 对卷积使用关键字参数mode = 'same'可确保其与信号位于相同的阵列形状上.
• 反卷积具有n = len(signal) - len(gauss) + 1个点. 因此，为了使它也驻留在相同的原始数组形状中，我们需要将其在两侧扩展为s = (len(signal)-n)/2.

• The filter should be shorter than the signal
• The filter should be much bigger than zero everywhere (here > 0.013 is good enough)
• Using the keyword argument mode = 'same' to the convolution ensures that it lives on the same array shape as the signal.
• The deconvolution has n = len(signal) - len(gauss) + 1 points. So in order to let it also reside on the same original array shape we need to expand it by s = (len(signal)-n)/2 on both sides.

当然，仍然欢迎对此问题进行进一步的发现，评论和建议.

Of course, further findings, comments and suggestion to this question are still welcome.

这篇关于了解科学解卷积的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！