Estimating characteristic phase and delay from broadband interaural time difference tuning curves

Abstract

Characteristic delay and characteristic phase are shape parameters of interaural time difference tuning curves. The standard procedure for the estimation of these parameters is based on the measurement of delay curves measured for tonal stimuli with varying frequencies. Common to all procedures is the detection of a linear behavior of the phase spectrum. Hence a reliable estimate can only be expected if sufficiently many relevant frequencies are tested. Thus, the estimation precision depends on the given bandwidth. Based on a linear model, we develop and implement methods for the estimation of characteristic phase and delay from a single broadband tuning curve. We present two different estimation algorithms, one based on a Fourier-analytic interpretation of characteristic delay and phase, and the other based on mean square error minimization. Estimation precision and robustness of the algorithms are tested on artificially generated data with predetermined characteristic delay and phase values, and on sample data from electrophysiological measurements in birds and in mammals. Increasing the signal-to-noise ratio or the bandwidth increases the estimation accuracy of the algorithms. Frequency band location and strong rectification also affect the estimation accuracy. For realistic bandwidths and signal-to-noise ratios, the minimization algorithm reliably and robustly estimates characteristic delay and phase and is superior to the Fourier-analytic method. Bandwidth-dependent significance thresholds allow to assess whether the estimated characteristic delay and phase values are meaningful shape parameters of a measured tuning curve. These thresholds also indicate the sampling rates needed to obtain reliable estimates from interaural time difference tuning curves.

Appendices

Appendix A: A linear weighting model

We assumed that the input neurons in the linear model are tuned exclusively to their best frequency. The jth input neuron ascends only to the integrator neuron if its best frequency ω _j is contained in the stimulus. Therefore, the tonal ITD tuning curve of the jth input neuron is described by the product of the delta function, an amplituded function and a cosine term (see Eq. (1)).

However, the approache by the delta function does not correspond to the cochlear frequency tuning. Neurons also respond to frequencies in a neighborhood of their best frequency. Thus Eq. (1) is idealized and somewhat unrealistic. To include such response behavior it is possible to replace the delta function by a weighting function m _j (⋅). For instance, the Gaussian function with mean ω _j and the half of the half width as standard deviation can be used as weighting function. Consequently, the response of the jth input neuron to a tonal stimulus can be modeled as

$$ \text{TC}_{j}(t,\omega)= m_{j}(\omega)\cdot a(\omega)\cdot \cos\left(2\pi[\omega(t-\text{CD})-\text{CP}]\right). $$

(23)

Due to the fact that tonal as well as broadband ITD tuning curves are real valued it is necessary to assume that the weighting function m _j (⋅) is positive and even, i.e. m _j (ω)≥0 and m _j (−ω)=m _j (ω) for all j∈{1⋯M}.

The response of the integrator neuron to a white noise stimuli is again given by Eq. (4). By the use of formula (23) and the identity $\displaystyle K(\omega ):= \sum\limits_{j=1}^{M} m_{j}(\omega )$ we obtain that

$$\begin{array}{@{}rcl@{}} \text{bTC}(t)\! &\,=\,&\! \sum\limits_{k=1}^{L}\sum\limits_{j=1}^{M}\text{TC}_{j}(t,\omega_{k})\\ &=& \sum\limits_{k=1}^{L}\sum\limits_{j=1}^{M} m_{j}(\omega_{k})a(\omega_{k})\cos\left(2\pi\!\left[\omega_{k}(t\,-\,\text{CD})\,-\,\text{CP}\right]\right)\\ &=& \sum\limits_{k=1}^{L} K(\omega_{k})a(\omega_{k})\cos\left(2\pi\left[\omega_{k}(t-\text{CD})-\text{CP}\right]\right)\\ &=& \sum\limits_{k=1}^{L} {\Delta}_{\omega} c(\omega_{k})\cos(2\pi\left[\omega_{k}(t-\text{CD})-\text{CP}\right]) \end{array} $$

(24)

where $\displaystyle c(\omega )= {\Delta }_{\omega }^{-1} a(\omega )K(\omega )$. As we have seen in Section 2.2 Eq. (24) can be interpretated as a Riemann sum approximation of the integal in Eq. (5). Therefore, the Hilbert Transform based estimation algorithm as well as the MSE fit algorithm is also applicable to a model in volving a frequency weighting function instead of the delta function.

Appendix B: The maximal absolute value of $\mathcal {I}_{+}$

In the Section 2.2.1 we claimed that the maximal absolute value of the function $\mathcal {I}$ is at the time t=CD. In this section we will proof this statement. However, we initially show the following theorem.

Theorem 1

Let c be a non-negative function with ${\int }_{-\infty }^{\infty }|c(\omega )|\:\mathrm {d}\omega <\infty$ and c ≠ 0 then the function $f:\mathbb {R} \rightarrow \mathbb {C}$ given by

