# Mathematical question that comes out of using bilinear transform

So this is related to the Cookbook and I tried solving it maybe two decades ago, gave up, and was reminded of the unsolved problem. But it's pretty damn straight forward, but I still got slogged down in the muck.

This is a simple Bandpass filter (BPF) with resonant frequency $\Omega_0$ and resonance $Q$:

$$H(s) = \frac{\frac{1}{Q}\frac{s}{\Omega_0}}{\left(\frac{s}{\Omega_0}\right)^2 + \frac{1}{Q}\frac{s}{\Omega_0} + 1}$$

At the resonant frequency

$$|H(j\Omega)| \le H(j\Omega_0) = 1$$

and the upper and lower bandedges are defined so that

$$|H(j\Omega_U)|^2 = \left|H\left(j\Omega_0 2^{BW/2} \right)\right|^2 = \tfrac12$$

$$|H(j\Omega_L)|^2 = \left|H\left(j\Omega_0 2^{-BW/2} \right)\right|^2 = \tfrac12$$

We call these the "half-power bandedges". Because we're audio, we define bandwidth in octaves, and in the analog world this bandwidth in octaves, $BW$, is related to $Q$ as:

$$\frac{1}{Q} = \frac{2^{BW} - 1}{\sqrt{2^{BW}}} = 2 \sinh\left( \tfrac{\ln(2)}{2} BW \right)$$

We're using bilinear transform (with pre-warped resonant frequency) which maps:

\begin{align} \frac{s}{\Omega_0} & \leftarrow \frac{1}{\tan(\omega_0/2)} \, \frac{1-z^{-1}}{1+z^{-1}} \\ \\ \frac{j \Omega}{\Omega_0} & \leftarrow \frac{j\tan(\omega/2)}{\tan(\omega_0/2)} \\ \end{align}

letting $z=e^{j \omega}$ and $s=j\Omega$.

The resonant angular frequency of the analog filter is $\Omega_0$, and with frequency warping compensation done to the resonant frequency in the realized digital filter, when $\omega=\omega_0$ (the user-defined resonant frequency), then $\Omega=\Omega_0$.

So if analog angular frequency is

$$\frac{\Omega}{\Omega_0} = \frac{\tan(\omega/2)}{\tan(\omega_0/2)}$$

then it's mapped to the digital angular frequency as

$$\omega = 2 \arctan\left( \frac{\Omega}{\Omega_0} \, \tan(\omega_0/2) \right)$$

Now, the upper and lower bandedges in the analog world are

$$\Omega_U = \Omega_0 2^{BW/2}$$ $$\Omega_L = \Omega_0 2^{-BW/2}$$

and in digital frequency domain are

\begin{align} \omega_U &= 2 \arctan\left( \frac{\Omega_U}{\Omega_0} \, \tan(\omega_0/2) \right) \\ &= 2 \arctan\left(2^{BW/2} \, \tan(\omega_0/2) \right) \\ \end{align}

\begin{align} \omega_L &= 2 \arctan\left(\frac{\Omega_L}{\Omega_0} \, \tan(\omega_0/2) \right) \\ &= 2 \arctan\left(2^{-BW/2} \, \tan(\omega_0/2) \right) \\ \end{align}

Then the actual difference, in log frequency of the bandeges (which is the actual bandwidth in the digital filter) is:

\begin{align} bw & = \log_2(\omega_U) - \log_2(\omega_L) \\ & = \log_2\Bigg(2\arctan\left( 2^{ BW/2} \, \tan(\omega_0/2) \right) \Bigg) - \log_2\Bigg(2\arctan\left( 2^{-BW/2} \, \tan(\omega_0/2) \right) \Bigg) \ \end{align}

or

$$\ln(2) \, bw = \ln\left(\arctan\left( e^{\ln(2)BW/2} \, \tan(\omega_0/2) \right) \right) - \ln\left(\arctan\left( e^{-\ln(2)BW/2} \, \tan(\omega_0/2) \right) \right)$$

This has a functional form of

$$f(x) = \ln\left(\arctan\left(\alpha \, e^{x} \right) \right) - \ln\left(\arctan\left( \alpha \, e^{-x}\right) \right)$$

where $f(x) \triangleq \ln(2) \, bw$, $x \triangleq \frac{\ln(2)}{2}BW$ and $\alpha \triangleq \tan(\omega_0/2)$

