Inverse FIR computation

Started by November 12, 2004
Hello all,

I am trying to compute the inverse of an FIR filter: something that
would reverse the effect of the first filter.

I just want to make sure my steps are correct:

- take the FFT of the FIR
- compute complex 1/f (in the frequency domain)
- invert back into the time domain. (IFFT)
- circular shift to adjust for the aliasing
- window the data

Does that seem right?
--Shafik

"Shafik" <shafik@u.arizona.edu> wrote in message
> Hello all, > > I am trying to compute the inverse of an FIR filter: something that > would reverse the effect of the first filter. > > I just want to make sure my steps are correct: > > - take the FFT of the FIR > - compute complex 1/f (in the frequency domain) > - invert back into the time domain. (IFFT) > - circular shift to adjust for the aliasing > - window the data > > Does that seem right? > --Shafik >
You can only do this if the zeros of your FIR filter lie within the unit circle.Otherwise you will get a divergent power series when you do the reciprocal. Eg if F(z^-1) = 1-2z^-1 and you compute 1/F it diverges. You can of course replect them within the unit circle first and then divide out. Why are you using FFTs - surely a simple recursion will perform the inverse. Tom
Shafik wrote:
> Hello all, > > I am trying to compute the inverse of an FIR filter: something that > would reverse the effect of the first filter. > > I just want to make sure my steps are correct: > > - take the FFT of the FIR > - compute complex 1/f (in the frequency domain)
Do you mean 1 over frequency, or 1 over frequency response? In the first case it won't work, in the second case you'll have trouble where the frequency response goes to zero.
> - invert back into the time domain. (IFFT) > - circular shift to adjust for the aliasing
Huh? Aliasing shouldn't be an issue if you're working with data that's already sampled.
> - window the data
Why?
> > Does that seem right? > --Shafik >
You'll never fully reverse the effects of a filter that has a frequency response that goes to zero or that drops below the noise level. The best you'll ever do is come close. Assuming I couldn't find a good result for this in a book, then as a first cut I would modify your step where you compute 1/H(w) with a threshold, and set my inverse function to be / 1/H(w) H(w) > threshold H_i(w) = | \ 0 H(w) < threshold This will give you an output signal that for some noise level is pretty close to the LMS estimate of your input signal. If you really, truly wanted to have an exactly optimal LMS estimate you'd need to know the signal strength and noise level, and you'd probably want to have a stationary input signal if you wanted the computations to finish before you keeled over from old age, and you'd probably end up with an H_i(w) that's pretty close to the one I gave you above, but with a nice swoopy curve that connected the 1/H(w) part with the 0 part. The above equation should be pretty good if you can set the threshold to more or less equal the SNR, and more or less applies even for nasty signals, nasty noise, nasty etc... -- Tim Wescott Wescott Design Services http://www.wescottdesign.com

Number 6 wrote:

> Why are you using FFTs - surely a simple recursion will perform the inverse.
Which one might that be? The one that finds the LMS inverse (arbitrary length, no time domain aliasing) in O(N^2) involves solving a Toeplitz matrix within the framework of the Levinson-Durbin algorithm but that's usually programmed iteratively rather than recursively. It's not simple, however, and I'm not sure there is recursive formulation. I've posted the Matlab code to do an FIR inversion or a division of two FIR's this way several times. I'll post it again here with a subject line that might help locate it with search engines: "FIR Inversion and Division With Levinson-Durbin, Matlab Code" Bob -- "Things should be described as simply as possible, but no simpler." A. Einstein
"Shafik" <shafik@u.arizona.edu> wrote in message news:<1100301484.477678.55210@c13g2000cwb.googlegroups.com>...
> Hello all, > > I am trying to compute the inverse of an FIR filter: something that > would reverse the effect of the first filter. > > I just want to make sure my steps are correct: > > - take the FFT of the FIR
Yes, you could do that, but you don't need to.
> - compute complex 1/f (in the frequency domain)
Eh... what does "f" mean here? Frequency? The spectrum coefficient at frequency f?
> - invert back into the time domain. (IFFT)
If you do the FFT thing, you would need an IFFT somewhere.
> - circular shift to adjust for the aliasing
Why would you expect aliasing to occur?
> - window the data
Why would you need this?
> Does that seem right?
No, it doesn't. What you attempt to do here, is to make a filter that has a filter domain transfer function that goes as H'(w) = 1/ H(w) where H(w) is the transfer function of the original filter and H'(w) is the transfer function of the inverse filter. There are several problems with this. First, if H(w) has zeros on the unit circle, H'(w) will have a pole on the unit circle, and you will have serious problems. Second, the filter spec of H'(w) is only partially specified. H'(w)== 1/H(w) for the Fourier frequencies w_n = 2 pi n/ N where n is integer and N is the FFT length, but in the intervals between the Fourier frequencies, both H(w) and H'(w) are unspecified, so you have no control over exactly "how inverse" H'(w) actually is. You do have some alternatives, though. For a minimum phase FIR, you can use a lattice representation for your FIR filter. In that case, you get the inverse filter for free. If your FIR filter actually is minimum phase, there exists a recursive method for generating the lattice coefficients. The same recursion detects whether the filter is minimum phase or not. If you find yourself with a mixed phase FIR filter, it is possible to use Wiener filters to design a filter that is "quasi inverse". My experience is that these filters introduce an amount of artifacts with the data, so one ought to be a bit careful when using them. Rune
> I am trying to compute the inverse of an FIR filter: something that > would reverse the effect of the first filter. > > I just want to make sure my steps are correct: > > - take the FFT of the FIR > - compute complex 1/f (in the frequency domain) > - invert back into the time domain. (IFFT) > - circular shift to adjust for the aliasing > - window the data
Shafik, Essentially, the method you describe is indeed correct, although the circular shifting is not necessary. As others have mentioned, you must beware of zeros in the frequency response, otherwise an ill-conditioned division will result when you perform the arithmetic. There are various schemes around this, and one of the most common (and simple) is to add an offset term to the denominator, obviating the risk of this instability. However, if the offset is too large the reconstruction effect of the filter becomes diluted. Also, if you are only interested in a linear phase filter, then the inversion is simply based on the magnitude of the response since there are only cosine terms in the expansion. However, if you are interested in both magnitude and phase, you need to perform a complex inversion. For your interest, the software package "Signal Wizard" will take any frequency response (as a text file) and generate the inverse FIR coefficients for you automatically. It works pretty well, and we regularly use it to restore audio recordings, either off-line (using wave files) or in real time (this is only possible if you also have the hardware module). Take a look at http://www2.umist.ac.uk/dias/pag/signalwizard.htm This website allows you to download a demo version of the software. Patrick Gaydecki

Patrick Gaydecki wrote:

> Essentially, the method you describe is indeed correct, although the > circular shifting is not necessary.
Circular shifting will be required if the FIR is mixed phase because the front of the inverse will be at the end of the result. You must also be aware that if the inverse is longer than the FIR itself then the result will be time aliased, i.e. the long inverse will be wrapped round and round and added into the result. This will generally be the case. This can be accomodated by adequated zero padding before doing the FFT and windowing afterward to smoothly truncate it to something shorter or by using a time domain approach that finds an LMS inverse of a specified size without aliasing. The latter is the Levinson-Durbin inverse. Bob -- "Things should be described as simply as possible, but no simpler." A. Einstein