Deconvolution of two FIRs

Question

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!

Answer 1

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.

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