如何使用FFT进行一维反卷积？

Question

如何使用FFT进行一维反卷积？

3

问题
我想使用卷积定理对两个测量数据A和B进行反卷积。我知道在进行卷积时，应该对数据进行零填充以防止循环卷积。但是，我不确定在反卷积中是否也需要进行零填充。问题
1.如何正确地根据卷积定理执行反卷积？
2.为什么以下示例无法正常工作？方法
因为A和B是已测量的，所以我创建了一个示例来进一步研究。我的想法是使用scipy.signal.convolve在模式same下创建B。

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import convolve
from scipy.fftpack import next_fast_len

# A, in the description above
A = np.array([1, 1, 1, 2, 1, 1])
# The result, I want to get from the deconvolution
kernel = np.array([0, 1, 2, 1, 0, 0]) 
#B, in the description above
B = convolve(kernel, data, mode='same') 

# Using the deconvolution theorem
f_A = np.fft.fft(A)
f_B = np.fft.fft(B)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)

dk 的结果为：

dk = array([ 2.28571429-9.25185854e-18j,  1.28571429+9.25185854e-18j,
       -0.71428571-9.25185854e-18j, -0.71428571+9.25185854e-18j,
        0.28571429-9.25185854e-18j,  1.28571429+9.25185854e-18j])

期望结果是：

dk = array([0, 1, 2, 1, 0, 0])

- beckstev

data和con是什么？你是不是指的是A和B？ - undefined

@meowgoesthedog，是的！对于转移错误我很抱歉。 - undefined

2

请参阅维纳去卷积。 - undefined

2个回答

0

# I had to modify the listed code for it to work under Python3. 
# I needed to upgrade to the scipy-1.4.1 and numpy-1.18.2
# and to avoid a TypeError: slice indices must be integers
# I needed to change / to // in the line marked below
import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import convolve 
from scipy.fftpack import next_fast_len 
# A, in the description above 
A = np.array([1, 1, 1, 2, 1, 1]) 
kernel=np.asarray([1,2,1]) 
paddedB = convolve(kernel, A, mode='full') 
print(paddedB)
paddedA=np.zeros(paddedB.shape[0]) 
# note // instead of / below
paddedA[kernel.shape[0]//2: kernel.shape[0]//2+A.shape[0]]=A[:] 
#pad both signal and kernel. Requires the size of the kernel 
# Using the deconvolution theorem 
f_A = np.fft.fft(paddedA) 
f_B = np.fft.fft(paddedB) # I know that you should use a regularization here 
r = f_B / f_A 
# dk should be equal to kernel 
dk = np.fft.ifft(r) 
# shift to get zero abscissa in the middle: 
dk=np.fft.fftshift(dk) 
print(dk)
# this gives:
#(py36) bash-3.2$ python decon.py
#[1 3 4 5 6 5 3 1]
#[ 1.11022302e-16+0.j -1.11022302e-16+0.j -9.62291355e-17+0.j
# 1.00000000e+00+0.j  2.00000000e+00+0.j  1.00000000e+00+0.j
# 9.62291355e-17+0.j -1.11022302e-16+0.j]

- Dr. Saturn

不要只是把答案粘贴到Stack Overflow上，你应该尝试回答如何得出这个特定的答案，并解释每个部分等等。不要只有一堆难以阅读的代码。 - undefined

嗨，Mark - 感谢你的建议 - 我必须说我很难甚至让这个帖子被接受 - 界面一直告诉我我的语法是非法的，所以我没有成功提供你正确建议的那种详细信息。我提供的代码是前一个回答者提供的代码的编辑版本。我尝试使用那段代码，但在Python 3中无法成功运行，除非我进行了我在介绍和注释中强调的更改。我认为其他人可能会觉得它有用。 - undefined

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- francis · Accepted Answer

实际上，由于核心是[1.0 2.0 1.0]，以2.0为中心（模糊和膨胀），因此核宽度为3。由于数组A在[0..5]上非空，完整的卷积数组paddedB在[-1..6]上非空。然而，函数scipy.signal.convolve(...,'same')返回一个截断的卷积数组B(0..5)=paddedB(0..5)。因此，与paddedB(-1)和paddedB(6)相关的信息丢失了，并且如果使用np.convolve()的same选项，则恢复核变得困难。

为了避免信息损失，输出的 paddedB 要进行填充，以包含卷积信号的支持度。支持度计算为函数 A 的支持度和核函数的支持度的 Minkowski sum。np.convolve() 的选项 full 直接计算 paddedB，不会有信息损失。

kernel=[1,2,1]
paddedB = convolve(kernel, A, mode='full')

使用卷积定理检索内核，需要将输入信号 A 填充以匹配函数 paddedB 的支持。

paddedA=np.zeros(paddedB.shape[0])
paddedA[kernel.shape[0]/2: kernel.shape[0]/2+A.shape[0]]=A[:]

# Using the deconvolution theorem
f_A = np.fft.fft(paddedA)
f_B = np.fft.fft(paddedB)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)
# shift to get zero frequency in the middle:
dk=np.fft.fftshift(dk)

