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

> 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

"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

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

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

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