了解 scipy 反卷积 [英] Understanding scipy deconvolve

问题描述

我正在尝试了解

编辑:注意有一个

一些重要的发现是:

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

当然，我们仍然欢迎对这个问题进一步的发现、评论和建议.

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

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

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.

The question would be:

• 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.

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

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

In the spirit of this question it would also help if someone was able to translate this example into python.

解决方案

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

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

Some important findings are:

• 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.

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