请注意使用函数 np.fft.fftshift() 来将零频率移动到中间。

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import convolve
from scipy.fftpack import next_fast_len

# A, in the description above
A = np.array([1, 1, 1, 2, 1, 1])

kernel=np.asarray([1,2,1])
paddedB = convolve(kernel, A, mode='full')
print paddedB

paddedA=np.zeros(paddedB.shape[0])
paddedA[kernel.shape[0]/2: kernel.shape[0]/2+A.shape[0]]=A[:]
#pad both signal and kernel. Requires the size of the kernel

# Using the deconvolution theorem
f_A = np.fft.fft(paddedA)
f_B = np.fft.fft(paddedB)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)
# shift to get zero abscissa in the middle:
dk=np.fft.fftshift(dk)

print dk

如果无法获得paddedB，而B是唯一可用的数据，则可以尝试通过用零填充B或平滑B的最后值来重建paddedB。这需要对核大小进行一些估计。

B = convolve(A,kernel, mode='same')
paddedB=np.zeros(A.shape[0]+kernel.shape[0]-1)
paddedB[kernel.shape[0]/2: kernel.shape[0]/2+B.shape[0]]=B[:]
print paddedB

最后，窗口可以应用于paddedA和paddedB，这意味着中间的值在估计核时更为重要。例如，Parzen / de la Vallée Poussin窗口：

import numpy as np 
import matplotlib.pyplot as plt 
from scipy.signal import convolve
from scipy.fftpack import next_fast_len
from scipy.signal import tukey
from scipy.signal import parzen

# A, in the description above
A = np.array([1, 1, 1, 2, 1, 1])

kernel=np.asarray([1,2,1])
paddedB = convolve(kernel, A, mode='full')
print paddedB


B = convolve(A,kernel, mode='same')
estimatedkernelsize=3
paddedB=np.zeros(A.shape[0]+estimatedkernelsize-1)
paddedB[estimatedkernelsize/2: estimatedkernelsize/2+B.shape[0]]=B[:]
print paddedB

paddedA=np.zeros(paddedB.shape[0])
paddedA[estimatedkernelsize/2: estimatedkernelsize/2+A.shape[0]]=A[:]

#applying window
#window=tukey(paddedB.shape[0],alpha=0.1,sym=True) #if longer signals, should be enough.
window=parzen(paddedB.shape[0],sym=True)
windA=np.multiply(paddedA,window)
windB=np.multiply(paddedB,window)


# Using the deconvolution theorem
f_A = np.fft.fft(windA)
f_B = np.fft.fft(windB)
# I know that you should use a regularization here 
r = f_B / f_A

# dk should be equal to kernel
dk = np.fft.ifft(r)
# shift to get the zero abscissa in the middle:
dk=np.fft.fftshift(dk)

print dk

然而，估计的核函数远非完美，因为A的大小较小：

[ 0.08341737-6.93889390e-17j -0.2077029 +0.00000000e+00j
 -0.17500324+0.00000000e+00j  1.18941919-2.77555756e-17j
  2.40994395+6.93889390e-17j  0.66720653+0.00000000e+00j
 -0.15972098+0.00000000e+00j  0.02460791+2.77555756e-17j]