# Analytical way to generate an IIR that is the square root response of another IIR

For simplicity take the transfer function of a second order IIR in $z$-domain:

$$H(z) = \frac {b_0 + b_1 z^{-1} + b_2 z^{-2}} {a_0 + a_1 z^{-1} + a_2 z^{-2}}$$

Is it possible to generate a new transfer function that yields the square root of the magnitude response of $H(z)$? Assuming the new transfer function is minimum phase. I know this would be possible using various numerical approaches such as global optimization or least squares, but has this specific problem been tackled or are their any approaches which would make the problem attainable without resorting to a heavy weight numerical approach? An approach that worked in the $s$-domain for an equivalent $H(s)$ would be equally useful.

• looks like you want a "pinking filter". maybe you should just try to answer that, rather than the general $\sqrt{ |H(z)| }$ question. – robert bristow-johnson Apr 25 '16 at 20:42
• @RBJ I was hoping there may have been some research done using pure mathematics to generate fractional powers of filter magnitudes as it has practical uses, say, when employing tail cancelling IIRs which square the magnitude. My gut instinct tells me that that a fourth order filter would provide a damn good approximation to the square root of a second order filter (based on something similar - but not - De Moivre's Theorem for fractional powers - which can be used to solve for fractional order filters), where something like, say, Taylor expansion is too naive. – keith Apr 26 '16 at 9:17

The only time you'll be able to do this exactly is when you know that the poles and zeros of $H(z)$ are in complex conjugate pairs inside and outside of the unit circle.

Note that for systems with real-valued impulse responses, this means you need four of each: two inside and two outside the unit circle.

The keyword to search on is spectral factorization.

This means that if you can write $H(z)$ as: $$H(z) = G(z) G(z^{-1})$$ where $G(z)$ is a minimum phase system, then $G$ just takes all the poles and zeros of $H$ that are inside the unit circle (and half those that are on the unit circle).

There exist approximate algorithms for finding a solution, but I am not well-versed in those.

• Thanks Peter, I think the first part of your answer is closest to answering my original question (although treating it as a spectral factorisation is not going to help in my situation). I have found a way to increase the order of my prototype filter without changing it's slope which is equivalent to what I want. – keith Apr 26 '16 at 13:19
• OK. Feel free to ask a follow up question (separately). Thanks for the check mark! – Peter K. Apr 26 '16 at 14:05

Let me just point out in addition to Peter K.'s more definitive answer, that there's an intuitive way to see that you can't do this in the general case. That is, there is no finite order transfer function that gives the square root of a given arbitrary transfer function. Let me indeed first consider the corresponding case on the $s$ -plane. Take $b_0 = b_1 = b_2 = a_0 = a_2 = 0$, so you get $$H(s) = \frac{1}{s},$$ i.e. an integrator. The desired square-root magnitude of this is $1/\sqrt{|s|} = |s|^{-\frac{1}{2}}$. However, the asymptotic response of any finite order filter is $|s|^n$, where $n$ is an integer.

So there is a counter-example, and actually a similar thing holds for any odd-order filter. In the case of even order filters, Peter's answer gives the method for checking if the square root can be taken.

Since there is a one-to-one map between rational functions on the $s$-plane and the $z$ -plane (one such map is the BLT), the same holds on the $z$ -plane too.

For an approximate solution, one can try Taylor expanding. Let $$H(z) = \frac{\sum_k b_k z^k}{\sum_l a_l z^l}$$ (note the change in notation). Now, the desired response is $\sqrt{|H(z)|}$. Let us try, for the sake of simplicity, approximate $\sqrt{H(z)}$ directly $$\sqrt{H(s)} = \frac{\sqrt{\sum_k b_k z^k}}{\sqrt{\sum_l a_l z^l}} = \frac{\sqrt{p_1(z)}}{\sqrt{p_2(z)}},$$ where $p_1(z)$ and $p_2(z)$ are the nominator and denominator polynomials respectively. Now expand this around some $z_0$ by writing $p_i(z) = p_i(z_0) + (p_i(z) - p_i(z_0)) = p_i(z_0)(1 + \hat{p}_i(z))$, where $$\hat{p}_i(z) = \frac{p_i(z) - p_i(z_0)}{p_i(z_0)}$$ and $i = 1, 2$, so that $\sqrt{p_i(z)} = \sqrt{p_i(z_0)} \sqrt{1 + \hat{p}_i(z)}$. Using the Taylor series of the square root we get $$\sqrt{p_i(z)} = \sqrt{p_i(z_0)}\sum_{n = 0}^\infty \frac{(-1)^n (2n)!}{(1 - 2n)(n!)^2 4^n} \hat{p}_i(z)^n.$$ We truncate to finite order, and we have an approximation in some region of the $z$-plane. However, there is a problem: the radius of convergence of the series expansion is $|\hat{p}_i(z)| < 1$. When $\hat{p}_i(s) = -1$, $p_i(z) = 0$, so the radius of convergence is limited by zeroes of $p_i(z)$.

In other words, the poles and zeroes of the response $H(z)$ obstruct the convergence of the series. It is easy to convince yourself that as a matter of fact, if there is at least one pole or zero inside the unit circle, the series won't converge everywhere on the unit circle with any choice of $z_0$. This is further made worse by the fact that you will probably want to choose $z_0$ real, to conserve the symmetry between $z$ and $z^\ast$, i.e. keep the system real.

Of course, you would have been happy with just replicating the square-root magnitude response of the system, and we tried to do a bit more. However, I believe you'll hit essentially the same problem, as the roots of the polynomial are not changed by transforming the phase. I can't prove that explicitly right now, though.

The one exception is if you don't have any poles, you can reflect all of your zeroes outside the unit circle to make a maximum phase version of your response, and then find a convergent Taylor series in this exact case.

So the message is that a Taylor expansion will almost surely not give you anything useful.

Any of the standard fitting algorithms would work around this by not replicating the response well near the poles and zeroes of the original transfer function, which should be fine at least for those poles and zeroes not on the unit circle.

• Taylor expansion would work perfectly if I wanted a filter of infinite order, so there are solutions. Can anyone do better than infinite order? Accuracy to around 1e-3dB across the whole frequency spectrum using an analytical approach as opposed to a brute force algorithm would be OK. – keith Apr 26 '16 at 9:25
• I'd think an iterative algorithm would give better results, but I can write a few words about Taylor expansions if I find the time. – Timo Apr 26 '16 at 9:28
• @keith Updated the answer with points about the Taylor expansion method. – Timo Apr 26 '16 at 13:28
• Thanks @Timo, that's very interesting and useful, I've given you an up vote in the way of thanks :-) I am making some progress with some mathematics around the notion of H(s^2). – keith Apr 26 '16 at 13:34
• @keith thanks! I think the $H(s^2)$ will hit the same problems as I described here, because the underlying problem is that $\sqrt{}$ is not an analytic function so it's not easy to approximate with rational functions. However, I might be mistaken, as I don't have a proof. I'd be interested to hear if you do get it to work! – Timo Apr 26 '16 at 14:02