$$ f(x) = {\int}_{-\infty}^{\infty}c(\omega) e^{2\pi i\omega (x-\text{CD})}\:\mathrm{d}\omega $$

(25)

has an unique maximal absolute value at x = CD. In particular

$$\lim\limits_{x\to \infty}f(x)=0.$$

Proof

At first we show the existence of the maximal absolute value x=CD. Without loss of generality we can set the CD to zero. Furthermore, we define $b(\omega ):=\sqrt {c(\omega )}$. Then

$$\begin{array}{@{}rcl@{}} f(x)\!\! &=&\!\! \!{\int}_{-\infty}^{\infty}|c(\omega)| e^{2\pi i\omega x}\:\mathrm{d}\omega \,=\,\! {\int}_{-\infty}^{\infty}\!|b(\omega)|^{2} e^{2\pi i\omega x}\:\mathrm{d}\omega\\ &=& {\int}_{-\infty}^{\infty}b(\omega)\overline{\left(M_{-x}b\right)(\omega)}\:\mathrm{d}\omega = \left<b,M_{-x}b\right>\end{array} $$

where (M _x b)(ω)=e ^2πiωx b(ω) for all $\omega \in \mathbb {R}$. By applying Plancherel’s Theorem for functions with a finite mean square error the following equations are true

$$\begin{array}{@{}rcl@{}} \left<b,M_{-x}b\right>&=&\left<\mathcal{F}^{-1}b, \mathcal{F}^{-1}\left(M_{-x}b\right)\right>=\left<\mathcal{F}^{-1}b, \mathcal{F}\left(M_{x}\hat{\hat{b}}\right)\right>\\ &=& \left<\mathcal{F}^{-1}b, T_{x}(\mathcal{F}^{-1}b)\right>\end{array} $$

where $T_{x}(\mathcal {F}^{-1}b)(\omega )=(\mathcal {F}^{-1}b)(\omega -x)$ for all $\omega \in \mathbb {R}$ and $\mathcal {F}^{-1}$ is the inverse Fourier Transform. It can be shown that the limit of $\left <{\mathcal {F}^{-1}b, T_{x}(\mathcal {F}^{-1} b)}\right >$ as x approaches ∞ is zero. Thus

$$ \lim\limits_{x\to\infty}f(x)=0. $$

(26)

By using the Cauchy-Schwarz inequality and the fact that $\|{\mathcal {F}^{-1}b}\|_{2}=\|{T_{x}(\mathcal {F}^{-1}b)}\|_{2}$ we are able to estimate as follows

$$\begin{array}{@{}rcl@{}} |f(x)|\!&=&\!|\left<\mathcal{F}^{-1}b, T_{x}(\mathcal{F}^{-1}b)\right>|\!\leq\! \|\mathcal{F}^{-1}b\|_{2}\!\cdot\!\|T_{x}(\mathcal{F}^{-1}b)\|_{2}\\ &=& \|\mathcal{F}^{-1}b\|_{2}^{2}=\left<\mathcal{F}^{-1}b,T_{0}(\mathcal{F}^{-1}b)\right>=f(0). \end{array} $$

Therefore the function f takes its maximal abolute value at x=CD.

Now we only have to show the unambiguity of this maximum value. Thus we assume that there is another value $\tilde x\in \mathbb {R}$ with $\tilde x\neq 0$ and $|{f(\tilde x)}| = |{f(0)}|$. Due to the inequality above we get that

$$|{\left<{\mathcal{F}^{-1}b,T_{\tilde x}(\mathcal{F}^{-1}b)}\right>}|=\|{\mathcal{F}^{-1}b}\|_{2}\cdot\|{T_{\tilde x}(\mathcal{F}^{-1}b)}\|_{2}.$$

However, this equality is true if and only if $T_{\tilde x}(\mathcal {F}^{-1}b)$ and $\mathcal {F}^{-1}b$ are linear dependent. Hence for all $\omega \in \mathbb {R}$

$$T_{\tilde x}(\mathcal{F}^{-1}b)(\omega)=\alpha \cdot \mathcal{F}^{-1}b(\omega)$$

where $\alpha \in \mathbb {C}$ with |α|=1. Furthermore, we can show by a complete induction that for all $k\in \mathbb {N}$

$$ T_{k\tilde x}(\mathcal{F}^{-1}b)(\omega)=\alpha^{k} \cdot \mathcal{F}^{-1}b(\omega). $$

(27)

Since |α|=1 and with Eq. (27)

$$\begin{array}{rcl} |f(k\tilde x)|&=&|\left<\mathcal{F}^{-1}b,T_{k\tilde x}(\mathcal{F}^{-1}b)\right>| = |\left<\mathcal{F}^{-1}b,\alpha^{k}\mathcal{F}^{-1}b\right>|\\ &=&|\alpha|^{k}\|\mathcal{F}^{-1}b\|_{2}^{2} = \|\mathcal{F}^{-1}b\|_{2}^{2} =f(0). \end{array}$$