What I want to do is invert $f(x)$ (but I know I can't do it exactly with a nice closed form). I have already done a first-order approximation and i want to bump it up to a third-order approximation. And this has become a sorta copulating female canine, even though it should be straight forward.

Now this has something to do with the Lagrange Inversion Formula and I only want to take it to one more term than I have.

We know from above that $f(x)$ is an odd-symmetry function:

$$f(-x) = -f(x)$$

This means $f(0)=0$ and all even-order terms of the Maclaurin series will be zero:

$$y = f(x) = a_1 x + a_3 x^3 + ...$$

The inverse function is also odd symmetry, goes through zero, and can be expressed as a Maclaurin series

$$x = g(y) = b_1 y + b_3 y^3 + ...$$

and if we know what $a_1$ and $a_3$ are of $f(x)$, then we have a good idea what $b_1$ and $b_3$ must be:

$$b_1=\frac{1}{a_1} \quad \quad \quad \quad b_3 = -\frac{a_3}{a_1^4}$$

Now, I am able to calculate the derivative of $f(x)$ and evaluate it at zero and i get

$$\\ a_1 = \frac{2 \alpha}{(1+\alpha^2) \arctan(\alpha)} = \frac{\sin(\omega_0)}{\omega_0/2}$$ $$\\ b_1 = \frac{(1+\alpha^2) \arctan(\alpha)}{2 \alpha} = \frac{\omega_0/2}{\sin(\omega_0)} \\$$

But I am having a bitch of a time getting $a_3$ and therefore $b_3$. Can someone do this? I'd even settle for a solid expression for the third derivative of $f(x)$ evaluated at $x=0$.

• Just to clarify: Your goal is to invert $$f(x) = \ln\left(\arctan\left(\alpha \, e^{x} \right) \right) - \ln\left(\arctan\left( \alpha \, e^{-x}\right) \right)$$, i.e. for a given $f(x)$, you want to find the $x$? In particular, you want to do it by polynomial expansion and you are looking for the 3rd coefficient (since the 2nd is zero do to odd-ity of the function). Right? – Maximilian Matthé Feb 18 '17 at 7:08
• So you want to know $BW$ given $bw$, i.e., you want to know what bandwidth of the analog filter you need to choose to obtain a desired bandwidth of the digital filter, right? – Matt L. Feb 18 '17 at 7:17
• yes, yes, and yes. – robert bristow-johnson Feb 18 '17 at 19:08
• @robertbristow-johnson I didn't read the question too carefully, but I did notice you are interested in $f'''(x)$ at $x=0$. Is it ok to use Mathematica or Wolfram Alpha to compute that? I get a pretty clean result: $\frac{4 (8-\pi^2) \alpha^3}{\pi^3}$. wolframalpha.com/input/… And if you remove the "evaluate at x=0" part, Wolfram spits out the couplating female canine in its full glory. – Atul Ingle Feb 19 '17 at 2:24
• Typo in my $f(x)$ there. The "clean" result is actually: $-(6 a^2)/((a^2 + 1)^2 atan(a)^2) + (2 a)/((a^2 + 1) atan(a)) + (16 a^5)/((a^2 + 1)^3 atan(a)) + (12 a^4)/((a^2 + 1)^3 atan(a)^2) - (16 a^3)/((a^2 + 1)^2 atan(a)) + (4 a^3)/((a^2 + 1)^3 atan(-1)(a)^3)$ wolframalpha.com/input/… – Atul Ingle Feb 19 '17 at 2:31

To complement my part to this question: Here is a somewhat shorted answer based upon a manual expansion of the odd function $f(x)$ \begin{align*} f(x)&=\ln\left(\arctan\left(\alpha e^x\right)\right)-\ln\left(\arctan\left(\alpha e^{-x}\right)\right)\\ &=f_1x+f_3x^3+O\left(x^5\right)\tag{1} \end{align*} into a series up to the third order. Some more details can be found at mathSE.

At first we focus at the left-hand term $$\ln\left(\arctan\left(\alpha e^x\right)\right)$$ of $f(x)$ and start with

Series expansion of $\arctan$:

We obtain \begin{align*} \arctan(\alpha e^x)&=\sum_{n=0}^\infty \frac{(-1)^n}{2n+1}\alpha^{2n+1} e^{(2n+1)x}\\ &=\cdots\\ &=\sum_{j=0}^\infty\frac{1}{j!}\sum_{n=0}^\infty(-1)^n(2n+1)^{j-1}\alpha^{2n+1}x^j\tag{2} \end{align*}

We now derive from (2) the coefficients up to $x^3$. Using the coefficient operator $[x^k]$ to denote the coefficient of $x^k$ in a series we obtain \begin{align*} [x^0]\arctan\left(\alpha e^x\right)&=\sum_{n=0}^\infty \frac{(-1)^n}{2n+1}\alpha^{2n+1}=\arctan \alpha\\ [x^1]\arctan\left(\alpha e^x\right)&=\sum_{n=0}^\infty (-1)^n\alpha^{2n+1}=\frac{\alpha}{1+\alpha ^2}\\ [x^2]\arctan\left(\alpha e^x\right)&=\frac{1}{2}\sum_{n=0}^\infty (-1)^n(2n+1)\alpha^{2n+1}\\ &=\cdots =\frac{\alpha}{2}\cdot\frac{d}{d\alpha}\left(\frac{\alpha}{1+\alpha^2}\right) =\frac{\alpha(1-\alpha^2)}{2\left(1+\alpha^2\right)^2}\\ [x^3]\arctan\left(\alpha e^x\right)&=\frac{1}{6}\sum_{n=0}^\infty (-1)^n(2n+1)^2\alpha^{2n+1}\\ &=\frac{\alpha^2}{6}\sum_{n=0}^\infty (-1)^n(2n+1)(2n)\alpha^{2n-1} +\frac{\alpha}{6}\sum_{n=0}^\infty (-1)^n(2n+1)\alpha^{2n}\\ &=\cdots =\left(\frac{\alpha^2}{6}\cdot\frac{d^2}{d\alpha^2}+\frac{\alpha}{6}\cdot\frac{d}{d\alpha}\right)\left(\frac{\alpha}{1+\alpha^2}\right)\\ &=\cdots =\frac{\alpha^5-6\alpha^3+\alpha}{6(1+\alpha^2)^3} \end{align*}

We conclude \begin{align*} \arctan\left(\alpha e^x\right)&=\arctan(\alpha)+\frac{\alpha}{1+\alpha^2}x+\frac{\alpha(1-\alpha^2)}{2\left(1+\alpha^2\right)^2}x^2\\ &\qquad+\frac{\alpha^5-6\alpha^3+\alpha}{6(1+\alpha^2)^3}x^3+O(x^4)\tag{3} \end{align*}



Powers in logarithmic series:

In order to derive the coefficients of the logarithmic series \begin{align*} \ln\left(\arctan(\alpha e^x)\right)&=\sum_{n=1}^\infty\frac{(-1)^{n+1}}{n}(\arctan(\alpha e^x)-1)^n \end{align*} we write the expression (3) as \begin{align*} \arctan\left(\alpha e^x\right)&=a_0+a_1x+a_2x^2+a_3x^3+O(x^4) \end{align*} and we consider \begin{align*} \ln\left(\arctan(\alpha e^x)\right)&=\sum_{n=1}^\infty\frac{(-1)^{n+1}}{n}((a_0-1)+a_1x+a_2x^2+a_3x^3)^n+O(x^4)\tag{4} \end{align*}

We now set $A(x)=(a_0-1)+a_1x+a_2x^2+a_3x^3$ and extract the coeffcients of $x^0$ to $x^3$ from \begin{align*} \left(A(x)\right)^n&=((a_0-1)+a_1x+a_2x^2+a_3x^3)^n\\ &=\sum_{j=0}^n\binom{n}{j}(a_0-1)^j\left(a_1x+a_2x^2+a_3x^3\right)^{n-j}\\ &=\sum_{j=0}^n\binom{n}{j}(a_0-1)^j\sum_{k=0}^{n-j}\binom{n-j}{k}a_1^kx^k\left(a_2x^2+a_3x^3\right)^{n-j-k}\tag{5}\\ \end{align*}

Series expansion of logarithm:

We calculate using (6) the coefficients of $\ln(\arctan(\alpha e^x))$ in terms of $a_j, 0\leq j\leq 3$

\begin{align*} [x^0]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^0]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^0](a_0-1)^n\\ &=\ln(a_0-1)\\ [x^1]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^1]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^0]a_1n(a_0-1)^{n-1}\\ &=a_1\sum_{n=0}^\infty(-1)^n(a_0-1)^n\\ &=\frac{a_1}{a_0}\\ [x^2]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^2]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}\left(a_2n(a_0-1)^{n-1}+\frac{1}{2}n(n-1)a_1^2(a_0-1)^{n-2}\right)\\ &=\cdots\\ &=\left(a_2+\frac{a_1^2}{2}\frac{d}{da_0}\right)\left(\frac{1}{a_0}\right)\\ &=\frac{a_2}{a_0}-\frac{a_1^2}{2a_0^2}\\ [x^3]&\ln(\arctan(\alpha e^x))=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}[x^3]A(x)\\ &=\sum_{n=1}^\infty \frac{(-1)^{n+1}}{n}\left(na_3(a_0-1)^{n-1}+a_1a_2n(n-1)(a_0-1)^{n-2}\right.\\ &\qquad\qquad\qquad\qquad\left.+\frac{1}{6}n(n-1)(n-2)a_1^3(a_0-1)^{n-3}\right)\\ &=\cdots\\ &=\left(a_3+a_1a_2\frac{d}{da_0}+\frac{a_1^3}{6}\frac{d^2}{da_0^2}\right)\left(\frac{1}{a_0}\right)\\ &=\frac{a_3}{a_0}-\frac{a_1a_2}{a_0^2}+\frac{a_1^3}{3a_0^3}\tag{7} \end{align*}



