python - 如何使用 FFT 进行一维反卷积?
问题描述
问题
我正在尝试使用卷积定理对两个测量数据A和B进行反卷积。我知道对于卷积,您应该对数据进行零填充以防止循环卷积。但是,如果零填充对于反卷积也是必不可少的,我很困惑。
问题
1.如何正确执行基于卷积定理的反卷积?
2. 为什么下面的例子不起作用?
方法
因为A和B是测量的,所以我创建了一个示例以供进一步研究。这个想法是通过使用in mode创建B。scipy.signal.convolve
same
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')
返回一个主干卷积数组B(0..5)=paddedB(0..5)
。因此,与和相关的信息会丢失,如果paddedB(-1)
paddedB(6)
使用选项same
,np.convolve()
恢复内核将变得困难。
为了避免信息丢失,输出paddedB
将被填充以包含卷积信号的支持,计算为函数 A 的支持和内核的支持的Minkowski 和。直接计算而不丢失信息的选项full
。np.convolve()
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]
推荐阅读
- python - Jupyter notebook 抛出 SyntaxError: invalid syntax on correct input
- mysql - Mysql select查询其中列有两个值
- python - 从 VBA 调用 Python 脚本
- game-engine - GdScript 获取碰撞时的节点路径
- javascript - 将 int 推送到数组时自动将其转换为字符串
- elasticsearch - 格式化elasticsearch结果
- c# - 删除和查找数组中的元素
- ios - 我想在 getPicture - ionic ios 的警报中更改“您的使用信息”的文本
- npm - 如何为现有的 Lerna mono-repo 生成每个包的 package-lock.json
- react-native - ReactNative:顶部的应用程序名称