Alternative derivation for a Thiran filter

Question

Following this reasoning (link on ee.se) I'm trying to derive the transfer function for a Thiran filter, but I get stuck. I know of the original paper, I just thought I'd go this way. This is what I'm doing:

• Construct the symbolic transfer function (2nd order, for simplicity) with unity numerator, which will result in an unscaled magnitude. This is not a problem, and it can be resolved by using the sum of the denominator coefficients.
• Derive the group delay and, according to the logic in the linked answer, equate all the same power coefficients to form a system of equations. This is where it fails with the roots coming out as either nonsense, or complex.

$$\begin{align} H(z)&=\dfrac{1}{1+az^{-1}+bz^{-2}} \tag{1} \\ \tau_{gd}(\omega)&=-\dfrac{2b\cos(2\pi f)+a(3b+1)\cos(\pi f)+a^2+2b^2}{2b\cos(2\pi f)+2a(b+1)\cos(\pi f)+a^2+b^2+1} \tag{2} \\ &\begin{cases} a(3b+1)&=2a(b+1) \\ a^2+2b^2&=a^2+b^2+1 \end{cases} \tag{3} \\ &\Rightarrow \\ &\begin{cases} a=x,\; b=1 \quad\text{("x" can be anything)}\\ a=0,\;b=-1 \end{cases} \tag{4} \end{align}$$

Both solutions in (4) are wrong, and it doesn't matter whether I'm using a constant in the numerator, or some $z^{-n}$ (which makes sense), they're still wrong. Then I thought I'd adapt (1) to include a proper numerator, with two, overlapped zeroes at -1:

$$\begin{align} H(z)&=\dfrac{1+2z^{-1}+z^{-2}}{1+az^{-1}+bz^{-2}} \tag{5} \\ \tau_{gd}(\omega)&=-\dfrac{a(b-1)\cos(\pi f)+b^2-1}{2b\cos(2\pi f)+2a(b+1)\cos(\pi f)+a^2+b^2+1} \tag{6} \\ &\begin{cases} -a(b-1)&=2a(b+1) \\ 1-b^2&=a^2+b^2+1 \end{cases} \tag{7} \\ &\Rightarrow \\ &\begin{cases} a=0,\; b=0 \\ a=\pm j\dfrac{\sqrt2}{3},\;b=-\dfrac13 \\ \end{cases} \tag{8} \end{align}$$

So I thought I'd evaluate both (2) and its derivative at 0 Hz, then equate the former to 1, and the latter to 0, to form the system of equations, but the derivative of (2) at 0 Hz is a blunt zero, no term (the numerator has $\sin(n\omega)$ terms, everywhere).

And yet another problem is that there are infinitely many derivatives for the group delay if the variable is $\omega$, since it's a $\cos(\omega)$, so should I derive w.r.t. the whole $\cos(\omega)$ (e.g. replace $\cos(\omega)=x$)? After all, I'm already forming the system of equations based on the coefficients of $\cos(\omega)$. So I did, for both (2) and (6), and that's where I got the first real solutions, but still wrong. For the first (both in raw and float):

[[a=0,b=-1/3],[a=-(3*sqrt(10)+10)/5,b=(2*sqrt(10)+5)/5],[a=(3*sqrt(10)-10)/5,b=-(2*sqrt(10)-5)/5],[a=-2,b=1]]
[[a=0.0,b=-0.3333333333333333],[a=-3.897366596101028,b=2.264911064067352],[a=-0.1026334038989724,b=-0.2649110640673518],[a=-2.0,b=1.0]]

and for the second:

[[a=0,b=0],[a=-(2*sqrt(7)+4)/3,b=(sqrt(7)+2)/3],[a=(2*sqrt(7)-4)/3,b=-(sqrt(7)-2)/3],[a=-2,b=1]]
[[a=0.0,b=0.0],[a=-3.097167540709727,b=1.548583770354863],[a=0.4305008740430605,b=-0.2152504370215302],[a=-2.0,b=1.0]]