Series expansion of $f(x)$:

Now it's time to harvest. We finally obtain with (3) and (7) respecting that $f(x)$ is odd

\begin{align*} \color{blue}{f(x)}&\color{blue}{=\ln\left(\arctan\left(\alpha e^x\right)\right)-\ln\left(\arctan\left(\alpha e^{-x}\right)\right)}\\ &=\cdots\\ &=\frac{2a_1}{a_0}x+2\left(\frac{a_3}{a_0}-\frac{a_1a_2}{a_0^2}+\frac{a_1^3}{3a_0^3}\right)x^3+O(x^5)\\ &\color{blue}{=\frac{2\alpha}{(1+\alpha^2)\arctan(\alpha)}x +\frac{\alpha}{3(1+\alpha^2)^3\arctan(\alpha)}\Bigg(\alpha^4-6\alpha^2+1}\\ &\qquad\color{blue}{\left.-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right)x^3+O(x^5)} \end{align*}

• Markus, while you are correct about $O(x^4)$, since we know $f(x)$ has odd-symmetry and that the even-order terms are zero, i think you can say this expansion is good to $O(x^5)$. – robert bristow-johnson Feb 20 '17 at 18:35
• @robertbristow-johnson: Yes, of course. Updated accordingly. :-) – Markus Scheuer Feb 20 '17 at 21:05
• Great effort! Trying to read this detailed and lengthy answer I couldn't see how you could isolate $O(x^4)$, in the equation (4), outside of the logarithm? The infinite series already includes every power of $x$, so what does the isolated $O(x^4)$ term mean there? – Fat32 Feb 20 '17 at 22:26
• Of course I got the feeling of what you want to mean there but then the proper notation could be something like this: $$\ln\left(\arctan(\alpha e^x)\right) ~=~ \sum_{n=1}^\infty\frac{(-1)^{n+1}}{n}((a_0-1)+a_1x+a_2x^2+a_3x^3 + O_1(x^4))^n ~=~ T_0 + T_1 x + T_2 x^2 + T_3 x^3 + O_2(x^4)$$ where I used $T$ to stay away from all of your other notation. And note that I have used $O_1$ and $O_2$ to distinguish between those two sets of coefficents. So now your equation (4) and this above line are not exactly the same. I don't think however it will effect any of your further progress. – Fat32 Feb 21 '17 at 16:18
• @Fat32: You might want to look at big-O notation – Markus Scheuer Feb 21 '17 at 17:00

