Sunday, July 22, 2018

audio - 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$.



Answer



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


No comments:

Post a Comment

periodic trends - Comparing radii in lithium, beryllium, magnesium, aluminium and sodium ions

Apparently the of last four, $\ce{Mg^2+}$ is closest in radius to $\ce{Li+}$. Is this true, and if so, why would a whole larger shell ($\ce{...