Consequently,

$$\lim_{\tilde x\to \infty}|f(k\tilde x)|=f(0).$$

This is inconsistent with Eq. (26) as f(0)≠0 (otherwise a≡0). Thus the maximal absolute value is unambiguous.

□

In accordance with Eq. (7) the function $\mathcal {I}_{+}$ is defined by

$$\begin{array}{@{}rcl@{}} \mathcal{I}_{+}(t)&=& {\int}_{0}^{\infty}\frac{a(\omega)}{2}\cdot e^{-2\pi i(\omega\text{CD} +\text{CP})}\cdot e^{2\pi i\omega t} \:\mathrm{d}\omega \\ &=& {\int}_{0}^{\infty}\frac{a(\omega)}{2}e^{-2\pi i \text{CP}}\cdot e^{2\pi i\omega (t-\text{CD})} \:\mathrm{d}\omega \\ &=& {\int}_{-\infty}^{\infty}c(\omega)\cdot e^{2\pi i\omega (t-\text{CD})} \:\mathrm{d}\omega \end{array} $$

(28)

with $c(w)=\frac {a(\omega )}{2}\cdot e^{-2\pi i \text {CP}}\cdot \chi _{[0,\infty )}(\omega )$. This function c satisfies the condition from Theorem 1, so that the maximal absolute value of the function $\mathcal {I}_{+}$ is at the time t = CD.

Appendix C: $\mathcal {I}_{+}$ is an analytic signal

A complex signal z with z(t)=x(t)+i⋅y(t) is an analytic signal if and only if the imaginary part is equal to the Hilbert Transform of the real part of the complex signal, i.e.

$$ y(t) = x(t) \ast \frac{1}{\pi t}=\lim\limits_{\varepsilon \downarrow 0 }{\int}_{|{\tau}|>\varepsilon}\frac{x(t-\tau)}{\pi\tau}\:\mathrm{d}\tau=: \mathcal{H}(x(t)). $$

(29)

In Eq. (28) we have seen that the function $\mathcal {I}_{+}$ can be written as

$$\begin{array}{@{}rcl@{}} \mathcal{I}_{+}(t)&=& {\int}_{0}^{\infty}\frac{a(\omega)}{2}e^{-2\pi i \text{CP}}\cdot e^{2\pi i\omega (t-\text{CD})} \:\mathrm{d}\omega \\ &=& {\int}_{-\infty}^{\infty}G(\omega)\cdot e^{2\pi i\omega t} \:\mathrm{d}\omega \end{array} $$

(30)

where $G(w)=\frac {a(\omega )}{2}\cdot e^{-2\pi i (\omega \cdot \text {CD}+\text {CP})}\cdot \chi _{[0,\infty )}(\omega )$. Due to the Fourier inversion theorem we know that on the one hand $\widehat {\mathcal {I}_{+}}=G$ and on the other hand that for alle $\omega \in \mathbb {R}$

$$\begin{array}{@{}rcl@{}} G(\omega)=\widehat{\text{bTC}}(\omega)\cdot H(\omega) \end{array} $$

(31)

where H is the Heaviside distribution with H(ω)=1 for ω∈[0,∞) and H(ω)=0 in the other case. By applying the convolution theorem to expression (31) it turns out that $\mathcal {I}_{+}$ is the convolution of the broadband ITD tuning function and the inverse Fourier Transform of the Heaviside distribution, i.e.

$$ \mathcal{I}_{+}(t)= (\text{bTC}\ast(\mathcal{F}^{-1}H))(t). $$

(32)

The inverse Fourier Transform of the Heaviside distribution is given by the following expression $\mathcal {F}^{-1}H = \frac {1}{2}\delta _{0}+\varphi$ where δ ₀ is the dirac delta function and φ is the Cauchy principal value which is defined as

$$\displaystyle\varphi(g)=\lim\limits_{\varepsilon \downarrow 0 }{\int}_{|{x}|>\varepsilon}\frac{g(x)}{x}\:\mathrm{d}x.$$

One can show that

$$\begin{array}{@{}rcl@{}} \mathcal{I}_{+}(t)&=& (\text{bTC}\ast(\mathcal{F}^{-1}H))(t)\\ &=& \left(\text{bTC}\ast\left(\frac{1}{2}\delta_{0}+\frac{i}{2\pi}\varphi\right)\right)(t)\\ &=& \frac{\text{bTC}(t)}{2} + i\cdot\lim\limits_{\varepsilon \downarrow 0 }{\int}_{|{x}|>\varepsilon}\frac{\text{bTC}(t-x)}{2\pi x}\:\mathrm{d}x\\ &=& \frac{\text{bTC}(t)}{2}+i\cdot\mathcal{H}\left(\frac{\text{bTC}(t)}{2}\right). \end{array} $$

Consequently, the function $\mathcal {I}_{+}$ is an analytic signal.