How to use the FFT for a 1D deconvolution?

Question

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])

See [Wiener Deconvolution](https://en.wikipedia.org/wiki/Wiener_deconvolution). — meowgoesthedog, Oct 25 '19 at 12:28

score 2 · Accepted Answer · edited Mar 30 '20 at 16:45

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]

thanks for your detailed answer, but why do you use `fftshift` after `ifft`? Although we are not in the frequency domain, you did mentionen to "shift to get zero frequency in the middle". Furthermore, the lenght of `dk` is in your example larger as the one of `A`. Is it possible to equal them in size? — beckstev, Oct 28 '19 at 09:45
may I ask you a question about an answer in this post: https://stackoverflow.com/a/54875879/10364180? Unfortunetly, I am not allowed to comment on this post, and I do not know how to contact you elsewhere. — beckstev, Oct 28 '19 at 14:45
The outcome of ifft is something near `[2 1 0 0 ... 1]`, because the middle of the kernel is at x=0, and x=0 is at index 0 for np.fft. The function `np.fftshift()` is applied to retreive the kernel at the middle of the frame, Indeed, you're right, we're not in the frequency domain, so it's not " ashift to get the central frequency in the middle", but "a shift to get absciss x=0 in the middle of the frame". — francis, Oct 28 '19 at 17:44
It is not mandatory, but a non-centered kernel results in the convoluted signal being shifted to the left or right compared to the original signal. It is visible in the [integral defining the convolution](https://en.wikipedia.org/wiki/Convolution). Kernels designed to blur a signal are therefore centered on x=0. — francis, Oct 29 '19 at 17:44

score 0 · Answer 2 · answered Mar 30 '20 at 04:33

# 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]

Instead of just pasting the answer into stack overflow, you should attempt to answer the question of how you came to that particular answer and explain each section etc. not just have a wall of unreadable code — Mark Davies, Mar 30 '20 at 08:03
HI, Mark - Thanks for the suggestion - I have to say that I struggled to get even this posting accepted - the interface kept telling me that my syntax was illegal, so I was not successful in providing the kind of detail that you rightly suggested. The code I provided was an edited version of the code provided by the previous respondent. I attempted to use that code but could not successfully run it in Python 3 without making the changes I highlighted by my introduction and comments. I thought others might find it useful — Dr. Saturn, Mar 31 '20 at 14:12

How to use the FFT for a 1D deconvolution?

2 Answers2