0
$\begingroup$

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:

  1. 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:

  1. 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 :)

$\endgroup$

1 Answer 1

1
$\begingroup$

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.
$\endgroup$
2
  • $\begingroup$ Dear Hilmar, thank you for your detailed answer :) From your comments, I got two new questions: 1. How could I find out if my filter is minimum phase? I know what a minimum phase filter is but I don't know how I could identify it looking at my time-domain data or at the FRF 2. When I compute the inverted IRF by inverting the FRF and using the ifft(), am I approximating already the IIR? How could I make a longer approximation of that IIR? I thought about "upsampling" the inverted FRF for that. I am very interested to see if there is any time-aliasing because of the IIR approximation. $\endgroup$
    – jfresnoz
    Commented Dec 26, 2021 at 22:36
  • $\begingroup$ 1) Log magnitude and phase of a minimum phase transfer function are Hilbert transforms of each other. You can Hilbert transform the log magnitude and compare to your measured phase. 2) Just make your FFT long enough. Any high-Q features or "close to zero" features in your original transfer function will be a problem. $\endgroup$
    – Hilmar
    Commented Dec 27, 2021 at 23:16

Your Answer

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

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