Ttgjyuk

I am currently reading and trying to understand a paper (Kulkarni and Colburn, 2004) that utilizes system identification methods to approximate head-related transfer functions.

The general approach is to

1. Compute an autoregressive (all-pole-) estimate of the transfer function using the autocorrelation method for linear prediction.


2. Use the AR estimate as a starting point to compute a pole-zero-representation of the system transfer function iteratively.


3. Evaluate the result of the estimation process on a logarithmic scale (Error in dB).

For the iterative procedure, the authors are proposing a cost function

$hatC = frac12piint_-pi^pi|H(e^jomega)A(e^jomega) - B(e^jomega)|^2 domega $,

where $hatH$ is the system transfer function, $A$ is the DTFT of the recursive coefficients and $B$ the DTFT of the transversal coefficients.
I understand this approach originates from this paper (Kalman, 1958).

This cost function is then extended for the iterative process as

$hatC_i = frac12piint_-pi^pi|fracH(e^jomega)A_i(e^jomega)A_i-1(e^jomega) - fracB_i(e^jomega)A_i-1(e^jomega)|^2 domega $,

where the index $i$ denotes variables corresponding to the $i$th iteration.
The iterative modification originates from this paper (Steiglitz and McBride, 1965).

In order to find a solution on a decibel scale, a weighting function $W$ is introduced:

$hatC_i = frac12piint_-pi^pi |W_i(e^jomega)|^2 |fracH(e^jomega)A_i(e^jomega)A_i-1(e^jomega) - fracB_i(e^jomega)A_i-1(e^jomega)|^2 domega $,

which is introduced in the paper as

$W_i(e^jomega) = fraclogleft( H(e^jomega)right) - logleft( fracB_i-1(e^jomega)A_i-1(e^jomega)right)$.

I understand this weighting function is the squared error in logarithmic scale between true and approximated transfer function, divided by the first cost function.

However, i have trouble understanding the process of arriving at the iterative cost function for several reasons. I would like to ask the following questions:

• Why is the first cost function $hatC$ preferred to, say, $frac12piint_-pi^pi |H(e^jomega) - B(e^jomega)/A(e^jomega)|^2$ in the first place? what does it do?


• In the iterative cost function $hatC_i$, what is the purpose of dividing by the recursive proportion of the last transfer function estimate?


• For what reason is a weighting function introduced for logarithmic error minimization, rather than just using its numerator as a cost function directly?

I would really appreciate any help or pointers into the right direction.

The chosen cost function is the mean squared error, i.e., the integral over a squared magnitude of the difference between frequency responses. The function

$$E(e^jomega)=H(e^jomega)-fracB(e^jomega)A(e^j omega)tag1$$

depends on frequency, so you can't minimize it directly, unless you want to minimize it for exactly one frequency $omega$, which is of course pointless. You can choose several error measures depending on $E(e^jomega)$ given in $(1)$. Two common choices are

$$varepsilon_1=max_omegaW(omega)|E(e^jomega)|tag2$$

and

$$varepsilon_2=int_0^piW(omega)|E(e^jomega)|^2domegatag3$$

with some positive weighting function $W(omega)$. Note that $varepsilon_1$ and $varepsilon_2$ given by $(2)$ and $(3)$ do not depend on frequency.

The authors of the paper you refer to chose to minimize the weighted mean square error given by $(3)$. However, instead of using the linear difference $(1)$, they chose to minimize the average logarithmic difference, i.e., the average error on a dB scale.

The problem with the minimization of $(3)$ is that for IIR filters, it results in a non-linear optimization problem, which is much harder to solve directly, and which might also have locally optimal solutions that are far from the global optimum. The cost function $hatC$ in your question is linear in the filter coefficients. The point of the iteration is now to solve a sequence of linear minimization problems (which is simple, just solve a system of linear equations) in order to compute the solution of the originally non-linear optimization problem (the minimization of $(3)$). Note that if convergence is achieved, then $A_i-1(e^jomega)=A_i(e^jomega)$, so the cost function $hatC_i$ is identical to the cost function of the original non-linear problem. Yet, only linear minimization problems are solved in each iteration.

I'm not sure I completely understand your last question, but the weighting function is there to change to problem from minimizing the average squared difference between frequency responses to the average squared differences between the logarithms of the magnitude responses. The necessary weighting function is unknown, but if the procedure converges - note that there is no guarantee that it does in all cases - the final weight function is such that the mean squared logarithmic error is minimized. Note that this latter problem is highly non-linear, whereas the proposed procedure only solves linear subproblems.