Computing $\rho$

Our goal is to find the allpass coefficient $\rho$ such that the frequency mapping

$$a(\omega )= \mbox{angle}\left\{{\cal A}_{-\rho }(e^{j\omega }) \right\}$$ (5)

best approximates the Bark scale $b(\omega )$ for a given sampling rate f_s. (Note that the frequencies $\omega$, $a(\omega )$, and $b(\omega )$ are all expressed in radians per sample, so that a frequency of half of the sampling rate corresponds to a value of $\pi$.)

Using squared frequency errors to gauge the fit between $a(\omega )$ and its Bark warped counterpart, the optimal mapping parameter $\rho ^*$ may be written as

$$\rho ^*= \mbox{Arg}\left[\min_{\rho }\left\{\left\Vert\,a(\omega )- b(\omega )\,\right\Vert\right\}\right],$$ (6)

where $\left\Vert\,\cdot\,\right\Vert$ represents the L₂ norm. (We use the superscript `$\ast$' to denote optimality in some sense.) Unfortunately, the frequency error

$$\epsilon _{\mbox{\tiny A}}\mathrel{\stackrel{\mathrm{\Delta}}{=}}a(\omega )- b(\omega )$$ (7)

is nonlinear in $\rho$, and its norm is not easily minimized directly. It turns out, however, that a related error,

$$\epsilon _{\mbox{\tiny C}}\mathrel{\stackrel{\mathrm{\Delta}}{=}}e^{ja(\omega )}- e^{jb(\omega )},$$ (8)

has a norm which is more amenable to minimization. The first issue we address is how the minimizers of $\left\Vert\,\epsilon _{\mbox{\tiny A}}\,\right\Vert$ and $\left\Vert\,\epsilon _{\mbox{\tiny C}}\,\right\Vert$ are related.

Figure 2: Frequency Map Errors

Denote by $\zeta$ and $\beta$ the complex representations of the frequencies $a(\omega )$ and $b(\omega )$ on the unit circle,

$$\zeta = e^{ja(\omega )}, \qquad \beta = e^{jb(\omega )}.$$ (9)

As seen in Fig.2, the absolute frequency error $\vert\epsilon _{\mbox{\tiny A}}\vert$ is the arc length between the points $\zeta$ and $\beta$, whereas $\vert\epsilon _{\mbox{\tiny C}}\vert$ is the chord length or distance:

$$\vert\epsilon _{\mbox{\tiny C}}\vert = 2\sin(\vert\epsilon _{\mbox{\tiny A}}\vert/2).$$ (10)

The desired arc length error $\epsilon _{\mbox{\tiny A}}$ gives more weight to large errors than the chord length error $\epsilon _{\mbox{\tiny C}}$; however, in the presence of small discrepancies between $\zeta$ and $\beta$, the absolute errors are very similar,

$$\vert\epsilon _{\mbox{\tiny C}}\vert \approx \vert\epsilon _{\mbox{\tiny A}}\vert, \quad \mbox{when } \vert\epsilon _{\mbox{\tiny A}}\vert\ll1.$$ (11)

Accordingly, essentially the same $\rho ^*$ results from minimizing $\left\Vert\,\epsilon _{\mbox{\tiny A}}\,\right\Vert$ or $\left\Vert\,\epsilon _{\mbox{\tiny C}}\,\right\Vert$ when the fit is uniformly good over frequency.

The error $\epsilon _{\mbox{\tiny C}}$ is also nonlinear in the parameter $\rho$, and to find its norm minimizer, an equation error is introduced, as is common practice in developing solutions to nonlinear system identification problems [17]. Consider mapping the frequency $z=e^{j\omega }$ via the allpass transformation ${\cal A}_{-\rho }(z)$,

$$\zeta = {z- \rho \over 1- z\rho }.$$ (12)

Now, multiply Eq.(12) by the denominator $(1-z\rho )$, and substitute $\zeta =\beta +\epsilon _{\mbox{\tiny C}}$ from Eq.(8), to get

$$(\beta + \epsilon _{\mbox{\tiny C}}) (1 - z\rho ) = z- \rho .$$ (13)

Rearranging terms, we have

$$(\beta - z) - (\beta z- 1) \rho = \epsilon _{\mbox{\tiny E}},$$ (14)

where $\epsilon _{\mbox{\tiny E}}$ is an equation error defined by

$$\epsilon _{\mbox{\tiny E}}\mathrel{\stackrel{\mathrm{\Delta}}{=}}(z\rho - 1) \epsilon _{\mbox{\tiny C}}.$$ (15)

Figure 3: Geometric Interpretation of Equation Error

Referring to Fig.3, note that the equation error is the difference between the unit circle chord $(\beta -z)$ and the $\rho$ circle chord $(\beta z-1)\rho$. The average of the input and Bark warped angles $\left[b(\omega )+\omega \right]/2$ bisects both these chords, and therefore the chords are parallel. The equation error may then be interpreted as the difference in chord lengths, rotated to the angle $\left[b(\omega )+\omega \right]/2 + \pi/2$. This suggests defining a rotated equation error which is real valued. Multiply Eq.(14) by $(z\beta )^{-{1\over2}}/2j$ to obtain

$$\sin\left[\frac{b(\omega )-\omega }{2}\right] - \rho \sin\left[\frac{b(\omega )+\omega }{2}\right] = \epsilon _{\mbox{\tiny R}},$$ (16)

where the rotated equation error $\epsilon _{\mbox{\tiny R}}$ is defined by

$$\epsilon _{\mbox{\tiny R}}\mathrel{\stackrel{\mathrm{\Delta}}{=}}{1\over 2j} (z\beta )^{-{1\over2}} \epsilon _{\mbox{\tiny E}} = {1\over 2j} (z\beta )^{-{1\over2}} (z\rho - 1) \epsilon _{\mbox{\tiny C}}.$$ (17)

The rotated equation error $\epsilon _{\mbox{\tiny R}}$ is linear in the unknown $\rho$, and its norm minimizer is easily expressed in closed form. Denote by $\omega_{k}, k=1,\ldots,K$ a set of frequencies corresponding to Bark frequencies $b(\omega_{k})$, and by $\mbox{\boldmath$d$}$ and $\mbox{\boldmath$s$}$ the columns

$$\mbox{\boldmath$d$} \left[\matrix{ \sin\{[b(\omega_{1})-\omega_{1}]/2\} \cr \vdots \cr \sin\{[b(\omega_{K})-\omega_{K}]/2\} }\right]$$ = (18)

$$\mbox{\boldmath$s$} \left[\matrix{ \sin\{[b(\omega_{1})+\omega_{1}]/2\} \cr \vdots \cr \sin\{[b(\omega_{K})+\omega_{K}]/2\} }\right].$$ = (19)

Then Eq.(16) becomes

$$\mbox{\boldmath$d$}- \mbox{\boldmath$s$}\rho = \mbox{\boldmath$\epsilon$}_{\mbox{\tiny R}},$$ (20)

where $\mbox{\boldmath$\epsilon$}_{\mbox{\tiny R}}$ is a column of rotated equation errors. Eq.(20) is now in the form of a standard least squares problem. As is well known [17,34], the solution $\rho ^*$ which minimizes the weighted sum of squared errors, $\mbox{\boldmath$\epsilon$}_{\mbox{\tiny R}}^\top \mbox{\boldmath$V$}\mbox{\boldmath$\epsilon$}_{\mbox{\tiny R}}$, the matrix $\mbox{\boldmath$V$}$ being an arbitrary positive-definite weighting, may be obtained by premultiplying both sides of Eq.(20) by $\mbox{\boldmath$s$}^\top \mbox{\boldmath$V$}$ and solving for $\rho$, noting that at $\rho =\rho ^*$, $\mbox{\boldmath$s$}^\top \mbox{\boldmath$V$}\mbox{\boldmath$\epsilon$}_{\mbox{\tiny R}}=0$ by the orthogonality principle. Doing this yields the optimal weighted least-squares conformal map parameter

$$\rho ^*= {\mbox{\boldmath$s$}^\top \mbox{\boldmath$V$}\mbox{\boldmath$d$}\over \mbox{\boldmath$s$}^\top \mbox{\boldmath$V$}\mbox{\boldmath$s$}} .$$ (21)

If the weighting matrix $\mbox{\boldmath$V$}$ is diagonal with kth diagonal element $v(\omega_{k})>0$, then the weighted least-squares solution Eq.(21) reduces to

$$\begin{eqnarray*} \rho ^*&=& \frac{\sum_{k=1}^K v(\omega _k) \sin\left[\frac{b(\omega_{k})-\omega_{k}}{2}\right] \sin\left[\frac{b(\omega_{k})+\omega_{k}}{2}\right] }{\sum_{k=1}^K v(\omega _k)\sin^2\left[\frac{b(\omega_{k})+\omega_{k}}{2}\right]} \\ &=& \frac{\sum_{k=1}^{K} v(\omega _k) \left\{\cos\left[b(\omega_{k})\right]- \cos(\omega_{k})\right\}}{% \sum_{k=1}^{K} v(\omega_{k}) \left\{\cos\left[b(\omega_{k}) + \omega_{k}\right]- 1\right\}}, \end{eqnarray*}$$

where we have used Equations Eq.(18) and Eq.(19), and the trigonometric identities

$$\begin{eqnarray*} \cos(A+B) - \cos(A-B) &=& -2\sin(A)\sin(B) \\ \cos(A)-1 &=& -2\sin^2(A/2) \end{eqnarray*}$$

to simplify the numerator and denominator, respectively.

It remains to choose a weighting matrix $\mbox{\boldmath$V$}$. Recall that we initially wanted to minimize the sum of squared chord-length errors $\mbox{\boldmath$\epsilon$}_{\mbox{\tiny C}}^\top \mbox{\boldmath$\epsilon$}_{\mbox{\tiny C}}$. The rotated equation error $\epsilon _{\mbox{\tiny R}}$ is proportional to the chord length error $\epsilon _{\mbox{\tiny C}}$, viz. Eq.(17), and it is easily verified that when $\mbox{\boldmath$V$}$ is diagonal with kth diagonal element

$$v(\omega_{k}) = {1\over 1 + \rho ^2 - 2\rho \cos\omega_{k}},$$ (22)

the chord length error and the weighted equation error coincide. Thus, with this diagonal weighting matrix, the solution Eq.(21) minimizes the chord-length error norm.

Note that the desired weighting depends on the unknown map parameter $\rho$. To overcome this difficulty, we suggest first estimating $\rho ^*$ using $\mbox{\boldmath$V$} = \mbox{\boldmath$I$}$, where $\mbox{\boldmath$I$}$ denotes the identity matrix, and then computing $\rho ^*$ using the weighting Eq.(22) based on the unweighted solution. This is analogous to the Steiglitz-McBride algorithm for converting an equation-error minimizer to the more desired ``output-error'' minimizer using an iteratively computed weight function [16].