Where the solutions involve zero, discard, where they're -2 and 1, again, discard, and the rest are real, but still wrong (I won't add the plots, but you can take my word for it, or test for yourselves).

But this works for the Bessel 2nd order, where:

$$\begin{align} B(s)&=\dfrac{b}{s^2+as+b} \tag{9} \\ \tau_{gd}^B(\omega)&=\dfrac{a\omega^2+ab}{\omega^4+(a^2-2b)\omega^2+b^2} \tag{10} \\ &\begin{cases} a&=a^2-2b \\ ab=b^2 \end{cases} \tag{11} \\ &\Rightarrow \\ &\begin{cases} a=0,\;b=0 \\ a=1,\;b=0 \\ a=3,\;b=3 \end{cases} \tag{12} \end{align}$$

But then, something funny happened when I accidentally typed 2 instead of 1 in diff(d(f),f,2), given (2). Without the 3rd argument it defaults to the 1st differentiation, but it was in plot and I was changing the 3rd arg. dynamically between 0 and 1, when the typo happened. And then I copy-pasted that to use with solve(), and this is what came up by solving:

[[b=1,a=-2],[b=1/5,a=-4/5]]

Which are exactly the coefficients for the Thiran filter, except (1) has 1 in the numerator while Thiran has 2/5 (gain normalization). Then I tried the same thing with (5), but the results came out either both zero, or -2 and 1, so then I remembered that the numerator adds its own delay, and so I equated (2) with 2, not 1, and the results came out exactly as above. So I took it a step further and used the 3rd order variant for (1) and, sure enough, these were the solutions:

[[c=%r14,b=1-2*%r14,a=%r14-2],[c=-2/3,b=7/3,a=-8/3],[c=-1/14,b=3/7,a=-1]]

where the last set has the Thiran coefficients.

So now I have to wonder: is the double differentiation for the group delay needed in the digital domain, to avoid the $\sin()$ problem? Doesn't this sort of break the (N-1) derivatives condition? What is the explanation?

---

The question remains, but I'd just like to add one very beautiful thing, in case it went unnoticed: the fact that you can add a numerator, other than a simple gain normalizing constant, makes this the equivalent of an inverse Chebyshev with flat group delay. $\scriptsize\text{Or maybe a sort of a Gaussian with zeroes}$. Anyway, here's how a fairly random transfer function looks like for a 4th order Thiran denominator (unity delay), combined with two 2nd order numerators that give zeroes to form a CloseEnough®™ stopband equiripple:

$$\dfrac{(1+0.5z^{-1}+z^{-2})(1+1.7z^{-1}+z^{-2})}{1-\dfrac87z^{-1}+\dfrac{9}{14}z^{-2}-\dfrac{4}{21}z^{-3}+\dfrac{1}{42}z^{-4}} \tag{13}$$

Blue is the magnitude in dB and red is 30x the negative group delay (to fit in with the magnitude). The only downside to this is the extra delay introduced by the numerator, which is 1 for each 2nd order, thus the overall group delay will be 2+Thiran.

Answer 1

I think I found the problem. I dropped it then but, now it irked me again and, as I looked over the past scribblings, I realized I compared the method to the analog version and I didn't follow it. Now it struck me: the group delay for the analog version is always double the order. If differentiated, the function will always be zero when evaluated at DC, due to the odd powers in the numerator. Therefore, it turns out that even the analog version can't be solved this way and it needs to have a 2nd order differentiation so that its evaluation at DC is non null. This is the reason for the digital filters, too, since they also have the same characteristics: the group delay is a function of $\cos(\omega)$, differentiated it results in a product of $\cos(\omega)\sin(\omega)$ which is null at DC. A second differentiation "cures" it.

I won't type all the formulas because they're a lot to type, feel free to do it yourselves but, for example, (2) in the OP, differentiated once, results in two terms of the form $(a+b\cos(x))(c\sin(x))+\sin(x)$, which is zero at DC. A second differentiation yields 4 terms involving at least one "clean" $a+b\cos(x)$. Evaluation at DC is no longer zero.

