4
$\begingroup$

Basic questions: What's the "correct" way to deconvolve two causal FIRs in the frequency domain (i.e. using the FFT), neither of which may be minimum phase but both may be considered to have stable exact inverses? By "correct", I mean, how does one avoid time-aliasing errors? If there is no "correct" way, what is the best way?

Some context: Given two causal FIRs, it is easy to avoid time-aliasing when convolving them in the frequency domain by first appropriately padding the IRs with trailing zeros (i.e. performing linear convolution via the FFT as discussed in standard DSP texts). However, it's not clear to me how to do something similar for deconvolution when the FIRs may not be minimum phase.

Deconvolving a minimum phase FIR out of a causal FIR (minimum phase or not) seems straightforward since we can compute the exact inverse of the minimum phase FIR (which we know will be causal and stable), and then linearly convolve this inverse with the other FIR.

However, since the exact inverse of a causal FIR that's not minimum phase is non-causal (assuming the inverse is BIBO stable - we force the inverse to have a finite length), performing linear convolution using the FFT after computing the inverse does not make sense since padding a non-causal inverse with trailing zeros doesn't make sense.

Related example with somewhat baffling results: In the figures below, I linearly convolve $x1$ with $x2i~(i = 1,2,3)$ to produce $yi$. Then, I deconvolve $x1$ out of each $yi$ to try and reproduce each $x2i$. To do this, I pad $x1$ with trailing zeros until the length matches that of $yi$, compute the exact inverse of the padded $x1$, and then circularly convolve the inverse with $yi$. The relevant difference between the plots across columns are the signals $x2i$ which have a different number of trailing zeros ($i$ denotes the number of trailing zeros). Can anyone provide some explanation as to why the results seem reasonable in columns 1 and 3 but not in column 2? Also, this might seem to be a reasonable way to deconvolve two signals and yet the results obtained don't always appear satisfactory (as this example shows). Any thoughts?

Note: In the plot titles, I have used $*$ to denote linear convolution and $/$ to denote deconvolution (all done in the frequency domain as described above).

Thanks!

Example

$\endgroup$
1
  • $\begingroup$ I found this question interesting but I don't know the answer. ( I'm a DSP permanent beginner. ). If you find or found out the answer from some other venue, would you mind posting it. Thanks. $\endgroup$
    – mark leeds
    Apr 23, 2019 at 4:05

1 Answer 1

4
$\begingroup$

I'm not sure why the problem occurs for you with the x22 example, but I believe the problem is because the deconvolution is not being performed over the right length.

Deconvolution of an FIR (inverting an FIR) requires, in general, convolution using an infinite duration impulse response. So just using the default lengths in performing the deconvolution using an FFT might not generate the correct result.

I've implemented the same problem in R, but the last (additional) row of plots shows the results obtained if I do the deconvolution over a longer signal than the default. This generates the right results.

enter image description here

As you can see, my implementation generates bizarre results for the x21 example rather than the x22 example in yours.


R Code Below

x1 <- c(0,0,0,0,0,0,0,1,1,1,1,1,1,1,0,0)
x21 <- c(0,1,1,1,0)
x22 <- c(0,1,1,1,0,0)
x23 <- c(0,1,1,1,0,0,0)

my_conv <- function(x,y)
{
  out_length <- length(x)+length(y)-1
  x_pad <- c(x,rep(0,out_length-length(x)))
  y_pad <- c(y,rep(0,out_length-length(y)))

  z <- Re(fft(fft(x_pad) * fft(y_pad),inverse=TRUE))

  return(z/out_length)
}

my_deconv <- function(x,y, out_length=0)
{
  if (out_length == 0)
  {
    out_length <- length(x)+length(y)-1
  }
  x_pad <- c(x,rep(0,out_length-length(x)))
  y_pad <- c(y,rep(0,out_length-length(y)))

  z <- Re(fft(fft(x_pad) / fft(y_pad),inverse=TRUE))

  return(z/out_length)
}

par(mfrow=c(5,3))
plot(x1)
plot(x1)
plot(x1)
plot(x21)
plot(x22)
plot(x23)
y1 <- my_conv(x1,x21)
plot(y1)
y2 <- my_conv(x1,x22)
plot(y2)
y3 <- my_conv(x1,x23)
plot(y3)
plot(my_deconv(y1,x1))
plot(my_deconv(y2,x1))
plot(my_deconv(y3,x1))
plot(my_deconv(y1,x1,50))
plot(my_deconv(y2,x1,50))
plot(my_deconv(y3,x1,50))
$\endgroup$
5
  • 1
    $\begingroup$ thanks peter. that's interesting and educational. $\endgroup$
    – mark leeds
    Apr 23, 2019 at 23:19
  • $\begingroup$ @Peter K. Thank you. I understand that doing the deconvolution over a longer signal can remove the artifacts. While your approach (padding by an arbitrary 50 zeros) works for the example I provided, 50 zeros may not be enough for a different signal. This can be a problem if we don't know what the output should be. The core of my question was: Is there a standard way to do deconvolution given ANY two FIRs? If padding with zeros is the solution, can the number of zeros to be padded be known apriori without knowing what the outcome of the deconvolution should be? $\endgroup$
    – Rahul
    Apr 24, 2019 at 12:55
  • $\begingroup$ @Peter K. In other words, just like we know that padding ANY two signals x and y to length(x) + length(y) - 1 or more allows linear convolution with no aliasing, is there a similar "rule", or even just a rule of thumb, for deconvolution? $\endgroup$
    – Rahul
    Apr 24, 2019 at 12:56
  • 1
    $\begingroup$ @Rahul The problem is that deconvolution of an FIR means convolution by an IIR. That means, sometimes even 50 samples won't be enough. I'd suggest looking at the impulse response of the inverted FIR and choosing a length that captures 99% of the inverse impulse response's energy. But I'm not even sure that will work... $\endgroup$
    – Peter K.
    Apr 25, 2019 at 1:10
  • 1
    $\begingroup$ @Peter K. Thanks! I appreciate the time you've taken to answer my question. $\endgroup$
    – Rahul
    Apr 28, 2019 at 14:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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