Problem
I am trying to deconvolve two measured data A and B using the convolution theorem. I know that for a convolution, you should zero pad your data to prevent circular convolution. However, I am confused if zero padding is also essential for a deconvolution.
Question
1.How do I perform a deconvolution based on the convolution theorem correctly?
2. Why does the following example not work?
Approach
Because A and B are measured, I created an example for further investigations. The idea is to create B by using scipy.signal.convolve in mode 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)
The result for dk is:
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])
Expected is:
dk = array([0, 1, 2, 1, 0, 0])
Indeed, since the kernel is [1.0 2.0 1.0] centered on 2.0 (a blurring and inflating), the kernel width is 3. Since the array A is non null on [0..5], the full convolved array paddedB is non null on [-1..6]. Nevertheless, function scipy.signal.convolve(...,'same') returns a trunkated convolved array B(0..5)=paddedB(0..5). Hence, information related to paddedB(-1) and paddedB(6) are lost, and recovering the kernel is made difficult if option same of np.convolve() is used.
To avoid that loss of information, the output paddedB is to be padded to contain the support of the convolved signal, computed as the Minkowski sum of the support of function A and the support of the kernel. The option full of np.convolve() directly computes paddedB without loss of information.
kernel=[1,2,1]
paddedB = convolve(kernel, A, mode='full')
To retreive the kernel using the convolution theorem, the input signal A is to be padded to match the support of function 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)
Notice the use of fucntion np.fft.fftshift() to get the zero frequency in the middle.
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
If obtaining paddedB is not possible and B is the only available data, you may try to reconstruct paddedB by padding B with zeros, or smoothing of the last values of B. It requires some estimate for the size of the kernel.
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
Finally, a window may be applied to both paddedA and paddedB, meaning that values in the middle are more significant as the kernel is to be estimated. For instance a 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
Nevertheless, the estimated kernel is far from being perfect, as the size of A is small:
[ 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]