2
$\begingroup$

Summary: From discrete convolution theorem, it is understandable that we need 2N-1 point DFT of both sequences in order to avoid circular convolution. If we need to do deconvolution of a given experimental M-point signal using the DFT, how do we avoid the circular effects. Do we still need to zero pad the signal to 2M-1 points for DFT deconvolution or not. Some sources say, without explanation, that zero padding is not needed in deconvolution.


It is a well known property of convolution, that convolution of two continuous functions is a multiplication in the frequency domain i.e., $y$(t)=$a$*$b$ becoming a multiplication after Fourier transform: F(y(t))= F($a$)F($b$). This was well known by early 1900s and clearly mentioned in 1941 as discussed in previous questions here.

In practical applications, deconvolution is also desirable. In deconvolution, two functions are divided in the Fourier domain to recover the original function, say $a$, if $y$(t) and $b$(t) are known. For example, if we wish to recover $a$, we can divide F($y$(t)) by F($b$) and do an inverse transform to get $a$. It may not be a mathematically rigorous way, but it is a popular technique in spectroscopy from an empirical perspective.

In majority of the instrumental applications, the function is a discrete and problem is to perform deconvolution using Discrete Fourier Transform. Now the point is that in the DFT, the multiplication of the DFT of two functions yield a circular convolution, not a linear one.

However, one can extract the true linear convolution by ensuring that the length of the DFT equal to 2N-1, where N is the length of the data points contained in the discrete function.

Illustration Suppose $x_{1}=[\underline{1}, 2,3]$ and $x_{2}=[\underline{1}, 1,1] .$ We can compute the linear convolution as $$ x_{3}[n]=x_{1}[n] * x_{2}[n]=**[\underline{1}, 3,6,5,3]** $$ If we instead compute $$ x_{3}[n]=\operatorname{IDFT}_{N}\left(\operatorname{DFT}_{N}\left(x_{1}[n]\right) \cdot \operatorname{DFT}_{N}\left(x_{2}[n]\right)\right) $$ we get $$ x_{3}[n]=\left\{\begin{array}{ll} {[\underline{6}, 6,6]} & N=3 \\ {[\underline{4}, 3,6,5]} & N=4 \\ **{[\underline{1}, 3,6,5,3]} & N=5** \\ {[\underline{1}, 3,6,5,3,0]} & N=6 \end{array}\right. $$

Linear convolution and DFT convolution match when the length of the DFT is at least 2N-1, which is 5 in this case.

Question: Do we need to apply the same length restriction of 2N-1 size of data when wish to divide in the frequency domain or not, in other words trying to do a linear deconvolution using DFT?

  1. Suppose this is a observed sampled signal with 256 points (left, Window 1). We also know what distorted it by linear convolution. Its length was 256 points.

What people do is that they do a 256 point DFT on the observed signal. DFT is also performed on the distorted function, 256 points.

a) Divide the DFT(observed signal)/DFT(distorting function);

enter image description here

b) Take the inverse DFT, we get

enter image description here

I searched plenty of places, books and papers but could not find a clear requirement about the length of DFT in the case of division/deconvolution. Thanks

