IIR filter relaxation time computation

Question

Let's define the relaxation time of a filter as the delay needed for the impulse response of a filter to reach and remain under an amplitude threshold.

Is there a closed-form solution to the computation of an estimate or boundary of the relaxation time of a discrete recursive (IIR) filter from its poles/zeros or from the coefficients of its difference equation, except by computing a reasonable length of its response and checking the absolute value of the response?

Answer 1

You know that calculating the impulse response from the transfer function involves partial fraction decomposition, whether it's in $s$ or $z$. After that you simply apply the inverse $z$-transform to find out the response.

For example, let's say you have this transfer function:

$$H(z)=\dfrac{2+3z^{-1}+4z^{-2}}{7-6z^{-1}+5z^{-2}} \tag{1}\label{1}$$

The coefficients are bogus, but they're chosen so that the poles reside within the ROC. The partial fraction expansion will look like this:

$$H(z)=\dfrac{r_1}{1-p_1z^{-1}}+\dfrac{r_2}{1-p_2z^{-1}} \tag{2}\label{2}$$

From which the impulse response will be:

$$\begin{align} h[n]&=r_1p_1^n+r_2p_2^n \\ {}&=2|p|^n\big[\Re(r)\cos(n\angle{p})-\Im(r)\sin(n\angle{p})\big] \tag{3}\label{3} \\ {}&=2|p|^n|r|\sin(n\angle{p}+\angle{jr}) \tag{4}\label{4} \\ {}&=2|r|\exp(\log(p)n)\sin(n\angle{p}+\angle{jr}) \tag{5}\label{5} \end{align}$$

Notice that $\eqref{4}$ and $\eqref{5}$ use $jr$ as argument for the angle, that's because $\eqref{3}$ has $\pm\Im\sin\pm\Re\cos$, which is reversed compared to the usual angle rotation, $\Re\sin\pm\Im\cos$, which means the argument needs to be rotated (hence the $j$).

For numerical values, the poles and the residues can be found with Octave's residuez():

num = [2 3 4];
den = [7 -6 5];
[r, p] = residuez(num, den)
R =
  -0.2571 + 0.6136i
  -0.2571 - 0.6136i
P =
   0.4286 - 0.7284i
   0.4286 + 0.7284i

The angles for r and for 1j*r are different, but that can be circumvented by either using the first value in the vector for both r and p, or the second. With this you know that the impulse response will be:

% predefined constants
magR = abs(r(1));
angR = arg(1j*r(1));
logP = log(magP);
magP = abs(p(1));
argP = arg(p(1));
n = [0:20];
% the actual impulse response
h = 2*magR*exp(logP*n).*sin(n*angP+angR);
h(1) = num(1)/den(1);

The first value in the vector will be the direct value, given by the evaluation of $H(0)$ (no delay impulse). You can plot this side-by side with the result from impz() (for verification), the "continuous-time" version (for fine-grain viewing), and the decaying curve (to see the match):

t = linspace(0, 20, 1001);
g = 2*magR*exp(logP.*t).*sin(t*angP + angR);
c = 2*magR*exp(logP.*t);
plot(n,h,"or",n,impz(num,den)(1:21),"xb",t,g,"",t,c)

The discrete values overlap (the x with the o points), while the "analog" version runs through all of them, with the exponential correctly defining the decay.

---

The simplification above (with the indices for r and p) may work for this 2nd order example, but for higher orders you need a sum. I got too greedy with simplifications, so the following should address all orders:

num = [2 3 4];
den = [7 -6 5 -4 3];
[r, p] = residuez(num, den);
n = [0:25];
h = real(sum(  abs(r).*exp( n.*log( abs(p) ) ) ...
  .*sin( n.*arg(p) + arg(1j*r) )  ));
h(1) = num(1)/den(1);
t = linspace(0, 25, 1001);
g = sum( r.*p.^t );
c1 = sum( abs(r).*exp(t.*log(abs(p))) );
c2 = sum( abs(r).*abs(p).^t );
plot( n, h, "or", ...
      n, impz(num, den)(1:length(n)), "xb", ...
      t, g, "g", ...
      t, abs(g), ":k", ...
      t, c1, ".-c", ...
      t, c2, "-.m" )
grid on;

For the sake of computation you can predefine absR=abs(r), argP=arg(p), etc. This time I've modified the graphs to include:

• h -- the discrete impulse response as calculated above, red "o";
• impz() -- the discrete impulse response as given by the built-in function (as reference), blue "x";
• g -- the "continuous-time" approach (densely sampled), green line;
• abs(g) -- the absolute value if g for better comparing the decaying exponential envelope, black dotted line;
• c1 -- the calculated decaying envelope, with $\text{e}^{-x}$, cyan dot-dash line;
• c2 -- the "original" envelope, with $p_k^n$, magenta dash-dot line.

The same sum is applied to $\eqref{3}$, $\eqref{4}$, and $\eqref{5}$. The degree of $H(z)$ has been increased by adding elements to the denominator's vector (here a 4th order). And the picture shows an overlapping h with impz(), thus identical responses, a g that goes through all the points in h and impz(), and overlapping envelopes providing the boundary for abs(g):

To conclude on point, the main difference beteen the discrete and the analog exponential decays:

• the $s$-domain has $\text{e}^{\Re(p)t}$;
• the $z$-domain has $\text{e}^{\log(|p|)n}$.

If the $s$-domain is stable then its real part is negative, and if the $z$-domain is within the ROC then the absolute value will be <1 and, as such, the logarithm will be negative.