Using Wolfram Alpha, $f'''(x)$ at $x=0$ evaluates to: \begin{align} \\ f'''(0) = & -\frac{6 \alpha^2}{(\alpha^2 + 1)^2 (\arctan(\alpha))^2} \ + \ \frac{2 \alpha}{(\alpha^2 + 1) \arctan(\alpha)} \\ & \quad \quad + \frac{16 \alpha^5}{(\alpha^2 + 1)^3 \arctan(\alpha)} \ + \ \frac{12 \alpha^4}{(\alpha^2 + 1)^3 (\arctan(\alpha))^2} \\ & \quad \quad \quad \quad - \frac{16 \alpha^3}{(\alpha^2 + 1)^2 \arctan(\alpha)} \ + \ \frac{4 \alpha^3}{(\alpha^2 + 1)^3 (\arctan(\alpha))^3} \\ \\ &= \frac{2 (\alpha^4 - 6 \alpha^2 + 1) \alpha}{(\alpha^2 + 1)^3 \arctan(\alpha)} + \frac{6 (\alpha^2 - 1) \alpha^2}{(\alpha^2 + 1)^3 (\arctan(\alpha))^2} + \frac{4 \alpha^3}{(\alpha^2 + 1)^3 (\arctan(\alpha))^3} \\ \end{align}

http://www.wolframalpha.com/input/?i=evaluate+d3%2Fdx3++(+ln+(arctan+(a+exp(+x)))+-+ln+(arctan(a+exp(-+x)))+)+at+x%3D0

We can also double check if this matches up with Markus's answer here.

His coefficient of $x^3$ comes out to be

$$\frac{\alpha}{3(1+\alpha^2)^3\arctan(\alpha)}\left(\alpha^4-6\alpha^2+1-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right).$$

If we multiply that by 6 and rearrange some factors we get:

$$\frac{2\alpha(\alpha^4-6\alpha^2+1)}{(1+\alpha^2)^3\arctan(\alpha)} - \frac{6\alpha^2(1-\alpha^2)}{(1+\alpha^2)^3(\arctan(\alpha))^2} + \frac{4 \alpha^3}{(1+\alpha^2)^3 (\arctan(\alpha))^3}$$

which matches up!

• Atul, it appears that your simplified answer is not consistent with Markus's answer at the math SE. it should be the case that $$f'''(x)\bigg|_{x\to0} \ = \ 3! \, a_3 = 6 a_3$$ i don't think every term in your f'''(0) is consistent with Markus. could be that Markus is wrong. – robert bristow-johnson Feb 20 '17 at 7:54
• @robertbristow-johnson I think they match up. – Atul Ingle Feb 20 '17 at 14:41
• they do now. i think Markus must have had an error. he did his answer the good old-fashioned way. – robert bristow-johnson Feb 20 '17 at 17:58
• Atul, you will get your bounty. but i explored the rules about bounty and they don't let me split it, but they let me award it twice, but one at a time. so since Markus has less rep than you here at dsp.se and since he grunged out an answer without the help of a computer, i'm awarding his bounty first. then i will place another bounty on this question and then i'll award it to you. it says i need to wait 23 hours. dunno who's gonna get my "check mark" yet. – robert bristow-johnson Feb 22 '17 at 7:56
• @robertbristow-johnson sorry for the late response. The coefficients are $2/3, -2/15, -16/945, -2/945$ for $\omega_0^2, \omega_0^4, \omega_0^6, \omega_0^8$ respectively. wolframalpha.com/input/… – Atul Ingle Feb 27 '17 at 1:12

The problem as posed in the question appears to have no closed-form solution. As mentioned in the question and shown in other answers, the result can be developed into a series, which can be accomplished by any symbolic math tool such as Mathematica. However, the terms become quite complicated and ugly, and it is unclear how good the approximation is when we include terms up to third order. Since we can't get an exact formula, it might be better to compute the solution numerically, which, unlike with the approximation, will give an (almost) exact result.

However, this is not what my answer is about. I suggest a different route which gives an exact solution by changing the problem formulation. After thinking about it for a while it turns out that it is the specification of the center frequency $\omega_0$ and the specification of the bandwidth as a ratio (or, equivalently, in octaves) which causes the mathematical intractability. There are two ways out of the dilemma:

1. specify the bandwidth of the discrete-time filter as a difference of frequencies $\Delta\omega=\omega_2-\omega_1$, where $\omega_1$ and $\omega_2$ are the lower and upper band edges of the discrete-time filter, respectively.
2. prescribe the ratio $\omega_2/\omega_1$, and instead of $\omega_0$ prescribe one of the two edge frequencies $\omega_1$ or $\omega_2$.

In both cases, a simple analytical solution is possible. Since it is desirable to prescribe the bandwidth of the discrete-time filter as a ratio (or, equivalently, in octaves), I'll describe the second approach.

Let's define the edge frequencies $\Omega_1$ and $\Omega_2$ of the continuous-time filter by

