如何使用 FFT 进行一维反卷积？
How to use the FFT for a 1D deconvolution?
问题
我正在尝试使用卷积定理对两个测量数据 A 和 B 进行反卷积。我知道对于卷积，您应该对数据进行零填充以防止循环卷积。但是，如果零填充对于反卷积也是必不可少的，我感到困惑。
问题
1.How我是否正确地根据卷积定理进行了反卷积？
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])
确实，由于内核是 [1.0 2.0 1.0] 以 2.0 为中心（模糊和膨胀），内核宽度是 3。由于数组 A 在 [0..5] 上是非空的，完整的卷积数组 paddedB 在 [-1..6] 上是非空的。然而，函数 scipy.signal.convolve(...,'same') returns 一个主干卷积数组 B(0..5)=paddedB(0..5)。因此，与 paddedB(-1) 和 paddedB(6) 相关的信息丢失了 ，恢复内核变得困难 if option same of np.convolve() 用于.
为避免信息丢失，将填充输出 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 或平滑最后一个值来重建 paddedB 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 可以应用于 paddedA 和 paddedB，这意味着中间的值更重要，因为要估计内核。例如 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]
# 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]
问题
我正在尝试使用卷积定理对两个测量数据 A 和 B 进行反卷积。我知道对于卷积，您应该对数据进行零填充以防止循环卷积。但是，如果零填充对于反卷积也是必不可少的，我感到困惑。
问题
1.How我是否正确地根据卷积定理进行了反卷积？
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])
确实，由于内核是 [1.0 2.0 1.0] 以 2.0 为中心（模糊和膨胀），内核宽度是 3。由于数组 A 在 [0..5] 上是非空的，完整的卷积数组 paddedB 在 [-1..6] 上是非空的。然而，函数 scipy.signal.convolve(...,'same') returns 一个主干卷积数组 B(0..5)=paddedB(0..5)。因此，与 paddedB(-1) 和 paddedB(6) 相关的信息丢失了 ，恢复内核变得困难 if option same of np.convolve() 用于.
为避免信息丢失，将填充输出 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 或平滑最后一个值来重建 paddedB 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 可以应用于 paddedA 和 paddedB，这意味着中间的值更重要，因为要估计内核。例如 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]
# 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]