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:
- 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.
- 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:
- 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$).
- 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)$.
- 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:
- 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.
- 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.
- 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.