$\endgroup$
9
  • $\begingroup$ I'd suggest you stand back and look again at the problem from which this issue arises -- convolution is a smoothing operation, suppressing high frequencies, so deconvoluton will obviously amplify high frequencies, so (numerical or other) noise. The deconvolution problem is ill-posed so is typically attacked with methods of regularisation rather than a direct inversion. This probably accounts for the paucity of references to it you noted. $\endgroup$ – J.J. Green Jan 15 at 9:18
  • $\begingroup$ If I understand correctly, your question is: we know $x_3 = x_1 \star x_2$ (of length $n+k-1$) and $x_1$ (of length $n$), how do we find $x_2$ (of length $k$) using DFT. The answer is simple: pad your $x_1$ and $x_3$ with zeroes so that both have length $N \geqslant n + k - 1$, calculate DFT, divide, calculate inverse DFT, and finally unpad (is that a word? Edit: is "trim" the right word?) to get $x_2$. $\endgroup$ – Mateusz Kwaśnicki Jan 15 at 9:49
  • $\begingroup$ @MateuszKwaśnicki, This is exactly the problem I am exploring. But let us come to an experimental problem. Suppose we have an experimental data with 2000 points and we also know what caused the distortion in the signal by linear convolution. In order to do DFT based deconvolution, should one do a DFT deconvolution with 2*2000-1 or just 2000 points. $\endgroup$ – M. Farooq Jan 15 at 14:35
  • $\begingroup$ If you know $x_3$ and $x_3$ has length $n+k-1 = 2000$, then you are good to go with $N = 2000$, why not. However, if you are working with real-life (inaccurate) data, then you should really stick to what J.J. Green wrote in his comment above. $\endgroup$ – Mateusz Kwaśnicki Jan 15 at 14:43
  • $\begingroup$ @MateuszKwaśnicki, I agree that DFT based deconvolution is not the perfect one but it is still very popular. Consider this as the experimental scenario. You will have an experimental data with N=2000, and a discrete "distorting function" which ruined the original signal by linear convolution, whose length is also N=2000. We wish to recover the true undistorted data. I think we should go for 2*2000-1 during DFT division in order to do recover it by linear deconvolution not circular deconvolution. $\endgroup$ – M. Farooq Jan 15 at 14:50
3
$\begingroup$

If you work in the frequency domain there is no need to zero-pad the data to achieve a minimum length. In the time domain the resulting signal will then be periodically extended beyond the boundaries.


Perhaps this could be the "reputable source" requested in the update: Direct frequency-domain deconvolution when the signals have no spectral inverse --- [open access]

Section 4 gives the algorithm: the time-domain signals $h(n)$ and $y(n)$ of length $N$ are padded by zeros to length $2N$, giving $h_1(n)$ and $y_1(n)$. The Fourier transforms $Y_1(k)$ and $H_1(k)$ have length $2N$, however, only the $N$ values at even $k$ are needed for the division/deconvolution.

So the algorithm requires a pair of $2N$ Fourier transforms but the inverse Fourier transform has length $N$.

$\endgroup$
10
  • $\begingroup$ Thanks, Let us say, we have a single Lorentzian peak as a signal with 1024 points. When we wish to decrease its width, one can take the DFT of the signal we can divide it with the DFT of another Lorentzian of smaller width than the signal. It also needs 1024 points. The key point is do we need to do zero pad both to 2x1024 points? And after inverse we keep the first 1024 points. In reality, I see no difference with or without zero padding. $\endgroup$ – M. Farooq Jan 17 at 22:13
  • 1
    $\begingroup$ an alternative to zero-padding in the time domain is interpolation in the frequency domain, so there is no need to zero pad; zero-padding adds no new information, it is just an efficient way to interpolate the intermediate frequencies you want to use in order to avoid the circular convolution; you might find this discussion instructive. $\endgroup$ – Carlo Beenakker Jan 17 at 22:27
  • 1
    $\begingroup$ if you wish one more reference: "Zero padding is only suitable for convolution and not for deconvolution.", see Deconvolution of time series in the laboratory --- zero padding interpolates the transform, and if you are not interested or do not need an interpolated signal in the time domain, there is no need to zero-pad the signal in the frequency domain before you transform back to the time domain; does this answer your question? $\endgroup$ – Carlo Beenakker Jan 19 at 11:47
  • 1
    $\begingroup$ that paper only does zero-padding in the time domain, not in the frequency domain; edge effects appear when you convolute in the time domain, when you deconvolute by division in the frequency domain there are no edge effects, so there is no need to zero-pad in the frequency domain. $\endgroup$ – Carlo Beenakker Jan 19 at 18:18
  • 1
    $\begingroup$ this is answered in the first paragraph of section 4 of the IEEE paper : the signal $y$ has length $N_1\equiv N$, the response function $h$ has length $N_2\equiv M\geq N$ and the time-domain padding should increase the length of both to $2N_1-1$, assuming that $N_1 is a power of two. $\endgroup$ – Carlo Beenakker Jan 20 at 10:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.