Compute inverse IRF from FRF impact testing

Question

I have a measured FRF by impact testing and I have computed the IRF by ifft(). Now I want to compute the inverse IRF. I have read many posts in this forum about inverting FIR filters. As far as I understood, I must check where are my poles and zeros to see if my inverted FIR filter will be causal and stable. My initial idea was just inverting my FRF and computing again the corresponding IRF (inverse) by ifft(). Then, my poles become zeros and vice versa. When I do this I see that the inverse looks non-causal and quite noisy. Therefore, I want to check the zeros and poles. Here comes my first question:

1. How do you compute the zplot in Matlab from an FIR filter coming from a measurement?

My second question:

2. In this thread (https://www.dsprelated.com/showthread/comp.dsp/24411-1.php) I read that the inverted FIR filter might be longer (answer from Patrick Gaydecki). Why?

My third question, related to the second question:

3. Somewhere else I read that my inverted FIR might become an IIR. Is that true? And Why? Maybe that's why my inverted FIR becomes "longer".

I know there are other methods to compute the inverse FIR, (e.g. LMS inverse), but I want to try this method first.

Thank you in advance for your answers :)

Answer 1

When I do this I see that the inverse looks non-causal and quite noisy.

That's very common for a measured frequency responses unless you have very good signal to noise ratio at ALL frequencies including DC and Nyquist, which is rarely the case. If there are areas where $X(f)$ is small, the inverse is very large and that typically amplifies the noise at these frequencies a lot and that dominates the impulse response

How do you compute the zplot in Matlab from an FIR filter

An FIR filter has all poles at $z=0$. The zeros are simply the root of the filter coefficients. Let's look at a simple example, an FIR filter with the taps h = [1 .75 .5 .25];

The zeros are the simply roots of the coefficient polynomial, i.e.

>> h = [1 0.75 .5 0.25];
filZeros = roots(h)
filZeros =
     -0.60583 +          0i
    -0.072085 +    0.63833i
    -0.072085 -    0.63833i

The transfer function is

$$H(z) = \frac{h_0+h_1z^{-1}+h_2z^{-2}+h_3z^{-3}}{1} $$

And the inverse is simply

$$G(z) = H^{-1}(z) = \frac{1}{h_0+h_1z^{-1}+h_2z^{-2}+h_3z^{-3}} $$

This is an IIR filter with only poles (i.e. all the zeros at $z = 0$) a so-called allpole filter. So yes, the inverse of an FIR is an all-pole IIR and will always be "longer".

This may or may not be stable. If the zeros of the original FIR are all inside the unit circle, than the poles of the inverse will be too and hence the inverse is stable. In this case both filters are minimum phase.

How do you compute the zplot in Matlab from an FIR filter coming from a measurement?

Here is where things get tricky. Unfortunately, calculating the roots of a polynomial is numerically difficult and prone to large errors especially if the order is large and/or there is some noise (as there is on all real world measurement).

There are couple of techniques to try

1. If the system you are measuring isn't minimum phase, add some bulk delay to the inverse to make it causal
2. Try to "repair" the bad SNR spots in the frequency response by interpolation, extrapolation and or suitable smoothing
3. Try least square error solutions in the time domain.
4. Use general purpose IIR fitters in the frequency domain to get an approximation of the inverse with a limited amount of poles.