And now, it works just fine for both versions: form a system of equations with all derivatives at DC zero, while $\tau_{\mathrm{gd}}(0)=1$, then solve. For Thiran, reusing (2) in the OP:

$$\begin{align} \tau^{(2)}_{\mathrm{gd}}(0)&=-\dfrac{a+2b}{1+a+b} \tag{14} \\ \left.\dfrac{\mathrm{d}^2\tau^{(2)}_{\mathrm{gd}}(\omega)}{\mathrm{d}\omega}\right|_{\omega=0}&=-\dfrac{2(a+4b+ab)(a+2b)}{(1+a+b)^3}+\dfrac{a+8b+3ab}{(1+a+b)^2} \tag{15} \\ &\begin{cases} (14)=1 \\ (15)=0 \end{cases} \\ \Rightarrow\qquad &a=\left[-2,\, -\dfrac45\right]\; ,\quad b=\left[1,\;\dfrac15\right] \tag{16} \end{align}$$

Out of which the second values of (16) are the ones for Thiran. The visual proof:

For the sake of consistency, here is for the 3rd order:

$$\begin{align} H^{(3)}(z)&=\dfrac{1+a+b+c}{1+az^{-1}+bz^{-2}+cz^{-3}} \tag{17} \\ \tau^{(3)}_{\mathrm{gd}}(\omega)&=-\dfrac{3c\cos(3\omega)+(4ac+2b)\cos(2\omega)+(5bc+3ab+a)\cos(\omega)+a^2+2b^2+3c^2}{2c\cos(3\omega)+2(ac+b)\cos(2\omega)+2(bc+ab+a)\cos(\omega)+a^2+b^2+c^2+1} \tag{18} \\ \tau^{(3)}_{\mathrm{gd}}(0)&=-\dfrac{a+2b+3c}{1+a+b+c} \tag{19} \\ \left.\dfrac{\mathrm{d}^2\tau^{(3)}_{\mathrm{gd}}(\omega)}{\mathrm{d}\omega}\right|_{\omega=0}&=-\dfrac{\big(8(ac+b)+2(a+ab+bc+9c)\big)(a+2b+3c)}{(1+a+b+c)^3}+\dfrac{8(2ac+b)+a+3ab+5bc+27c}{(1+a+b+c)^2} \tag{20} \\ \left.\dfrac{\mathrm{d}^4\tau^{(3)}_{\mathrm{gd}}(\omega)}{\mathrm{d}\omega}\right|_{\omega=0}&=\mathrm{a\; too\; large\; to\; type\; formula...} \tag{21} \\ &\begin{cases} (19)=1 \\ (20)=0 \\ (21)=0 \end{cases} \\ \Rightarrow\qquad &a_3=-1\; ,\quad b_3=\dfrac37\; ,\quad c_3=-\dfrac{1}{14} \tag{22} \end{align}$$

I've only shown the needed solution out of three but, it works for all orders. Visual proof:

And, for the sake of completion, for 2nd order analog Bessel (fortunately a bit simpler):

$$\begin{align} H^\mathrm{a}(s)&=\dfrac{b}{s^2+as+b} \tag{23} \\ \tau^\mathrm{a}_{\mathrm{gd}}(\omega)&=\dfrac{a(\omega^2+b)}{\omega^4+(a^2-2b)\omega+b^2} \tag{24} \\ &\begin{cases} \tau^a_{\mathrm{gd}}(0)=\dfrac{a}{b}=1 \\ \left.\dfrac{\mathrm{d}^2\tau^\mathrm{a}_{\mathrm{gd}}(\omega)}{\mathrm{d}\omega}\right|_{\omega=0}=\dfrac{2a(3b-a^2)}{b^3}=0 \end{cases} \tag{25} \\ \Rightarrow\quad &a=[0,3]\; ,\quad b=[0,3] \tag{26} \end{align}$$

This, too, works for 3rd order and above.