$$|H(j\Omega_1)|^2=|H(j\Omega_2)|^2=\frac12\tag{1}$$

with $\Omega_2>\Omega_1$, where $H(s)$ is the transfer function of a second-order band pass filter:

$$H(s)=\frac{\Delta\Omega s}{s^2+\Delta\Omega s+\Omega_0^2}\tag{2}$$

with $\Delta\Omega=\Omega_2-\Omega_1$, and $\Omega_0^2=\Omega_1\Omega_2$. Note that $H(j\Omega_0)=1$, and $|H(j\Omega)|<1$ for $\Omega\neq\Omega_0$.

We use the bilinear transform to map the edge frequencies $\omega_1$ and $\omega_2$ of the discrete-time filter to the edge frequencies $\Omega_1$ and $\Omega_2$ of the continuous-time filter. Without loss of generality we can choose $\Omega_1=1$. For our purposes the bilinear transform then takes the form

$$s = \frac{1}{\tan\left(\frac{\omega_1}{2}\right)}\frac{z-1}{z+1}\tag{3}$$

corresponding to the following relationship between continuous-time and discrete-time frequencies:

$$\Omega=\frac{\tan\left(\frac{\omega}{2}\right)}{\tan\left(\frac{\omega_1}{2}\right)}\tag{4}$$

From $(4)$ we obtain $\Omega_2$ by setting $\omega=\omega_2$. With $\Omega_1=1$ and $\Omega_2$ computed from $(4)$, we obtain the transfer function of the analog prototype filter from $(2)$. Applying the bilinear transform $(3)$, we get the transfer function of the discrete-time band pass filter:

$$H_d(z)=g\cdot\frac{z^2-1}{z^2+az+b}\tag{5}$$

with

\begin{align}g&=\frac{\Delta\Omega c}{1+\Delta\Omega c+\Omega_0^2c^2}\\a&=\frac{2(\Omega_0^2c^2-1)}{1+\Delta\Omega c+\Omega_0^2c^2}\\b&=\frac{1-\Delta\Omega c+\Omega_0^2c^2}{1+\Delta\Omega c+\Omega_0^2c^2}\\c&=\tan\left(\frac{\omega_1}{2}\right)\end{align}\tag{6}

Summary:

The bandwidth of the discrete-time filter can be specified in octaves (or, generally, as a ratio), and the parameters of the analog prototype filter can be computed exactly, such that the specified bandwidth is achieved. Instead of the center frequency $\omega_0$, we specify the band edges $\omega_1$ and $\omega_2$. The center frequency defined by $|H_d(e^{j\omega_0})|=1$ is an outcome of the design.

The necessary steps are as follows:

1. Specify the desired ratio of band edges $\omega_2/\omega_1$, and one of the band edges (which is of course equivalent to simply specifying $\omega_1$ and $\omega_2$).
2. Choose $\Omega_1=1$ and determine $\Omega_2$ from $(4)$. Compute $\Delta\Omega=\Omega_2-\Omega_1$ and $\Omega_0^2=\Omega_1\Omega_2$ of the analog prototype filter $(2)$.
3. Evaluate the constants $(6)$ to obtain the discrete-time transfer function $(5)$.

Note that with the more common approach where $\omega_0$ and $\Delta\omega=\omega_2-\omega_1$ are specified, the actual band edges $\omega_1$ and $\omega_2$ are an outcome of the design process. In the proposed solution, the band edges can be specified and $\omega_0$ is an outcome of the design process. The advantage of the latter approach is that the bandwidth can be specified in octaves and the solution is exact, i.e., the resulting filter has exactly the specified bandwidth in octaves.

Example:

Let's specify a bandwidth of one octave, and we choose the lower band edge as $\omega_1=0.2\pi$. This gives an upper band edge $\omega_2=2\omega_1=0.4\pi$. The band edges of the analog prototype filter are $\Omega_1=1$ and from $(4)$ (with $\omega=\omega_2$) $\Omega_2=2.2361$. This gives $\Delta\Omega=\Omega_2-\Omega_1=1.2361$ and $\Omega_0^2=\Omega_1\Omega_2=2.2361$. With $(6)$ we get for the discrete-time transfer function $(5)$

$$H_d(z)=0.24524 \cdot \frac{z^2-1}{z^2-0.93294z+0.50953}$$

which achieves exactly a bandwidth of 1 octave, and the specified band edges, as shown in the figure below:

Numerical solution of the original problem:

From the comments I understand that it is important to be able to exactly specify the center frequency $\omega_0$ for which $|H_d(e^{j\omega_0})|=1$ is satisfied. As mentioned before it is not possible to get an exact closed-form solution, and a series development produces quite unwieldy expressions.

For the sake of clarity I would like to summarize the possible options with their advantages and disadvantages:

1. specify the desired bandwidth as a frequency difference $\Delta\omega=\omega_2-\omega_1$, and specify $\omega_0$; in this case a simple closed-form solution is possible.
2. specify the band edges $\omega_1$ and $\omega_2$ (or, equivalently, the bandwidth in octaves, and one of the band edges); this also leads to a simple closed-form solution, as explained above, but the center frequency $\omega_0$ is an outcome of the design and cannot be specified.
3. specify the desired bandwidth in octaves and the center frequency $\omega_0$ (as asked in the question); no closed form solution is possible, nor is there (for the time being) any simple approximation. For this reason I think it's desirable to have a simple and efficient method for obtaining a numerical solution. This is what is explained below.

