IIR filter relaxation time computation

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?

• If you have the transfer function then you can derive its time response (similar to this, for example), and the result will be a product or sum of exponentials and sines/cosines. The exponentials will give you what you're looking for. Jun 3, 2020 at 13:54
• I'm not sure I'm following your thought. I don't want to expand the impulse response. If I derive for example the z-domain expression for a 2nd order IIR filter, I obtain an expression with a 4th order denominator, i.e. I've doubled the poles. Jun 3, 2020 at 16:01
• If you know the 3 dB bandwidth you can use the rule of thumb that the 90%/10% fall time is approximately 0.35/BW in Hz. This is true for a first order system but reasonably accurate for any IIR system due to its dominant pole dominating the long term response. Jun 3, 2020 at 17:00
• So that said, you can identify the dominant pole and then use a simpler equation based on that. (Or if there isn't one dominant pole you can identify the dominant pole/zeros and use a simpler equation based on those) Jun 3, 2020 at 17:01
• @moala Well, that was a bit of a brainfart. I was thinking in terms of the inverse Z-transform, but what you get is not like what you get in s domain (what I expected, an $e^{-t}\cos t$-like response). You get a sum of coefficients with powers of the indices, something like $\sum_k A_k\,p_k^k\,u_k$, where $A$, $p$, and $u$ are defined, but $k$ is increasing. Maybe it can be deduced, mathematically, but I don't know how to (though, if it can be done, I wouldn't mind seeing it). Jun 3, 2020 at 18:41

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.

• I doubt OP is here, anymore, but I've updated the answer with a more generic approach (and a conclusion). Jun 30 at 23:13
• dunno if it's correct but +1 for effort & apparent quality Jul 1 at 15:55
• @OverLordGoldDragon Well, I say it is but, if it isn't, that's why I provided the code. Thank you for the vote of confidence, though. Jul 1 at 21:13
• Thanks, I'll have a look at some time in the future. Jul 17 at 22:41
• Oh, I just found it. I think this also answers my question Jul 19 at 9:37