Python反卷积给出了意想不到的结果答案

【问题标题】：Python deconvolution giving unexpected resultPython反卷积给出了意想不到的结果
【发布时间】：2021-11-23 23:27:29
【问题描述】：

下面我绘制了信号（寿命衰减），我试图从脉冲响应函数（即除法器 (IRF)）中反卷积。所以我应该让衰减更尖锐一点。这是我查看的主题的示例，它提供了我需要的内容：

请不要为我的代码，我只使用分频器的峰值 (IRF)，而不是图像上显示的整个数组序列。

我正在使用以下代码来做到这一点：

IRF = IRF * (max(decay)/max(IRF))
# replace 0s to avoid error message 
IRF = np.where(IRF == 0, 0.1, IRF)    
decay = np.where(decay == 0, 0.1, decay)  
# take only the quotient part of the result 
deconv = scipy.signal.deconvolve(decay, IRF)[0]
# "padding" the deconvolved signal so it has the same size as the original signal 
s = int((len(decay)-(len(deconv)))/2)  ## difference on each side 
deconv_res = np.zeros(len(decay))   
end = int(len(decay)-s-1)  # final index
deconv_res[s:end] = deconv
deconv = deconv_res
# convolved normalized to decay height for plotting 
deconv_n = deconv * (max(decay)/max(deconv))

IRF是一个float64数组，信号是一个uint16数组。

我承认我对反卷积的数学不太熟悉，所以我盲目地尝试不同的东西，比如尝试不同的除法器函数，但没有任何结果接近预期。我得到的最后一个结果是这样的（参见原始信号的图以及它试图去卷积的信号..）

如果我不明白 scipy.deconvolve 中的某些内容，谁能给我一些想法，这种奇怪行为的原因可能是什么，或者甚至是一些可以帮助我的高级阅读材料？或者，如果您认为这个问题更多的是物理问题而不是编码问题，那么提出这样一个问题的更好地方？

【问题讨论】：

标签： python scipy signal-processing deconvolution

【解决方案1】：

问题在于 deconvolve 是一种多项式除法，它将输出信号分解为 $conv(h, x) + r$，如果您的信号有噪声，它可能会产生奇怪的结果。此外，如果脉冲响应中的第一个样本很小，它往往会产生指数增长的输出。

我会为这个问题做的是 FFT 的划分。

N = 2**(ceil(log2(len(IRF) + len(decay)))
filtered = ifft(fft(decay, N) / fft(IRF, N))

【讨论】：

仍然必须处理除以零的问题。例如，将脉冲与 sinc 卷积。在频域中，这是一个（周期性）复指数乘以一个矩形窗口。用 sinc 去卷积以获得脉冲相当于除以具有零的矩形窗口。为什么这是一个错误是有道理的，对于反转周期函数的窗口没有唯一的答案。
这适用于一般情况。在这种情况下，由于滤波器的响应很短，并且过渡很陡，因此频域中不应该有零点：) 让我们看看作者是怎么说的。
通过转换到频域并返回，它工作得非常好（考虑到嘈杂的数据），非常感谢。我选择了你的答案作为解决方案。如果您不介意，您能否向我解释为什么它不像这里使用高斯滤波器的示例那样，是因为信号的形状不清晰吗？ stackoverflow.com/questions/40615034/…
谢谢。差异的解释更多地与噪声有关，如果前导系数很小并且滤波器或数据中存在噪声，则噪声会激发不稳定模式。 dsp.stackexchange.com/a/78325/55891。反卷积不会优化最小二乘误差，它会将您的信号分解为滤波部分和常数部分。