When $\omega_0$ is specified we use a form of the bilinear transform with a normalization constant that is different from the one used in $(3)$ and $(4)$:

$$\Omega=\frac{\tan\left(\frac{\omega}{2}\right)}{\tan\left(\frac{\omega_0}{2}\right)}\tag{7}$$

We define $\Omega_0=1$. Denote the specified ratio of band edges of the discrete-time filter as

$$r=\frac{\omega_2}{\omega_1}\tag{8}$$

With $c=\tan(\omega_0/2)$ we get from $(7)$ and $(8)$

$$r=\frac{\arctan(c\Omega_2)}{\arctan(c\Omega_1)}\tag{9}$$

With $\Omega_1\Omega_2=\Omega_0^2=1$, $(9)$ can be rewritten in the following form:

$$f(\Omega_1)=r\arctan(c\Omega_1)-\arctan\left(\frac{c}{\Omega_1}\right)=0\tag{10}$$

For a given value of $r$ this equation can be solved for $\Omega_1$ with a few Newton iterations. For this we need the derivative of $f(\Omega_1)$:

$$f'(\Omega_1)=c\left(\frac{r}{1+c^2\Omega_1^2}+\frac{1}{c^2+\Omega_1^2}\right)\tag{11}$$

With $\Omega_0=1$, we know that $\Omega_1$ must be in the interval $(0,1)$. Even though it's possible to come up with smarter initial solutions, it turns out that the initial guess $\Omega_1^{(0)}=0.1$ works well for most specs, and will result in very accurate solutions after only $4$ iterations of Newton's method:

