我试图理解scipy.signal.deconvolve https://docs.scipy.org/doc/scipy-0.15.0/reference/generated/scipy.signal.deconvolve.html.
从数学的角度来看,卷积只是傅立叶空间中的乘法,所以我期望
对于两个函数f
and g
:
Deconvolve(Convolve(f,g) , g) == f
在 numpy/scipy 中,情况要么不是这样,要么我遗漏了重要的一点。
尽管已经存在一些与 SO 上的反卷积相关的问题(例如here https://stackoverflow.com/questions/17063775/convolution-and-deconvolution-in-python-using-scipy-signal and here https://stackoverflow.com/questions/15483346/scipy-deconvolution-function)他们没有解决这一点,其他人仍然不清楚(this https://stackoverflow.com/questions/21990327/convolution-deconvolution-using-scipy) 或未答复 (here https://stackoverflow.com/questions/36204207/what-are-the-constraints-on-the-divisor-argument-of-scipy-signal-deconvolve-to-e)。还有两个关于信号处理SE的问题(this https://dsp.stackexchange.com/questions/8287/deconvolution-in-python and this https://dsp.stackexchange.com/questions/32181/how-does-scipy-signal-deconvolve-work)的答案对于理解 scipy 的反卷积函数如何工作没有帮助。
问题是:
- 如何重建原始信号
f
从复杂的信号中,
假设你知道卷积函数 g.?
- 或者换句话说:这个伪代码是如何实现的
Deconvolve(Convolve(f,g) , g) == f
翻译成numpy/scipy?
Edit: 请注意,这个问题的目的并不是为了防止数字不准确(尽管这也是一个开放式问题 https://stackoverflow.com/questions/36204207/what-are-the-constraints-on-the-divisor-argument-of-scipy-signal-deconvolve-to-e)但要了解卷积/反卷积如何在 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()
Edit: 请注意,有一个matlab例子 https://terpconnect.umd.edu/~toh/spectrum/Deconvolution.html,展示如何使用以下方法对矩形信号进行卷积/解卷积
yc=conv(y,c,'full')./sum(c);
ydc=deconv(yc,c).*sum(c);
本着这个问题的精神,如果有人能够将这个例子翻译成Python,它也会有所帮助。