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

Question

问题
我正在尝试使用卷积定理对两个测量数据 A 和 B 进行反卷积。我知道对于卷积，您应该对数据进行零填充以防止循环卷积。但是，我很困惑零填充对于反卷积是否也至关重要。

问题
1.如何正确地根据卷积定理进行反卷积？
2. 为什么下面的例子不起作用？

方法
由于测量了 A 和 B，因此我创建了一个示例以供进一步研究。这个想法是通过在same模式下使用scipy.signal.convolve来创建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]) 

Answer 1

事实上，由于内核是以 2.0(模糊和膨胀)为中心的 [1.0 2.0 1.0]，因此内核宽度为 3。由于数组 A 在 [0..5 ]，完整卷积数组 paddedB 在 [-1..6] 上不为空。尽管如此，函数scipy.signal.convolve(...,'same')返回一个主干卷积数组B(0..5)=paddedB(0..5)。因此，与 patedB(-1) 和 patedB(6) 相关的信息将会丢失，并且恢复内核会变得困难如果np.convolve() 的选项相同使用。

为了避免信息丢失，输出 patedB 将被填充以包含 support卷积信号的值，计算为Minkowski sum功能A的支持和内核的支持。 选项完整 np.convolve()直接计算 patedB，不会丢失信息。

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

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

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

如果无法获取pusedB并且B是唯一可用的数据，您可以尝试通过用零填充B来重建pusedB，或对 B 的最后值进行平滑。它需要对内核的大小进行一些估计。

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

最后，一个window可以应用于 patedA 和 poodleB，这意味着中间的值在要估计内核时更重要。例如 Parzen/de la Vallée Poussin window :

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]