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.