$$\Omega_1^{(n+1)}=\Omega_1^{(n)}-\frac{f(\Omega_1^{(n)})}{f'(\Omega_1^{(n)})}\tag{12}$$

With $\Omega_1$ obtained with a few iterations of $(12)$ we can determine $\Omega_2=1/\Omega_1$ and $\Delta\Omega=\Omega_2-\Omega_1$, and and we use $(5)$ and $(6)$ to compute the coefficients of the discrete-time filter. Note that the constant $c$ is now given by $c=\tan(\omega_0/2)$.

Example 1:

Let's specify $\omega_0=0.6\pi$ and a bandwidth of $0.5$ octaves. This corresponds to a ratio $r=\omega_2/\omega_1=2^{0.5}=\sqrt{2}=1.4142$. With an initial guess of $\Omega_1=0.1$, $4$ iterations of Newton's method resulted in a solution $\Omega_1=0.71$, from which the coefficients of the discrete-time can be computed as explained above. The figure below shows the result:

The filter was calculated with this Matlab/Octave script:

% specifications
bw = 0.5;    % desired bandwidth in octaves
w0 = .6*pi;  % resonant frequency

r = 2^(bw);  % ratio of band edges
W1 = .1;     % initial guess (works for most specs)
Nit = 4;     % # Newton iterations
c = tan(w0/2);

% Newton
for i = 1:Nit,
f = r*atan(c*W1) - atan(c/W1);
fp = c * ( r/(1+c^2*W1^2) + 1/(c^2+W1^2) );
W1 = W1 - f/fp
end

W1 = abs(W1);
if (W1 >= 1), error('Failed to converge. Reduce value of initial guess.'); end

W2 = 1/W1;
dW = W2 - W1;

% discrete-time filter
scale = 1 + dW*c + W1*W2*c^2;
b = ( dW*c/scale) * [1,0,-1];
a = [1, 2*(W1*W2*c^2-1)/scale, (1-dW*c+W1*W2*c^2)/scale ];



Example 2:

I add another example to show that this method can also deal with specifications for which most approximations will give non-sensical results. This is often the case when the desired bandwidth and the resonant frequency are both large. Let's design a filter with $\omega_0=0.95\pi$ and $bw=4$ octaves. Four iterations of Newton's method with an initial guess $\Omega_1^{(0)}=0.1$ result in a final value of $\Omega_1=0.00775$, i.e., in a bandwidth of the analog prototype of $\log_2(\Omega_2/\Omega_1)=\log_2(1/\Omega_1^2)\approx 14$ octaves. The corresponding discrete-time filter has the following coefficients and its frequency response is shown in the plot below:

b = 0.90986*[1,0,-1];
a = [1.00000   0.17806  -0.81972];


The resulting half power band edges are $\omega_1=0.062476\pi$ and $\omega_2 = 0.999612\pi$, which are indeed exactly $4$ octaves (i.e., a factor of $16$) apart.

• two initial comments (i haven't read this through, yet, Matt): first, i am interested in log frequency more so than linear frequency. for the analog BPF (or the digital BPF with resonant frequency much lower than Nyquist), there is perfect symmetry about the resonant frequency. – robert bristow-johnson Feb 20 '17 at 18:08
• and the second comment is this, while i thank you for apparently sticking with the notation of $s=j\Omega$ and $z=e^{j\omega}$, i wish you would stick to the notation that the analog and digital resonant frequencies are $\Omega_0$ and $\omega_0$, respectively, and the analog upper and lower bandedges are $\Omega_U$ and $\Omega_L$ respectively and likewise for the digital bandedges: $\omega_U$ and $\omega_L$. we know that, in log frequency, half of the bandwidth is above $\Omega_0$ and half is below. but, due to warping, that is not exactly true for the digital BPF filter. – robert bristow-johnson Feb 20 '17 at 18:14
• as i read this more, it is important to me that the resonant frequency gets mapped through the bilinear transform exactly. so i understand this approach, Matt, but i want to stick with exact mapping of $\omega_0$ and then tweak the $BW$ until $bw$ is what is specified. – robert bristow-johnson Feb 20 '17 at 18:33
• @robertbristow-johnson: OK, fair enough, you want an exact specification of $\omega_0$. That's possible if you specify $\Delta\omega$ as a linear difference (which you don't want, I understand). A neat solution is not possible with specified $\omega_0$ AND a bandwidth in octaves. – Matt L. Feb 20 '17 at 19:34
• @robertbristow-johnson: I've added a very simple numerical solution to my answer (4 Newton iterations). – Matt L. Feb 21 '17 at 13:24

okay, i promised to put up bounty and i will keep my promise. but i have to confess that i might renege a little bit on being satisfied with just the third derivative of $f(x)$. what i really want are the two coefficients for $g(y)$.

so i didn't realize that there was this Wolfram language as an alternative to mathematica or Derive and i didn't realize it could so easily compute the third derivative and simplify the expression.

and this Markus guy at the math SE posted this answer (which i thought was gonna have to be the amount of grunge i thought would be needed).

\begin{align*} y = f(x) &= \ln\left(\arctan\left(\alpha e^x\right)\right) - \ln\left(\arctan\left(\alpha e^{-x}\right)\right)\\ \\ &\approx a_1 x \ + \ a_3 x^3 \\ \\ &=\frac{2\alpha}{(1+\alpha^2)\arctan(\alpha)} x + \frac{\alpha}{3(1+\alpha^2)^3\arctan(\alpha)}\Bigg(\alpha^4-6\alpha^2+1\\ &\qquad\left.-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right) x^3 \\ \\ \end{align*}

so i put together the third-order approximation to the inverse:

\begin{align*} x = g(y) &\approx b_1 y \ + \ b_3 y^3 \\ \\ &= \frac{1}{a_1} y \ - \ \frac{a_3}{a_1^4} y^3 \\ \\ &= \frac{(1+\alpha^2)\arctan(\alpha)}{2\alpha} y - \frac{(1+\alpha^2)(\arctan(\alpha))^3}{48 \alpha^3}\Bigg(\alpha^4-6\alpha^2+1 \\ &\qquad\left.-\frac{3\alpha(1-\alpha^2)}{\arctan(\alpha)}+\frac{2\alpha^2}{\left(\arctan(\alpha)\right)^2}\right) y^3 \\ \\ &= \frac{(1+\alpha^2)\arctan(\alpha)}{2\alpha} y - \frac{(1+\alpha^2)(\arctan(\alpha))^3}{48 \alpha}\Bigg(\alpha^2-6+\alpha^{-2} \\ &\qquad\left.-\frac{3(1-\alpha^2)}{\alpha\arctan(\alpha)}+\frac{2}{\left(\arctan(\alpha)\right)^2}\right) y^3 \\ \\ &= y \left(\arctan(\alpha)\frac{\alpha + \alpha^{-1}}{2}\right) \Bigg(1 \ + \\ & \left((\arctan(\alpha))^2\left(1-\frac{\alpha^2+\alpha^{-2}}{6} \right) - \arctan(\alpha)\frac{\alpha-\alpha^{-1}}{2}-\frac{1}{3}\right) \frac{y^2}{4} \Bigg) \\ \\ \end{align*}

i was kinda hoping someone else would do this. recall $y=f(x) \triangleq \ln(2) \, bw$, $g(y) = x \triangleq \frac{\ln(2)}{2}BW$ and $\alpha \triangleq \tan(\omega_0/2)$

\begin{align*} x = g(y) &\approx y \left(\arctan(\alpha)\frac{\alpha + \alpha^{-1}}{2}\right) \Bigg(1 \ + \\ & \left((\arctan(\alpha))^2\left(1-\frac{\alpha^2+\alpha^{-2}}{6} \right) - \arctan(\alpha)\frac{\alpha-\alpha^{-1}}{2}-\frac{1}{3}\right) \frac{y^2}{4} \Bigg) \\ \\ \frac{\ln(2)}{2}BW &\approx (\ln(2) bw) \left(\arctan(\alpha)\frac{\alpha + \alpha^{-1}}{2}\right) \Bigg(1 \ + \\ & \left((\arctan(\alpha))^2\left(1-\frac{\alpha^2+\alpha^{-2}}{6} \right) - \arctan(\alpha)\frac{\alpha-\alpha^{-1}}{2}-\frac{1}{3}\right) \frac{(\ln(2) bw)^2}{4} \Bigg) \\ \\ \end{align*}

i have three convenient trig identities:

$$\frac12 \left(\alpha + \alpha^{-1} \right) = \frac12 \left(\tan(\omega_0/2) + \frac{1}{\tan(\omega_0/2)} \right) = \frac{1}{\sin(\omega_0)} \\$$

$$\frac12 \left(\alpha - \alpha^{-1} \right) = \frac12 \left(\tan(\omega_0/2) - \frac{1}{\tan(\omega_0/2)} \right) = -\frac{1}{\tan(\omega_0)} \\$$

$$\frac12 \left(\alpha^2 + \alpha^{-2} \right) = \frac12 \left(\tan^2(\omega_0/2) + \frac{1}{\tan^2(\omega_0/2)} \right) = \frac{1}{\sin^2(\omega_0)} + \frac{1}{\tan^2(\omega_0)} \\ = \frac{2}{\sin^2(\omega_0)} - 1$$

"finally" we got:

$$BW \approx bw \, \frac{\omega_0}{\sin(\omega_0)} \left(1 \ + \ \frac{(\ln(2))^2}{24} \left( 2(\omega_0^2 - 1) - \left(\frac{\omega_0}{\sin(\omega_0)}\right)^2 + 3\frac{\omega_0}{\tan(\omega_0)} \right) (bw)^2 \right)$$

this ain't so bad. fits on a single line. if someone sees an error or a good way to simplify further, please lemme know.

with the power series approximation from the comment above,

$$BW \approx bw \, \frac{\omega_0}{\sin(\omega_0)} \left(1 \ + \ (\ln(2))^2 \left( \tfrac{1}{36}\omega_0^2 - \tfrac{1}{180}\omega_0^4 - \tfrac{2}{2835}\omega_0^6 \right) (bw)^2 \right)$$

• also, i am not sure that Atul's answer for $f'''(0)$ and Markus's answer for $a_3$ are consistent. i wonder if someone might be able to straight that out in an answer that could get in on the bounty. – robert bristow-johnson Feb 20 '17 at 6:54
• I also found out about Wolfram's cloud notebook which is like Mathematica in your webbrowser. Go to sandbox.open.wolframcloud.com/app and type in 6*SeriesCoefficient[ Series[Log[ArcTan[a E^x]] - Log[ArcTan[a/E^x]],{x,0,5}],3] – Atul Ingle Feb 20 '17 at 15:07
• @AtulIngle, i incorporated the Markus's corrections into the inverse function. would you mind checking the result for $g(y)$? – robert bristow-johnson Feb 20 '17 at 18:28
• i would appreciate it if someone would check my substitution back to $g(y)$, particularly the factor that multiplies $y^2$. very soon i will return $\alpha$ back to $\tan(\omega_0/2)$ which will cause a whole other simplification and form. but i will hold off a little in case someone tells me my simplifications above are wrong. – robert bristow-johnson Feb 20 '17 at 19:57
• @ robert bistow-johnson I checked your final expression for g(y) using Mathematica, it looks right. – Atul Ingle Feb 20 '17 at 23:27

so here are some quantitative results. i plotted spec'd bandwidth $bw$ for the digital filter on the x-axis and the resulting digital bandwidth on the y-axis. there are five plots from green to red representing the resonant frequency $\omega_0$ normalized by Nyquist:

$\frac{\omega_0}{\pi} =$ [0.0002 0.2441 0.4880 0.7320 0.9759]

so the resonant frequency goes from nearly DC to nearly Nyquist.

here is no compensation (or pre-warping) at all for bandwidth:

here is the simple first-order compensation that the Cookbook has done all along:

here is the third-order compensation that we just solved here:

what we want is for all of the lines to lie directly on the main diagonal.

i had made a mistake in the third-order case and corrected it in this revision. it does look like the third-order approximation to $g(y)$ is a bit better than the first-order approximation for small $bw$.

so i diddled with the coefficient of the 3rd-order term (i wanna leave the 1st-order term the same), lessening it's effect. this is from multiplying just the 3rd-order term by 50%:

this is reducing it to 33%:

and this is reducing the 3rd-order term to 25%:

since the object of an inverse function is to undo the specified function, the point of this whole thing is to get the curves of the composite function to lie as close to the main diagonal as possible. it's not too bad for up to 75% Nyquist for resonant frequency $\omega_0$ and 3 octaves bandwidth $bw$. but not so much better to really make it worth it in the "coefficient cooking" code that is executed whenever the user turns a knob or slides a slider.

• How can the bandwidth become negative in the second and third plot?? – Matt L. Feb 21 '17 at 19:53
• it can't, which is why i am so far unimpressed with this third-order approximation to the real $x=g(y)$ which is the inverse function of $$f(x) = \ln\left(\frac{\arctan(\alpha e^x)}{\arctan(\alpha e^{-x})}\right)$$ i don't think that third-order approximation is an improvement over the first-order approximation that has existed for a couple decades. so what is plotted is $$f( \hat{g}(y) )$$ where $\hat{g}(y)$ is the approximation to the true inverse $g(y)$ where $$y = f( g(y) )$$ because $f(x)$ is bipolar (even though negative bandwidth is nonsensical) $f(\hat{g}(y))$ can go negative. – robert bristow-johnson Feb 21 '17 at 23:37
• oh, @MattL. the fact that $f(x)$ passes through the origin should not surprize you even if bandwidth is never really negative. that bandwidth mapping function is odd symmetry, so the first and second plot do not surprize me at all. but the third plot is disappointing. – robert bristow-johnson Feb 21 '17 at 23:40
• I was just wondering why you plotted the curves for negative bandwidths. But anyway, if I'm not mistaken then the series you use is a kind of Taylor series expansion at $bw=0$, right? So why would you even expect it to approximate the real behavior well at larger bandwidths if you only use two terms? – Matt L. Feb 22 '17 at 6:52
• i just wanted to make sure the functions are odd-symmetry and go through the origin real nicely. yes, this is all about Taylor (or more specificly, Maclaurin) series. you will notice, @MattL., that i think one term does rather nicely for all of the resonant frequencies that ain't terribly close to Nyquist. leaving the linear term unchanged, i diddled a little with the third-order term a little (stay tuned, i'll show the results) and it does pretty well. but not so much better than the first-order that i think i should bother changing it in the Cookbook. – robert bristow-johnson Feb 22 '17 at 7:27