# Wideband DOA Estimation through Projection Matrix Interpolation

J. Selva

###### Abstract

This paper presents a method to reduce the complexity of the deterministic
maximum likelihood (DML) estimator in the wideband direction-of-arrival (WDOA)
problem, which is based on interpolating the array projection matrix in the
temporal frequency variable. It is shown that an accurate interpolator like
Chebyshev’s is able to produce DML cost functions comprising just a few
narrowband-like summands. Actually, the number of such summands is far smaller
(roughly by factor ten in the numerical examples) than the corresponding
number in the ML cost function that is derived by dividing the spectrum into
separate bins. The paper also presents two spin-offs of the interpolation
method. The first is a fast procedure to compute one-dimensional search
estimators like Multiple Signal Classification (MUSIC), that exploits the
close relation between Chebyshev interpolation and the Discrete Cosine
Transform (DCT). And the second is a detection-estimation procedure for the
DML estimator. The methods in the paper are assessed in several numerical
examples.

##  I Introduction

The estimation of the directions of arrival (DOAs) to an array of sensors is a
fundamental problem in signal processing, which has been studied in multitude
papers during the last decades. Most of the existing estimators are based on
the so-called narrowband assumption, i.e, on considering that the signals’
bandwidths are much smaller than their corresponding central frequencies. As
is well known, this assumption greatly simplifies the modeling, because the
geometry of a given impinging wave just translates into a single complex
factor for each sensor-wave pair. Actually, the various estimation methods in
array processing can be viewed as effective ways to exploit this single-factor
structure, [1]. In practice, however, the signals’ bandwidths are often wide
relative relative to their central frequencies and, as a consequence, the
array geometry affects separate spectral components differently. Thus, in this
wideband case there is a factor for each combination of sensor, wave and
spectral component, and this fact complicates the estimation significantly.

In the literature, most techniques attempt to reduce this wideband case to the
narrowband one, by dividing the spectrum into separate bins, and then using a
narrowband model in each of them. Afterward, the main issue is how these
separate models should be combined. There are fundamentally two approaches for
this combination. It can either be coherent [2, 3, 4, 5] in the sense that it
incorporates the array snapshots from separate bins into a single covariance
matrix, or incoherent, meaning that a narrowband estimate is computed for each
bin first, and then these last estimates are combined through a method like a
weighted average, [6]. Besides, there exist methods that share features of the
coherent and incoherent approaches like the test of orthogonality of projected
subspaces (TOPS) in [7, 8]. The methods that employ this reduction to the
narrowband case have been surveyed in several references, [7], [9, Ch. 3] and
[10, Sec. 4.3].

We must additionally mention that the direct estimation through well-
established methodologies like the maximum likelihood (ML) principle has also
been attempted in a number of references. The deterministic and stochastic ML
estimators have been analyzed in [11] theoretically. An approximate ML
estimator has been proposed in [12, 13] and the ML estimator for near-field
sources has been analyzed in [14]. The main drawback of these direct
approaches is their complexity, given that the cost function is the sum of the
narrowband cost functions for each spectral bin, and the number of bins is
usually high.

In this paper, we propose to interpolate the projection matrix appearing in
the DML cost function using a Chebyshev interpolator, in order to reduce the
cost function’s complexity. Specifically, we will show that such interpolator
is able to exploit the smooth variation with the temporal frequency, in order
to yield a cost function composed of just a few narrowband-like components.
Though other interpolators are usable, we consider only the Chebyshev
interpolator for simplicity and due to its excellent performance.

The paper has been organized as follows. In the next sub-sections we introduce
the notation and recall several results on Chebyshev interpolation. Then, in
Sec. II we will introduce the signal model for the wideband DOA problem, and
in Sec. III the proposed interpolation method for the DML cost function. In
Sec. IV we will derive an interpolator for one-dimensional search estimators,
that will be employed in incoherent MUSIC (IC-MUSIC) and in the detection-
estimation procedure in Sec. V. Finally, we will assess the methods in the
paper in Sec. VI numerically.

###  I-A Notation

The notation will be the following.

  * •

We will write vectors and matrices in lower and upper bold face respectively.
So $y$ will denote a vector and $Y$ a matrix.

  * •

$I$ will denote the identity matrix.

  * •

$[\mbox{\boldmath$y$\unboldmath}]_{k}$ and
$[\mbox{\boldmath$Y$\unboldmath}]_{p,q}$ will represent the $k$th and $(p,q)$
elements of $y$ and $Y$ respectively.

  * •

The vector formed by appending a scalar $a$ to a column vector $y$ will be
denoted $[\mbox{\boldmath$y$\unboldmath};a]$.

  * •

’$\odot$’ will be the element-by-element product of two equal-size matrices or
vectors.

  * •

$\mbox{\boldmath$Y$\unboldmath}^{\dagger}$ will denote the pseudo-inverse of
matrix $Y$. As is well known,
$\mbox{\boldmath$Y$\unboldmath}^{\dagger}=(\mbox{\boldmath$Y$\unboldmath}^{H}\mbox{\boldmath$Y$\unboldmath})^{-1}\mbox{\boldmath$Y$\unboldmath}^{H}$
if $Y$ has full-column rank.

  * •

New symbols and functions will be introduced using the operator ’$\equiv$’.

###  I-B Background on Chebyshev interpolation

In this paper, the well-known Chebyshev interpolator will play a key role,
[15, Sec. 6.2]. In order to recall it, we first consider a generic function
${g}(x)$, defined in an interval $[a,b]$, and introduce the following
notation:

  * •

${T}_{P}(y)$ and ${U}_{P}(y)$ respectively denote the Chebyshev polynomials of
the first and second kinds and order $P$.

  * •

$y_{p}$ is the $p$th root of ${T}_{P}(y)$, $y_{p}\equiv-\cos(\pi(p-1/2)/P)$,
$p=1,\,2,\ldots,\,P$.

  * •

${y}(x)$ denotes the linear function mapping the interval $[a,b]$ onto
$[-1,1]$.

  * •

$x_{p}$ is the set of abscissas following $y_{p}={y}(x_{p})$.

  * •

${\delta}_{p}$ denotes the discrete delta function: ${\delta}_{0}=1$ and
${\delta}_{p}=0$ if $p\neq 0$.

Using these definitions, the Chebyshev interpolator of ${g}(x)$ is

$\tilde{}{g}(x)\equiv\sum_{k=0}^{P-1}c_{k}{T}_{k}({y}(x)),$ (1)

where the coefficients $c_{k}$ are given by

$c_{k}\equiv\frac{2-{\delta}_{k}}{P}\sum_{p=1}^{P}{g}(x_{P-p+1})\cos\Big{(}\frac{\pi
k}{P}\Big{(}p-\frac{1}{2}\Big{)}\Big{)}.$ (2)

Alternatively, the following formula expresses $\tilde{}{g}(x)$ in terms of
the samples ${g}(x_{p})$,

$\tilde{}{g}(x)=\sum_{p=1}^{P}{\phi}_{p}(x){g}(x_{p}),$ (3)

where

${\phi}_{p}(x)\equiv\frac{{T}_{P}(y(x))}{P{U}_{P-1}(y_{p})({y}(x)-y_{p})}.$
(4)

A key observation is that (2) is the type-2 DCT of the samples ${g}(x_{p})$,
except for a scale factor and an index reversal. This fact allows us to
oversample $\tilde{}{g}(x)$ using a procedure akin to the zero-padding FFT
method [16], but based on the DCT. Specifically, consider an integer ${R>P}$,
the roots of the $R$th Chebyshev polynomial
${y}^{\prime}_{r}\equiv-\cos(\pi(r-1/2)/R)$, ($r=1,\ldots,\,R$), and the
corresponding abscissas for $\tilde{}{g}(x)$, $x^{\prime}_{r}$,
[$y^{\prime}_{r}\equiv{y}(x^{\prime}_{r})$]. In order to compute the samples
$\tilde{}{g}(x^{\prime}_{r})$, $r=1,\,2,\ldots,R$, we may first compute the
coefficients $c_{k}$ through one type-2 DCT [Eq. (2)], then add zero padding
up to length $R$, and finally obtain the desired samples
$\tilde{}{g}(x^{\prime}_{r})$ through one type-3 DCT (inverse of type-2 DCT).

##  II Signal model

Consider a linear array of $M$ sensors at positions $d_{m}$, into which $K$
wideband signals impinge with angles of arrival $\theta_{k}$ relative to the
broadside, ($m=1,\,\ldots,\,M$, $k=1,\,\ldots,\,K$). We assume that these
signals are passband with spectra contained in $[f_{o}-B/2,f_{o}+B/2]$, where
$B>0$ and $f_{o}>B/2$. After demodulation from frequency $f_{o}$ to baseband,
the signal from the $m$th sensor follows the model

${x}_{m}(t)=\sum_{k=1}^{K}{e}^{-j2\pi
f_{o}\tau_{m}\gamma_{k}}{b}_{m}(\gamma_{k}){s}_{k}(t-\tau_{m}\gamma_{k})+{w}_{m}(t),$
(5)

where

  * •

$\tau_{m}$ is the delay between the array’s reference point and the $m$th
sensor at the propagation velocity,

$\tau_{m}\equiv\frac{d_{m}}{c},\;\textrm{(propagation velocity $c$)},$ (6)

  * •

$\gamma_{k}$ is the sine of the $k$th angle of arrival,

$\gamma_{k}\equiv\sin\theta_{k},$ (7)

  * •

${b}_{m}(\gamma_{k})$ is the $m$th sensor pattern,

  * •

${s}_{k}(t)$ is the lowpass equivalent of the $k$th impinging signal,

  * •

and the ${w}_{m}(t)$ are independent complex white noise processes of equal
variance.

Next, we assume that the receiver takes $N$ samples of each ${x}_{m}(t)$ with
a period $T$ fulfilling the Nyquist condition ($BT<1$), and then computes the
DFT. If these last samples are denoted $\tilde{}{x}(r/(NT))$ (integer $r$),
then for large $N$ they follow the model

$\begin{array}[]{l}\displaystyle\tilde{}{x}_{m}\Big{(}\frac{r}{NT}\Big{)}=\sum_{k=1}^{K}{e}^{-j2\pi(f_{o}+r/(NT))\tau_{m}\gamma_{k}}{b}_{m}(\gamma_{k})\tilde{}{s}_{k}\Big{(}\frac{r}{NT}\Big{)}\\\
\displaystyle\hfill+\tilde{}{w}_{m}\Big{(}\frac{r}{NT}\Big{)},\end{array}$ (8)

with $r=r_{1},\,r_{1}+1,\ldots,\,r_{2}$, where $r_{1}$ and $r_{2}$ are the
first and last frequency indices in which at least one signal has significant
power, and $\tilde{}{s}_{k}(r/(NT))$ and $\tilde{}{w}_{m}(r/(NT))$ have the
obvious definitions. In vector notation, this model can be written as

$\mbox{\boldmath$x$\unboldmath}_{r}=\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$s$\unboldmath}_{r}+\mbox{\boldmath$w$\unboldmath}_{r},$
(9)

where

$\begin{array}[]{l}\displaystyle[\mbox{\boldmath$x$\unboldmath}_{r}]_{m}\equiv\tilde{}{x}_{m}\Big{(}\frac{r}{NT}\Big{)},\,\,[\mbox{\boldmath$\gamma$\unboldmath}]_{k}\equiv\gamma_{k},\\\
\displaystyle[\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})]_{m,k}\equiv{e}^{-j2\pi(f_{o}+r/(NT))\tau_{m}\gamma_{k}}{b}_{m}(\gamma_{k}),\\\
\displaystyle[\mbox{\boldmath$s$\unboldmath}_{r}]_{k}\equiv\tilde{}{s}_{k}\Big{(}\frac{r}{NT}\Big{)},\,\,[\mbox{\boldmath$S$\unboldmath}]_{\cdot,r}\equiv\mbox{\boldmath$s$\unboldmath}_{r},\\\
\displaystyle[\mbox{\boldmath$w$\unboldmath}_{r}]_{m}\equiv\tilde{}{w}_{m}\Big{(}\frac{r}{NT}\Big{)}.\end{array}$
(10)

For the model in (9), the DML estimator of $S$ and $\gamma$ is given by the
arguments minimizing the cost function

${L}_{0}(\mbox{\boldmath$S$\unboldmath},\mbox{\boldmath$\gamma$\unboldmath})\equiv\sum_{r=r_{1}}^{r_{2}}\|\mbox{\boldmath$x$\unboldmath}_{r}-\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$s$\unboldmath}_{r}\|^{2}.$
(11)

As is well known, this cost function can be minimized in
$\mbox{\boldmath$s$\unboldmath}_{r}$ for fixed $\gamma$, using the pseudo-
inverse of
$\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})$. If
the resulting $\mbox{\boldmath$s$\unboldmath}_{r}$ is substituted into (11),
then we obtain the compressed cost function,

${L}_{1}(\mbox{\boldmath$\gamma$\unboldmath})\equiv\sum_{r=r_{1}}^{r_{2}}\|\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$x$\unboldmath}_{r}\|^{2},$
(12)

where
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$
denotes the projection matrix for the orthogonal complement of
$\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})$,

$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})\equiv\mbox{\boldmath$I$\unboldmath}-\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})^{\dagger}.$
(13)

At this point, we have that the number of summands in (12) is usually very
high (on the thousands), and each of them depends on a different projection
matrix
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$.
This implies that the direct minimization of
${L}_{1}(\mbox{\boldmath$\gamma$\unboldmath})$ would be too expensive
computationally. In the next section, we propose an interpolation method that
reduces the number of summands drastically.

##  III Proposed interpolation method

In the DML cost function in (12), the number of summands can be very large
though, additionally, the projection matrix
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$
varies smoothly with the frequency index $r$. We may exploit this last fact by
means of an interpolator for
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$
in order to reduce the number of summands. In this paper, we propose to employ
the Chebyshev interpolator in (3) due to its high accuracy for a small number
of input function values. Specifically, we propose to interpolate
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$
using the Chebyshev interpolator in (3) with linear mapping

${y}(x)=\frac{2x-r_{1}-r_{2}}{r_{2}-r_{1}},$ (14)

and interpolation abscissas

$\rho_{p}\equiv\frac{r_{2}-r_{1}}{2}x_{p}+\frac{r_{2}+r_{1}}{2},$ (15)

where ${x}_{p}$ was already defined in Sec. I-B. The resulting interpolator is

$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})\approx\sum_{p=1}^{P}\phi_{p}(r)\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath}).$
(16)

Let us insert this formula into (12) in order to obtain a new DML cost
function. We have

${L}_{1}(\mbox{\boldmath$\gamma$\unboldmath})\approx{L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath})\equiv\sum_{r=r_{1}}^{r_{2}}\Big{\|}\sum_{p=1}^{P}{\phi}_{p}(r)\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$x$\unboldmath}_{r}\Big{\|}^{2}\\\
=\sum_{r=r_{1}}^{r_{2}}{\mathrm{tr}}\Big{\\{}\sum_{p=1}^{P}{\phi}_{p}(r)\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$x$\unboldmath}_{r}\mbox{\boldmath$x$\unboldmath}_{r}^{H}\Big{\\}}\\\
=\sum_{p=1}^{P}{\mathrm{tr}}\Big{\\{}\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath})\sum_{r=r_{1}}^{r_{2}}{\phi}_{p}(r)\mbox{\boldmath$x$\unboldmath}_{r}\mbox{\boldmath$x$\unboldmath}_{r}^{H}\Big{\\}}.$
(17)

Finally, we define the correlation matrices

$\mbox{\boldmath$R$\unboldmath}_{p}\equiv\sum_{r=r_{1}}^{r_{2}}{\phi}_{p}(r)\mbox{\boldmath$x$\unboldmath}_{r}\mbox{\boldmath$x$\unboldmath}_{r}^{H},$
(18)

in order to write ${L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath})$ as

${L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath})=\sum_{p=1}^{P}{\mathrm{tr}}\\{\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$R$\unboldmath}_{p}\\}.$
(19)

This is the DML cost function proposed in this paper. For assessing later its
performance, let us introduce another cost function based on the usual
approach, in which the spectrum is divided into separate bins using a filter
bank. For this, divide the range $[r_{1},r_{2}]$ into $P_{b}$ bins of the form

$I_{b,p}\equiv[r_{1}+(p-1)d,r_{1}+pd[,\;1\leq p\leq P_{b},$ (20)

where $d\equiv(r_{2}-r_{1})/P_{b}$, let $\rho_{b,p}$ denote the $p$th bin’s
center, and consider the bin (or piecewise constant) interpolator given by

$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})\approx\sum_{p=1}^{P_{b}}\phi_{b,p}(r)\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{b,p},\mbox{\boldmath$\gamma$\unboldmath}),$
(21)

where $\phi_{b,p}(r)$ is the indicator function for the $p$th bin,

$\phi_{b,p}(r)\equiv\left\\{\begin{array}[]{ll}1&\>\>\mbox{if}\>\>\rho_{b,p}-d/2\leq
r<\rho_{b,p}+d/2\\\ 0&\text{otherwise.}\end{array}\right.$ (22)

Now, a derivation similar to that in (17) but using the bin interpolator in
(21) yields the cost function

${L}_{1}(\mbox{\boldmath$\gamma$\unboldmath})\approx{L}_{1,b}(\mbox{\boldmath$\gamma$\unboldmath})\equiv\sum_{p=1}^{P_{b}}{\mathrm{tr}}\\{\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{b,p},\mbox{\boldmath$\gamma$\unboldmath})\mbox{\boldmath$R$\unboldmath}_{b,p}\\},$
(23)

where $\mbox{\boldmath$R$\unboldmath}_{b,p}$ is the covariance matrix for the
$p$th bin,

$\mbox{\boldmath$R$\unboldmath}_{b,p}\equiv\sum_{r\in
I_{b,p}}\mbox{\boldmath$x$\unboldmath}_{r}\mbox{\boldmath$x$\unboldmath}_{r}^{H}.$
(24)

We will check in the numerical examples that, though (19) and (23) have a
similar form, the Chebyshev interpolator has produced a much smaller number of
summands; i.e, $P$ in (19) is much smaller than $P_{b}$ in (23).

Finally, note that in (19) each summand can be interpreted as a narrowband
cost function at the frequency $f_{o}+\rho_{p}/(NT)$. However, the expected
values of the covariance matrices $\mbox{\boldmath$R$\unboldmath}_{p}$ do not
necessarily admit a decomposition into a noise and signal subspace, given that
(18) is a wideband average of the outer products
$\mbox{\boldmath$x$\unboldmath}_{r}\mbox{\boldmath$x$\unboldmath}_{r}^{H}$;
(i.e, $r$ varies along the full frequency range $[r_{1},r_{2}]$).

##  IV Complexity reduction applied to one-dimensional search estimators:
beamformer and IC-MUSIC

For estimators involving a one-dimensional search in the DOA parameter
$\gamma$, like IC-MUSIC, an additional Chebyshev interpolator may
substantially reduce the total computational burden, by decreasing the number
of cost-function evaluations required. We may illustrate this technique for a
simple beamformer estimator, and then adapt it to other cases.

In a so-called beamformer estimator, the estimates are the abscissas of the
$K$ main local minima of ${L}_{1,c}(\gamma)$ [or ${L}_{1,b}(\gamma)$], where
$\gamma$ is a scalar parameter. The technique that we propose starts by
introducing a Chebyshev interpolator for
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\gamma)$
of the form in (3) with weight functions ${\eta}_{q}(\gamma)$ [akin to
$\phi_{p}(x)$], domain $[-1,1]$, and interpolation abscissas
$\tilde{\gamma}_{q}$. For a sufficiently large order $Q$, this interpolator is

$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\gamma)\approx\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\tilde{\gamma}_{q}).$
(25)

Next, we write
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\tilde{\gamma}_{q})$
as

$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\tilde{\gamma}_{q})=\mbox{\boldmath$I$\unboldmath}-\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}^{H},$
(26)

where $\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}$ is the normalized signature

$\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}\equiv\frac{\mbox{\boldmath$a$\unboldmath}(\rho_{p},\tilde{\gamma}_{q})}{\|\mbox{\boldmath$a$\unboldmath}(\rho_{p},\tilde{\gamma}_{q})\|}.$
(27)

Finally, we insert (26) and (27) into (19) to obtain an interpolation formula
for ${L}_{1,c}(\gamma)$,

${L}_{1,c}(\gamma)\approx\sum_{p=1}^{P}{\mathrm{tr}}\Big{\\{}\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\tilde{\gamma}_{q})\mbox{\boldmath$R$\unboldmath}_{p}\Big{\\}}\\\
=\sum_{p=1}^{P}{\mathrm{tr}}\Big{\\{}\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\big{(}\mbox{\boldmath$I$\unboldmath}-\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}^{H}\big{)}\mbox{\boldmath$R$\unboldmath}_{p}\Big{\\}}\\\
=\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\sum_{p=1}^{P}{\mathrm{tr}}\\{\mbox{\boldmath$R$\unboldmath}_{p}\\}-\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}^{H}\mbox{\boldmath$R$\unboldmath}_{p}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}.$
(28)

In this expression, the $p$-index sum is the value of
${L}_{1,c}(\tilde{\gamma}_{q})$, and it is only necessary to compute $Q$ such
values for fully specifying the formula. This fact allows us to compute the
beamformer estimates in the following steps:

  1. 1.

Compute ${L}_{1,c}(\tilde{\gamma}_{q})$, $q=1,\,2\ldots,\,Q$, [$p$-index sum
in (28)]. Note that the vectors $\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}$
are constant and, therefore, can be pre-computed.

  2. 2.

Upsample the previous step’s output to a number of samples $R>Q$, using the
DCT zero-padding technique in Sec. I-B.

  3. 3.

Locate the main local minima in the sample sequence from the previous step.

  4. 4.

Refine, if necessary, the local minima’s abscissas using a one-dimensional
optimization method like Newton’s. This steps involves no evaluations of
${L}_{1,c}(\gamma)$, given that we may instead evaluate its interpolator in
(28).

As can be readily inferred, this technique is applicable to any other one-
dimensional search estimator, and we will employ in this paper two of them.
The first is an extended beamformer that takes into account a fixed component
vector $\mbox{\boldmath$\gamma$\unboldmath}_{o}$. Its cost function is
${L}_{1,c}([\mbox{\boldmath$\gamma$\unboldmath}_{o};\gamma])$ and its
interpolator is obtained by replacing $\mbox{\boldmath$R$\unboldmath}_{p}$
with
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath}_{o})\mbox{\boldmath$R$\unboldmath}_{p}$
in (28),

${L}_{1,c}([\mbox{\boldmath$\gamma$\unboldmath}_{o},\gamma])\approx\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\\\
\cdot\sum_{p=1}^{P}{\mathrm{tr}}\\{\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath}_{o})\mbox{\boldmath$R$\unboldmath}_{p}\\}-\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}^{H}\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath}_{o})\mbox{\boldmath$R$\unboldmath}_{p}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}.$
(29)

The second estimator is IC-MUSIC, in which the one-dimensional search is
performed on the total pseudo-spectrum consisting of the sum of the usual
MUSIC pseudo-spectra in separate bins. This sum is the result of replacing
$\mbox{\boldmath$R$\unboldmath}_{b,p}$ with
$\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}$
in the cost function in (23), where $\mbox{\boldmath$U$\unboldmath}_{b,p}$
spans the signal subspace of $\mbox{\boldmath$R$\unboldmath}_{b,p}$,

${L}_{1,mu}(\gamma)\equiv\sum_{p=1}^{P_{b}}{\mathrm{tr}}\\{\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{b,p},\gamma)\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}\\}.$
(30)

Next, in order to obtain its interpolator, note that this expression has the
same form as ${L}_{1,c}(\gamma)$ in (19), but with
$\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}$
in place of $\mbox{\boldmath$R$\unboldmath}_{p}$. Therefore, if we repeat the
derivation in (28) but starting from (30), we obtain the last line of (28) but
with
$\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}$
in place of $\mbox{\boldmath$R$\unboldmath}_{p}$,

${L}_{1,mu}(\gamma)\\\
\approx\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\cdot\sum_{p=1}^{P}{\mathrm{tr}}\\{\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}\\}-\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}^{H}\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}$
(31)

Finally, we may simplify this last expression as follows:

${L}_{1,mu}(\gamma)\approx\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\Big{(}KP-\sum_{p=1}^{P}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}^{H}\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}\Big{)}\\\
\approx
KP-\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\sum_{p=1}^{P}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}^{H}\mbox{\boldmath$U$\unboldmath}_{b,p}\mbox{\boldmath$U$\unboldmath}_{b,p}^{H}\hat{}\mbox{\boldmath$a$\unboldmath}_{p,q}.$
(32)

This is the interpolator proposed for the IC-MUSIC pseudo-spectrum in (30). In
(32), we have used the approximation $\sum_{q=1}^{Q}{\eta}_{q}(\gamma)\approx
1$, given that the Chebyshev interpolator also approximates the constant
functions.

##  V Proposed DML detection-estimation method

In the sequel, we present a detection-estimation method for computing the DML
estimate, that resembles well-known procedures for the narrowband case, like
those in [17, 18] and [19, Sec. 4.6)]. Fundamentally, it is a combination of
two steps, one for detecting additional components, and another for refining a
given estimate. The first is a variant of the Generalized Likelihood Ratio
(GLR) test in [20], and the second is the modified variable projection method
(MVP) in [21] which converges in a small number of iterations.

The method operates on a given estimate
$\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$, repeating in turn a
detection step followed by an estimation step, until a statistical test fails.
The sub-indices $K$ and $\alpha$ in
$\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$ are the length of this same
vector and the number of iterations in the estimation step respectively. The
initial vector $\mbox{\boldmath$\gamma$\unboldmath}_{0,0}$ is empty. Both
steps are described in the next two sub-sections.

###  V-A Detection step

Given an iterate $\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$ (which can
be initially empty), the detection method decides whether to look for an
additional parameter $\gamma$ and, if so, selects as $\gamma$ the minimum of
${L}_{1,c}([\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha};\gamma])$. The
detection is based on a statistical test for
${L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha})$, assuming that
$\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$ is the true parameter vector.
More precisely, if $\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$ is the
true vector, then ${L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha})$
follows a $\chi^{2}$ distribution with $2(M-K)(r_{2}-r_{1}+1)$ degrees of
freedom. Thus, an additional component is sought if
${L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha})>A$, with

$A\equiv\frac{\sigma^{2}}{2}F^{-1}(1-P_{FA}),$ (33)

where $F$ is the $\chi^{2}$ cumulative distribution function, and $P_{FA}$ a
fixed false-alarm probability.

We may summarize this detection step as follows. If
${{L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha})<A}$, the iterative
process finishes and the final estimate is
$\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$, but if
${{L}_{1,c}(\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha})>A}$ then it
proceeds to the estimation step with the new vector
$\mbox{\boldmath$\gamma$\unboldmath}_{K+1,0}=[\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha};\gamma]$,
where

$\gamma=\arg\min_{\gamma^{\prime}}{L}_{1,c}([\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha};\gamma^{\prime}]).$
(34)

This minimum is located using the DCT zero-padding technique in Sec. IV.

###  V-B Estimation step

This step is an implementation of the MVP method in [21], already used in [18,
22] and [19, Sec. 4.6)] for the narrowband problem. For introducing it, let us
define first the following matrix of differentials with the same size as
$\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})$,

$[\mbox{\boldmath$D$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})]_{m,k}\equiv\frac{d}{d\gamma_{k}}[\mbox{\boldmath$A$\unboldmath}(r,\mbox{\boldmath$\gamma$\unboldmath})]_{m,k}.$
(35)

Its explicit expression can be readily computed from (10). Also, it is
convenient to define the following shorthand notation

$\displaystyle\mbox{\boldmath$A$\unboldmath}_{p,K,\alpha}$
$\displaystyle\equiv\mbox{\boldmath$A$\unboldmath}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}),\;\mbox{\boldmath$D$\unboldmath}_{p,K,\alpha}\equiv\mbox{\boldmath$D$\unboldmath}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}),$
(36)
$\displaystyle\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp,p,K,\alpha}$
$\displaystyle\equiv\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(\rho_{p},\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}).$
(37)

Given an iterate $\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$, the MVP
method refines it using the iteration

$\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha+1}=\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}-\mu\mbox{\boldmath$H$\unboldmath}_{K,\alpha}^{-1}\mbox{\boldmath$g$\unboldmath}_{K,\alpha},$
(38)

where usually ${\mu=1}$, though ${\mu<1}$ may be used to ensure a cost
function decrease, and $\mbox{\boldmath$g$\unboldmath}_{K,\alpha}$ and
$\mbox{\boldmath$H$\unboldmath}_{K,\alpha}$ are the gradient and approximate
Hessian of (19) at
$\mbox{\boldmath$\gamma$\unboldmath}=\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$,

$\displaystyle\mbox{\boldmath$g$\unboldmath}_{K,\alpha}\equiv-2\mathrm{Re}\\{\mathrm{diag}\\{\sum_{p=1}^{P}\mbox{\boldmath$A$\unboldmath}_{p,{K,\alpha}}^{\dagger}\mbox{\boldmath$R$\unboldmath}_{p}\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp,p,{K,\alpha}}\mbox{\boldmath$D$\unboldmath}_{p,{K,\alpha}}\\}\\}$
(39) $\displaystyle\mbox{\boldmath$H$\unboldmath}_{K,\alpha}\equiv
2\mathrm{Re}\\{\sum_{p=1}^{P}(\mbox{\boldmath$D$\unboldmath}_{p,{K,\alpha}}^{H}\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp,p,{K,\alpha}}\mbox{\boldmath$D$\unboldmath}_{p,{K,\alpha}}\\})^{T}$
$\displaystyle{}\hfill\odot(\mbox{\boldmath$A$\unboldmath}_{p,{K,\alpha}}^{\dagger}\mbox{\boldmath$R$\unboldmath}_{p}(\mbox{\boldmath$A$\unboldmath}_{p,{K,\alpha}}^{\dagger})^{H})\\}.$
(40)

These expressions can be efficiently computed from the Housholder QR
decompositions of the matrices
$\mbox{\boldmath$A$\unboldmath}_{p,{K,\alpha}}$, as shown in [19, Sec.
4.6.4b)]. Actually, it is not necessary to compute neither
$\mbox{\boldmath$A$\unboldmath}_{p,{K,\alpha}}^{\dagger}$ nor
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp,p,{K,\alpha}}$. The
iteration in (38) is repeated until there is no significant reduction in the
cost function’s value. Afterward, the execution proceeds to the detection step
again.

##  VI Numerical examples

We have performed several numerical examples following the signal model in
Sec. II, which are presented in the sequel. In them, the main parameters were
the following:

Central frequency. The signals’ central frequency was $f_{o}=2.4$ GHz.

Received signals. There were two simulation scenarios, that we term
“independent-signals” (IS) and “correlated-signals” (CS) scenarios. In the IS
scenario, the lowpass equivalents of the received signals were three linearly-
modulated signals ${s}_{o,k}(t)$ with raised-cosine modulating pulse (roll-off
0.2). The modulation sequences in the three were variance-one independent and
complex white. The corresponding lowpass signals at the sensor array reference
point were $a^{\prime}_{k}{s}_{o,k}(t-\tau^{\prime}_{k})$, with

Amplitudes $a^{\prime}_{k}$: |  $0.626+j0.7798$, $-0.4432-j0.552$,
---|---
|  $0.3138+j0.3908$,
Delays $\tau^{\prime}_{k}/(2f_{o})$: |  $0$, $0.6$, $37.53$.

In the CS scenario, the lowpass equivalents were generated as in the IS
scenario, except for the fact that the three initial signals ${s}_{o,k}(t)$
were the same one.

DFT length. ${N=2048}$ and the sampling period $T$ followed ${BT=0.8}$. The
frequency index range was $[r_{1},r_{2}]=[-819,818]$, where index 0 was
exactly placed at frequency $f_{o}$.

Sensor array. Uniform linear array with $M=10$ sensors and isotropic patterns,
[${b}_{m}(\gamma)=1$ for all sensors].

Angles of arrival. $-0.7895$, $-0.6816$, and $0.2734$ rads. Corresponding
parameters in $\gamma$: $-0.71$, $-0.63$, and 0.27.

Signal-to-noise ratio (SNR). Ratio of the signal and noise powers in (8),
averaged over all the values of $r$ and $m$ in that equation.

Detection false-alarm probability. In the example about the detection method,
we have used $P_{FA}=0.01$.

Estimators. We have evaluated three estimators:

  * •

_IC-MUSIC._ Incoherent MUSIC with interpolated pseudo-spectrum in (32). The
local minima were located using the DCT technique in Sec. I-B, followed by a
one-dimensional Newton method for each minimum.

  * •

_BinML._ DML estimator using the bin cost function in (23).

  * •

_ChebML._ Proposed DML implementation based on minimizing the cost function in
(19).

In both BinML and ChebML, the detection-estimation method in Sec. V was used
in all cases. However, for known number of signals $K$, the pair of detection
and estimation steps was executed exactly $K$ times, i.e, the test in Sec. V-A
was not employed.

One-dimensional search in BinML and ChebML. In all cases, it was $Q=50$ and
the Chebyshev grid was oversampled by factor $2$ using the DCT zero-padding
method in Sec. I-B.

Number of Monte Carlo trials. The estimators’ performances have been evaluated
in 100 Monte Carlo trials.

###  VI-A Projection matrix interpolation performance

Fig. 1 shows the error of the bin interpolator in (21) when used to
approximate the projection matrix
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$
for ${P_{b}=47}$ bins. This last number of bins ensures an error below $-50$
dB at all frequencies, as can be readily checked, and is the minimum $P_{b}$
fulfilling this condition. In this figure, the error is defined as the maximum
amplitude among the elements of
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})-\hat{}\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$,
where
$\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$
and
$\hat{}\mbox{\boldmath$P$\unboldmath}_{{}\scriptscriptstyle\perp}(r,\mbox{\boldmath$\gamma$\unboldmath})$
are the true and interpolated projection matrices respectively. The
oscillations in this figure are due to the fact that the error is negligible
close to the bins’ centers and increases away from them. Fig. 2 shows the same
error but for the Chebyshev interpolator and several interpolation orders $P$.
Note that ${P=4}$ ensures by large the $-50$ dB threshold while, as just
commented, the bin interpolator requires $P_{b}=47$. As can be readily seen,
the Chebyshev interpolation error decreases very fast with $P$. Actually,
$P=6$ already places the error below the $-100$ dB threshold. These error
curves have a $-\infty$ value at the positions corresponding to the Chebyshev
polynomial roots after the linear interval mapping
$[-1,1]\rightarrow[r_{1},r_{2}]$, as explained in Sec. (I-B).

Figure 1:  Interpolation error for the projection matrix using the bin
interpolator with 47 bins. This is the smallest number of bins ensuring the
$-50$ dB threshold for all frequencies. 234567 Figure 2:  Projection matrix
approximation error using the Chebyshev interpolator versus the frequency
index. 57911 Figure 3:  Projection matrix approximation error using the
Chebyshev interpolator versus the parameter difference $\gamma_{2}-\gamma_{1}$
for several interpolation orders $P$.

Fig. 3 shows the same error for the Chebyshev interpolator but versus the
parameter difference $|\gamma_{2}-\gamma_{1}|$, assuming two impinging waves.
Again, we can see the strong reduction in the error as $P$ increases. Also,
note that there is no relevant error increase for small parameter differences.

###  VI-B Performance of the BinML and ChebML estimators

Fig. 4(a) shows the RMS error of the IC-MUSIC estimator for $P=10$ and $50$
bins, and of the BinML estimator for $P=40$ to $70$ bins, in the IS scenario
and assuming a known number of signals $K$. Note that IC-MUSIC if far from the
Cramer-Rao (CR) bound and BinML attains it.

Fig. 4(b) shows the same error but for the ChebML estimator. Note that the
mismatch is negligible already for $P=5$. Thus, in this example ChebML
requires at least ten times less cost function summands than BinML, with the
corresponding savings in computational burden.

501040, 50, 6070 – – – | IC-MUSIC
---|---
——– | BinML
\+ + + | CR bound
(a) RMS error in dB for the estimation of $\gamma_{1}$ of estimators IC-MUSIC
(10 and 50 bins) and BinML (40 to 70 bins) in the IS scenario. Crosses (’+’)
mark the CR bound.

12345, 6,$\ldots$ (b) RMS error in dB for the estimation of $\gamma_{1}$ of
estimator ChebML ($P=1,\,2,\ldots$) in the IS scenario. Crosses (’+’) mark the
CR bound.

Figure 4:  RMS error in dB for the estimation of $\gamma_{1}$ of estimators
IC-MUSIC (a), BinML (a), and ChebML (b) in the IS scenario, together with CR
bound.

Figs. 5(a) and 5(b) are the equivalents of Figs. 4(a) and 4(b) but for the CS
scenario. Note that in this case the CR bound is higher due to the correlation
among impinging waves. Here, we may draw roughly the same conclusions as in
the IS case, except for the fact that BinML and ChebML require a somewhat
higher order at high SNRs, ($P=6$ or 7 for ChebML, $P_{b}=70$ or more for
BinML).

40, 5060, 705010 – – – | IC-MUSIC
---|---
——– | BinML
\+ + + | CR bound
(a) RMS error in dB for the estimation of $\gamma_{1}$ of estimators IC-MUSIC
(10 and 50 bins) and BinML (40 to 70 bins) in the CS scenario. Crosses (’+’)
mark the CR bound.

1234567 (b) RMS error in dB for the estimation of $\gamma_{1}$ of estimator
ChebML ($P=1,\,2,\ldots$) in the CS scenario. Crosses (’+’) mark the CR bound.

Figure 5:  RMS error in dB for the estimation of $\gamma_{1}$ of estimators
IC-MUSIC (a), BinML (a), and ChebML (b) in the CS scenario, together with CR
bound.

###  VI-C Performance of the detection-estimation scheme

Fig. 6 shows the performance of the detection scheme in Sec. V in the IS
scenario for ${P=6}$ (and unknown $K$). The curve in Fig. 6(b) is the
probability of detecting the three components in $\gamma$, and Fig. 6(a) shows
the RMS error but estimated for those realizations in which the three
components had been detected. The detector employed the actual noise variance
$\sigma^{2}$ in the test in (33), if the SNR was below 37 dB. For larger SNRs,
$\sigma^{2}$ is kept at the 37 dB level in the detection test, given that the
error in the approximation of the true wave parameters by
$\mbox{\boldmath$\gamma$\unboldmath}_{K,\alpha}$ becomes noticeable at high
SNRs. Note that the detector obtains the three components with probability
almost equal to one for SNRs above the $-5$ dB threshold.

(a) RMS error and CR bound (’+’) for the detection scheme in Sec. V.  (b)
Probability of detecting the three impinging waves for the detection scheme in
Sec. V.

Figure 6:  Performance of the detection scheme in Sec. V in the IS scenario.

##  VII Conclusions

We have presented an interpolation method for the wideband DOA problem, that
makes it possible to obtain DML cost functions composed of a small number of
narrowband-like components. The method is based on approximating the
projection matrix involved in the DML cost function using a Chebyshev
interpolator. Compared with the approach based on dividing the spectrum into
separate bins and then applying a narrowband model in each of them, the
proposed method produces a much smaller number of cost function summands, so
reducing the computational burden significantly. We have presented two spin-
offs of the method that combine the already commented interpolator with
another one of the same type. This second interpolator is used to reduce the
complexity of one-dimensional searches. These spin-off are an implementation
of the IC-MUSIC estimator and a detection-estimation scheme for the DML
estimator. The methods in the paper have been assessed in several numerical
examples.

## References

  * [1] Harry L. van Trees,  Detection, Estimation, and Modulation Theory. Part IV, Optimum array processing,  John Wiley & Sons, Inc, first edition, 2002.
  * [2] H. Wang and M. Kaveh,  “Coherent signal-subspace processing for the detection and estimation of angles of arrival of multiple wide-band sources,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 33, no. 4, pp. 823–831, Aug 1985.
  * [3] S. Valaee and P. Kabal,  “Wideband array processing using a two-sided correlation transformation,”  IEEE Transactions on Signal Processing, vol. 43, no. 1, pp. 160–172, Jan 1995.
  * [4] T. K. Yasar and T. E. Tuncer,  “Wideband DOA estimation for nonuniform linear arrays with Wiener array interpolation,”  in 2008 5th IEEE Sensor Array and Multichannel Signal Processing Workshop, July 2008, pp. 207–211.
  * [5] W. J. Zeng and X. L. Li,  “High-resolution multiple wideband and nonstationary source localization with unknown number of sources,”  IEEE Transactions on Signal Processing, vol. 58, no. 6, pp. 3125–3136, June 2010.
  * [6] M. Wax, Tie-Jun Shan, and T. Kailath,  “Spatio-temporal spectral analysis by eigenstructure methods,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 32, no. 4, pp. 817–827, Aug 1984.
  * [7] Yeo-Sun Yoon, L. M. Kaplan, and J. H. McClellan,  “TOPS: new DOA estimator for wideband signals,”  IEEE Transactions on Signal Processing, vol. 54, no. 6, pp. 1977–1989, June 2006.
  * [8] A. K. Shaw,  “Improved wideband DOA estimation using modified TOPS (mTOPS) algorithm,”  IEEE Signal Processing Letters, vol. 23, no. 12, pp. 1697–1701, Dec 2016.
  * [9] Sathish Chandran, Ed.,  Advances in Direction-of-Arrival Estimation,  Artech House, 2006.
  * [10] Engin Tuncer and Benjamin Friedlander, Eds.,  Classical and modern direction-of-arrival estimation,  Elsevier, 2009.
  * [11] M. A. Doron, A. J. Weiss, and H. Messer,  “Maximum-likelihood direction finding of wide-band sources,”  IEEE Transactions on Signal Processing, vol. 41, no. 1, pp. 411–414, Jan 1993.
  * [12] Lean Yip, Joe C. Chen, Ralph E. Hudson, and Kung Yao,  “Cramer-Rao bound analysis of wideband source localization and DOA estimation,”  in International Symposium on Optical Science and Technology. International Society for Optics and Photonics, 2002, pp. 304–316.
  * [13] L. Yip, C. E. Chen, R. E. Hudson, and K. Yao,  “DOA estimation method for wideband color signals based on least-squares joint approximate diagonalization,”  Proceedings of Sensor Array and Multichannel Signal Processing, pp. 104–107, 2008.
  * [14] Joe C. Chen, Ralph E. Hudson, and Kung Yao,  “Maximum-likelihood source localization and unknown sensor location estimation for wideband signals in the near-field,”  IEEE Transactions on Signal Processing, vol. 50, no. 8, pp. 1843–1854, 2002.
  * [15] John C. Mason and David C. Handscomb,  Chebyshev polynomials,  CRC Press, 2002.
  * [16] Donald Fraser,  “Interpolation by the FFT revisited–an experimental investigation,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 37, no. 5, pp. 665–675, May 1989.
  * [17] I. Ziskind and M. Wax,  “Maximum likelihood localization of multiple sources by alternating projection,”  IEEE Transactions on Acoustics, Speech, and Signal Processing, vol. 36, pp. 1553–1560, Oct. 1988.
  * [18] J. Selva,  “ML estimation and detection of multiple frequencies through periodogram estimate refinement,”  IEEE Signal Processing Letters, vol. 24, no. 3, pp. 249–253, March 2017.
  * [19] T. V. Ho, J. G. McWhirter, A. Nehorai, U. Nickel, B. Ottersten, B. D. Steinberg, P. Stoica, M. Viberg, Z. Zhu, S. Haykin, et al.,  Radar array processing, vol. 25,  Springer Science & Business Media, 2013.
  * [20] Mats Viberg, Björn Ottersten, and T. Kailath,  “Detection and estimation in sensor arrays using weighted subspace fitting,”  IEEE Transactions on Signal Processing, vol. 39, no. 11, pp. 2436–2449, Nov. 1991.
  * [21] L. Kaufman,  “A variable projection method for solving separable nonlinear least squares problems,”  BIT, vol. 15, pp. 49–57, 1975.
  * [22] J. Selva,  “An efficient Newton-type method for the computation of ML estimators in a Uniform Linear Array,”  IEEE Transactions on Signal Processing, vol. 53, no. 6, pp. 2036–2045, June 2005.


1Nagoya University, Japan, 2TARVO, Inc., Japan

# Electrolaryngeal Speech Intelligibility Enhancement Through Robust
Linguistic Encoders

###### Abstract

We propose a novel framework for electrolaryngeal speech intelligibility
enhancement through the use of robust linguistic encoders. Pretraining and
fine-tuning approaches have proven to work well in this task, but in most
cases, various mismatches, such as the speech type mismatch (electrolaryngeal
vs. typical) or a speaker mismatch between the datasets used in each stage,
can deteriorate the conversion performance of this framework. To resolve this
issue, we propose a linguistic encoder robust enough to project both EL and
typical speech in the same latent space, while still being able to extract
accurate linguistic information, creating a unified representation to reduce
the speech type mismatch. Furthermore, we introduce HuBERT output features to
the proposed framework for reducing the speaker mismatch, making it possible
to effectively use a large-scale parallel dataset during pretraining. We show
that compared to the conventional framework using mel-spectrogram input and
output features, using the proposed framework enables the model to synthesize
more intelligible and naturally sounding speech, as shown by a significant 16%
improvement in character error rate and 0.83 improvement in naturalness score.

Index Terms—  Intelligibility enhancement, electrolaryngeal speech, atypical
speech

##  1 Introduction

Voice conversion (VC) [1], the task known as changing the speaker information
while keeping linguistic information unchanged, has had rapid improvements in
the age of deep learning. One of its sub-applications, intelligibility
enhancement [2, 3, 4], has made way for atypical speakers to regain the
ability to speak like typical speakers. Atypical speakers have difficulties in
producing phoneme sounds and speak at a slower rate, making daily
communication a tedious task for them. One type of atypical speech,
electrolaryngeal (EL) speech, is produced by speakers diagnosed with a
disrupted larynx, the organ responsible for generating the source excitation.
While an electrolarynx [5] is used as a replacement for the larynx, the
resulting speech becomes unnatural due to the electrolarynx producing robotic-
like source excitation and being unable to produce natural pitch variation.
For pitch-based languages like Japanese, changing the pitch throughout a
sentence as well as the use of voiced/unvoiced sounds is essential to infer
the meaning of different words, making this an important task.

Several previous works in intelligibility enhancement have found that an
effective solution is to first learn the alignments between typical and
atypical speech through a parallel dataset. Since EL speakers speak at a
slower rate and are unable to pronounce some phoneme sounds, learning the
alignment between the two is important in this task. For example, [6] does
this by using a strong sequence model such as a Transformer [7, 8]. Due to
data scarcity and the data-hungry nature of Transformer-based models, several
works [4, 9] have emphasized the effectiveness of pretraining on a large-scale
typical speech dataset and fine-tuning it on a small-scale atypical speech
dataset. However, a major problem in this naive pretraining and fine-tuning
framework is that the typical and EL speech types are vastly different from
each other. Thus, although a simple pretraining and fine-tuning approach
brings improvements, there is a performance ceiling in such an approach. Our
previous work [10] resolved this by observing that fine-tuning first with
large-scale synthetic speech can effectively soften the mismatch between the
speech types and speakers, making the pretraining and fine-tuning approach
more effective. However, there is still a lot of room for improvement in
further reducing the speech type and speaker mismatches, as the synthesis
performance is still far from human-level speech.

We resolve the speech type and speaker mismatch issues encountered in
pretraining and fine-tuning approaches by introducing a new framework which
uses recognition, alignment, and synthesis modules. Specifically, we use
strong recognition modules containing dense linguistic information (such as
phonetic posteriorgrams [11] and HuBERT [12, 13] features) as input and output
features of the alignment module, effectively allowing the alignment module to
focus on solely learning linguistic features. Through the recognition module
allowing focus on modeling linguistic information, we effectively remove the
speech type and speaker mismatches occuring between each stage during
pretraining and fine-tuning, resulting in better performance compared to the
baseline. Moreover, with the use of a Diffusion-based [14] synthesis decoder
to generate the target speaker mel-spectrogram from the HuBERT output
features, we shift the burden of synthesizing to a target speaker to this
module due to its strong generation capabilities, improving the generation
quality of the waveform. Finally, this proposed framework optimizes the use of
parallel VC pretraining to further improve performance. Our contributions are
as follows:

  * •

We propose a novel framework for electrolaryngeal speech intelligibility
enhancement, composed of recognition, alignment, and synthesis modules. We
show that using this framework can synthesize speech with a 16% CER
improvement and a 0.83 higher naturalness score compared to the baseline.

  * •

We resolve the speech type mismatch issues by developing a linguistic encoder
robust to both EL and typical speech types. Through a unified representation
being used as inputs, the alignment module can focus on solely modeling the
linguistic features, resulting in significantly more intelligible speech.

  * •

We show the other important components of the framework, such as using HuBERT
output features and using parallel VC pretraining in ablation studies.

##  2 Conventional framework

We use our previous work [10] as our baseline, which uses the Transformer [7,
8] to transform the mel-spectrogram of an EL speech utterance into a mel-
spectrogram of a typical speech utterance. A pretraining technique using text-
to-speech (TTS) and autoencoder (AE) was used to efficiently learn linguistic
information from a large-scale typical speech data. The process is done by
training a TTS model with the target speaker. Then, an AE-style pretraining is
conducted by using the decoder of the TTS model as initialization parameters,
and reconstructing the target speaker by also using it as inputs. The decoder
parameters are frozen, such that the encoder is efficiently pretrained. The
network is first fine-tuned on the parallel synthetic EL and typical speech.
We found that fine-tuning first on synthetic EL speech (even with lots of
mispronunciations in synthesis) softens the speech type and speaker mismatches
when fine-tuning. Then, network is fine-tuned on the target EL and typical
speech data. Moreover, since the TTS pretraining technique uses text
information as inputs and models strong linguistic information [8], such a
speaker-independent pretraining style was beneficial in reducing the speech
type and speaker mismatches when fine-tuning.

Although bringing in large improvements, the framework is still limited as it
is still far from human-level speech. One main problem is that the use of mel-
spectrograms contains a lot of information related to the speech type and
speaker, which degrades the performance due to the speech type and speaker
mismatches between the datasets used in fine-tuning in each stage. One way to
resolve this issue is by using linguistic encoders, which has shown success in
several works in speech synthesis [1, 11]. By using a linguistic encoder to
extract dense linguistic information from speech and using these as the input
and output features, the focus during conversion can be on the linguistic-
related features, reducing the speech type and speaker mismatches. This has
been applied to intelligibility enhancement where works such as [15] use a
fine-tuned automatic speech recognition (ASR) model on the atypical speech;
however, this approach does not use the pretraining and fine-tuning framework
used in the majority of works. Thus, although the ASR model can effectively
extract linguistic features from the atypical speech, the pretraining on the
large-scale typical speech dataset becomes less effective as the ASR model
fine-tuned on the atypical speech cannot properly decode typical speech. We
further investigate how to develop such a robust linguistic encoder and its
performance when used as input and output features to improve intelligibility.

##  3 Proposed framework

An overview of the entire framework can be seen in Fig. 1. We detail the task
of each module of the proposed framework below.

Fig. 1: An overview of the proposed framework. The framework contains three
main modules to convert from EL to typical speech: the recognition, alignment,
and synthesis modules. Note that each module is trained separately.

###  3.1 Recognition module

The recognition module uses a linguistic encoder estimating the phonetic
posteriorgrams (PPGs) from the bottleneck features of an ASR encoder to
extract the linguistic information. In our previous work [16], we showed that
an effective approach to improving ASR model performance for EL speech is
through a three-stage training framework. First, the model was pretrained on a
large-scale typical speech dataset. Next, we fine-tuned the network on
synthesized EL speech in an intermediate fine-tuning stage. Due to the limited
data in training an EL speech synthesis model, the synthesized speech also
contained lots of mispronunciations. However, we found that since the model
used this stage to learn the EL speech characteristics instead of the
linguistic information, it was sufficient enough for the synthetic EL speech
to only represent the EL speech features. Finally, we fine-tuned the network
on the ground truth EL speech to learn the linguistic features and decode at a
high accuracy. We adopt this framework as the backbone of the recognition
module.

The goal of the linguistic encoder now is to be robust enough to remove the
speech type features from both typical and EL speech, while also accurately
extracting linguistic information. With a unified representation, the
performance of a pretraining and fine-tuning framework becomes robust to the
speech type mismatches. Although this has been an easy task in typical VC,
several ASR works have found developing speaker-independent models [17, 18,
19, 20] for atypical speakers a difficult task due to the high variance in
their speech. A naive approach to resolve this would be to simply fine-tune
the ASR model on both the EL and typical speech such that the model is not
only optimized for EL speech. However, similar to previous works, we found
that fine-tuning the ASR model on both types of speech at the same time causes
degradations.

To improve this, we simply introduce a speech type ID loss $L_{\text{SID}}$
during training. Since our previous work discovered that the intermediate
fine-tuning focuses on learning speech type identity features, we make the
network learn both speech types during this stage. Let
$X=\\{X_{\text{TYP}},X_{\text{EL}}\\}$ be the training data, which is composed
of a typical and an EL dataset $X_{\text{TYP}}$ and $X_{\text{EL}}$. The
speech type ID loss $L_{\text{SID}}$ identifies whether the speaker is a
typical or EL speaker and is optimized using a binary cross-entropy loss.
Since we use both EL and typical data, we mask the outputs from the typical
speech inputs during the calculation of the CTC/Attention losses
$L_{\text{ctc}}$ and $L_{\text{attn}}$ [21]. The masking avoids making the
model learn two highly variant types of speech, improving decoding
performance. Through this approach, we effectively optimize the ASR model for
EL speech, while also ensuring that it does not forget how to decode typical
speech. We show in Eq. 1 the detailed loss calculation during intermediate
fine-tuning. After the intermediate fine-tuning stage, we fine-tune on
$X_{\text{EL}}$ with the CTC/Attention losses $L_{\text{ctc}}$ and
$L_{\text{attn}}$ [21] as usual.

$L_{\text{ASR}}=L_{\text{SID}}(X)+L_{\text{ctc}}(X_{EL})+L_{\text{attn}}(X_{EL})$
(1)

###  3.2 Alignment module

The alignment module resolves the intelligibility enhancement aspect. To
improve intelligibility, the alignment module needs to fulfill two tasks.
First, due to the different temporal structure of EL speech, the model needs
to increase the speaking rate similar to a typical speaker. Next, since EL
speakers cannot produce certain phonemes, the alignment module also needs to
correct the phoneme pronunciation. Similar to the baseline described in
Section 2, we adopt the use of a Transformer [7, 8] sequence-to-sequence model
to resolve these issues. We also follow the same fine-tuning procedure with
synthetic data and then the target data due to its success.

We improve this framework by using the PPG features produced by the
recognition module as the inputs. These PPG features would further reduce the
mismatches in speech type and speakers during pretraining and fine-tuning, as
the linguistic encoder allows the alignment module to solely focus on modeling
linguistic information. To further reduce the burden on the alignment network,
we also use HuBERT as the output features. Aside from HuBERT providing dense
linguistic information, using a variant of HuBERT with soft features [13] has
also been found successful in removing speaker features and in cross-lingual
settings, which would further improve the performance of the synthesis module
described later.

Moreover, the TTS/AE pretraining described in Section 2 was initially used in
our baseline due to the unavailability of a large-scale parallel dataset;
however, with the release of [22], we first verify whether parallel VC is
indeed better. Although using parallel VC pretraining would directly model the
fine-tuning task, this would not contain the speaker-independent properties of
the TTS/AE pretraining, which might cause more degradations in the multiple
fine-tuning stages due to the speech type and speaker mismatches. However,
owing to the proposed framework focusing solely on linguistic features, we
remove this possibility.

###  3.3 Synthesis module

As our goal is to force the alignment module to focus only on modeling
linguistic information, the task of synthesizing into the target speaker is
placed on a synthesis module. Since the typical dataset used as a target
speaker is also limited in size, we use a Diffusion model [14] as the decoder
of this module, as this framework has been proven effective in synthesizing
speech in a target speaker even in few-shot settings [23]. To improve the few-
shot performance of the Diffusion model, similar to [23], we pretrain first on
a large-scale multi-speaker dataset with classifier-free guidance [24] and use
fixed speaker embeddings and HuBERT features as conditioning features. Then,
we adapt the model to the few-shot data for another set of iterations. To
train the model, we iteratively add noise for $N$ timesteps to the mel-
spectrogram and predict the noise at timestep $n$ during training by using the
noisy mel-spectrogram at $n-1$ as input along with the conditioning features.
During inference, we pass in Gaussian noise and predict the mel-spectrogram
after $N$ iterations. Finally, to synthesize the audio waveforms from the
predicted mel-spectrograms, we use HiFiGAN (V1) [25] as the vocoder.

##  4 Experimental Settings

###  4.1 Datasets

Table 1: Total duration (in mins) and number of utterances for each split used in our experiment. | Minutes | No. of Utterances
---|---|---
Dataset | Train | Dev | Test | Train | Dev | Test
TYPICAL | 4.38 | 2.14 | 2.17 | 116 | 40 | 40
EL | 5.77 | 2.96 | 2.99 | 116 | 40 | 40

The EL dataset is spoken in Japanese; thus, unless otherwise stated, the
described dataset is also in Japanese. For the recognition module, we followed
the same training framework as in [16] to train a linguistic encoder. We first
pretrained on a large-scale typical speech dataset containing around 2k hours
of speech data [26]. Next, we fine-tuned the network on a total of 27k
utterances of synthetic EL data and typical speech. Finally, we fine-tuned the
network on our privately acquired EL speech data. We evaluated the performance
of the linguistic encoder using the aforementioned EL data, and its parallel
counterpart spoken by a typical speaker. We used a 116/40/40 split for the
train, dev, and test data. We also conducted ablation studies on the
performance when using a larger dataset of 15 speakers, containing both their
simulated EL (by using an external electrolarynx) and typical speech.

For the alignment module, we first pretrained our network on HiFiCaptain [22],
a large-scale parallel dataset of typical speakers, containing around 18k
utterances in total. We used the female speaker as the source and the male
speaker as the target. Then, we used the same setup as in [10] where we first
fine-tuned on synthetic EL, synthesized from text from the JSUT [27] corpus.
We then fine-tuned the model on our target parallel EL and typical speech data
in Table 1. Note that compared to our previous work in [10], this current
split is different, as we composed the evaluation data with longer utterances
to show the effectiveness of the proposed method. We used the same 116/40/40
split used in the recognition module, so the test data is unseen by both the
recognition and alignment modules.

For the synthesis decoder, we used the JVS dataset [28], a dataset containing
30 hours of speech from 100 speakers, to pretrain the model before fine-tuning
it on the target typical speech. For the synthesis vocoder, we used a
pretrained model on VCTK [29], an English dataset with 44 speakers of around
40 hours in total.

Table 2: Resulting CER on both EL and typical speech with different training data setups. Training data and method | EL | TYPICAL
---|---|---
Speaker-dependent on typical speech | 77.1 | 4.3
Speaker-dependent on target EL speech [16] | 15.9 | 61.5
Fine-tuned on multiple EL and typical | 28.7 | 3.8
Fine-tuned on EL and typical speech | 18.1 | 16.6
Proposed fine-tuned on EL and typical with speech type ID loss only |  16.2 |  6.8
Proposed fine-tuned on EL and typical with speech type ID loss and masking |  13.8 |  6.0
Table 3: Objective and subjective evaluation results on the synthesized speech from different systems, along with the ground truth recorded speech. MOS is calculated with a 95% confidence interval. We detail the input and output features, along with the pretraining method used in the alignment module. (System) Description | Inputs | Outputs | Pretraining method | MCD ($\downarrow$) | CER ($\downarrow$) | F0 RMSE ($\downarrow$) | F0 CORR ($\uparrow$) | MOS ($\uparrow$)
---|---|---|---|---|---|---|---|---
(1) Baseline [10] | mel | mel | TTS/AE | 7.78 | 35.0 | 51.19 | 0.30 |  2.42 $\pm$ 0.17
(2) Baseline (ablation) | mel | mel | Parallel VC | 7.70 | 33.5 | 49.95 | 0.35 |  2.38 $\pm$ 0.18
(3) Proposed | PPG | HuBERT | TTS/AE | 7.45 | 32.2 | 50.39 | 0.28 |  2.90 $\pm$ 0.17
(4) Proposed | PPG | HuBERT | Parallel VC | 7.14 | 19.0 | 52.41 | 0.29 | 3.25 $\pm$ 0.15
(5) Proposed (ablation) | PPG | mel | Parallel VC | 7.54 | 29.1 | 51.16 | 0.37 |  2.78 $\pm$ 0.18
Ground truth | - | - | - | - | 4.3 | - | - |  4.85 $\pm$ 0.07

###  4.2 Model architecture

The recognition model used a Conformer architecture [30], which has 12 layers
for both the encoder and decoder, the same as in [16]. On the other hand, the
Transformer model in the alignment module has six layers for both the encoder
and decoder, the same as in [10]. The diffusion model of the synthesis decoder
was adopted from [31] and uses 512-dimension channels to predict the noise
between each timestep. To handle speech inputs, we replaced the text encoder
with the PPG encoder as the conditioning features. To integrate speaker
information, we used a pretrained WavLM
model111https://huggingface.co/microsoft/wavlm-base-plus-sv (which was fine-
tuned for speaker verification) as speaker embeddings and fused it to each
residual block using conditional layer normalization [32]. We set the number
of diffusion steps $N$ to 100. No changes were made in HiFiGAN (V1).

###  4.3 Evaluation metrics

For objective evaluations, we measured the synthesis quality through metrics
such as character error rate (CER), mel-cepstral distortion (MCD), log F0 root
mean square error (F0 RMSE), and log F0 correlation (F0 CORR). For CER, we
used the same Conformer model in Table 2 trained on the large-scale typical
speech. For subjective evaluations, we recruited 15 native Japanese speakers
to measure the naturalness of the synthesized speech using a 5-scale mean
opinion score (MOS) test222Demo: lesterphillip.github.io/icassp2024_el_sie.

##  5 Results and Discussion

###  5.1 Validating the recognition module

We first present how to develop a robust linguistic encoder. We investigate
different training setups as shown in Table 2. First, we see the difficulty in
using a speaker-independent model, as optimizing on either EL or typical
speech results in degradations in the other. Using a model optimized just on
EL speech would degrade the large-scale pretraining stage. Next, we see that
fine-tuning with multiple EL and typical speakers to make the model more
generalized, also does not have any effectiveness, caused by the high variance
between these speakers. Finally, we show that simply fine-tuning the model on
both EL and typical speech can be effective but not fully optimized, as there
is still a performance gap from the speaker-dependent setups.

We show that our proposed method of using a speech type ID loss and masking
the typical speech during CTC-Attention loss calculation makes the model learn
to decode both EL and typical speech. This is because the model learns how to
decode EL speech, while also not forgetting typical speech features learned
during pretraining through the speech type ID loss. To verify this, removing
the masking of typical speech results in slightly worse scores. Through this,
we can decode both EL and typical speech at an accuracy similar to the
speaker-dependent setups.

###  5.2 Comparison of input/output features

We show the effectiveness of the proposed linguistic encoder in this task. As
seen in Table 3 our proposed method of using the PPG/HuBERT features (Sys. 4)
can significantly improve the synthesized speech with a 16% improvement in CER
and 0.83 in naturalness score over Sys. 1, the baseline that uses mel-
spectrograms as inputs. This proves our initial hypothesis that the proposed
linguistic encoder can effectively remove speech type information while also
extracting accurate linguistic information. We also conducted a study by using
mel-spectrogram outputs. As shown in Sys. 5, using HuBERT instead of mel-
spectrograms as outputs helps further stabilizes the model, as it also
contains dense linguistic information similar to the PPGs. Aside from this,
Sys. 3 and 4 that used HuBERT features and the synthesis decoder had the top
naturalness scores, further showing the effectiveness and necessity of a
synthesis decoder over directly predicting the mel-spectrogram.

###  5.3 Comparison of pretraining techniques

In Section 3.2, we discussed that the TTS/AE pretraining also helps in
resolving the speech type and speaker mismatches during pretraining and fine-
tuning through its speaker-independent pretraining style [8]. However, upon
comparing the baseline techniques, Sys. 2 has slightly better scores than Sys.
1 except in MOS. Thus, we prove that TTS/AE is not sufficient to create a
speaker-independent property. Through the proposed approach in Sys. 4, we can
directly model the fine-tuning task by using parallel VC pretraining, while
also being able to implement a speaker-independent property by using
PPG/HuBERT as input and output features, which reduces the mismatches during
each fine-tuning stage. It is important to note that although Sys. 3 used both
speaker-independent training styles, since the input (PPG) and output (HuBERT)
features were different, the proposed method was not able to fully utilize the
effectiveness of AE pretraining. Moreover, we find that compared to the other
systems, Sys. 4 has the highest F0 RMSE score and the second lowest F0 CORR
score, showing that the proposed method truly allowed the alignment module to
focus on modeling linguistic information, but caused a small tradeoff in
modeling pitch.

##  6 Conclusions

We proposed the use of robust linguistic encoders that removes speech features
from both EL and typical speech. The major benefit that this brings is that it
creates a unified representation for both EL and typical speech, reducing the
speech type mismatches between each dataset in a pretraining and fine-tuning
framework. The proposed method allows the model to focus on modeling
intelligibility, where it outperforms the baseline with a 16% improvement in
CER and a 0.83 higher naturalness score.

Acknowledgements This work was partly supported by AMED under Grant Number
JP21dk0310114, Japan, and by JST CREST under Grant Number JPMJCR19A3.

##  7 References

## References

  * [1] Berrak Sisman et al. “An Overview of Voice Conversion and Its Challenges: From Statistical Modeling to Deep Learning” In _IEEE/ACM TASLP_ 29, 2021, pp. 132–157
  * [2] Robert L. MacDonald et al. “Disordered Speech Data Collection: Lessons Learned at 1 Million Utterances from Project Euphonia” In _Proc. Interspeech_ , 2021, pp. 4833–4837 DOI: 10.21437/Interspeech.2021-697
  * [3] Fadi Biadsy et al. “Parrotron: An End-to-End Speech-to-Speech Conversion Model and its Applications to Hearing-Impaired Speech and Speech Separation” In _Proc. Interspeech_ , 2019
  * [4] Rohan Doshi et al. “Extending Parrotron: An End-to-End, Speech Conversion and Speech Recognition Model for Atypical Speech” In _Proc. ICASSP_ , 2021, pp. 6988–6992 DOI: 10.1109/ICASSP39728.2021.9414644
  * [5] Mark I. Singer and Eric D. Blom “An Endoscopic Technique for Restoration of Voice after Laryngectomy” PMID: 7458140 In _Annals of Otology, Rhinology & Laryngology_ 89.6, 1980, pp. 529–533 DOI: 10.1177/000348948008900608
  * [6] Ming-Chi Yen et al. “Mandarin Electrolaryngeal Speech Voice Conversion with Sequence-to-Sequence Modeling” In _Proc. ASRU_ , 2021, pp. 650–657 DOI: 10.1109/ASRU51503.2021.9687908
  * [7] Ashish Vaswani et al. “Attention is all you need” In _Proc. NeurIPS_ 30, 2017
  * [8] Wen-Chin Huang et al. “Pretraining Techniques for Sequence-to-Sequence Voice Conversion” In _IEEE/ACM TASLP_ 29, 2021, pp. 745–755
  * [9] Zhehuai Chen et al. “Conformer Parrotron: a Faster and Stronger End-to-end Speech Conversion and Recognition Model for Atypical Speech” In _Proc. Interspeech_ , 2021
  * [10] Ding Ma et al. “Two-stage training method for Japanese electrolaryngeal speech enhancement based on sequence-to-sequence voice conversion” In _Proc. SLT_ , 2023, pp. 949–954 IEEE
  * [11] Songxiang Liu et al. “Any-to-Many Voice Conversion With Location-Relative Sequence-to-Sequence Modeling” In _IEEE/ACM TASLP_ 29, 2021, pp. 1717–1728 DOI: 10.1109/TASLP.2021.3076867
  * [12] Wei-Ning Hsu et al. “Hubert: Self-supervised speech representation learning by masked prediction of hidden units” In _IEEE/ACM TASLP_ 29 IEEE, 2021, pp. 3451–3460
  * [13] Benjamin Niekerk et al. “A Comparison of Discrete and Soft Speech Units for Improved Voice Conversion” In _Proc. ICASSP_ , 2022
  * [14] Jonathan Ho, Ajay Jain and Pieter Abbeel “Denoising diffusion probabilistic models” In _Proc. NeurIPS_ 33, 2020, pp. 6840–6851
  * [15] Disong Wang et al. “Learning Explicit Prosody Models and Deep Speaker Embeddings for Atypical Voice Conversion” In _Proc. Interspeech 2021_ , 2021, pp. 4813–4817 DOI: 10.21437/Interspeech.2021-285
  * [16] Lester Phillip Violeta et al. “Intermediate Fine-Tuning Using Imperfect Synthetic Speech for Improving Electrolaryngeal Speech Recognition” In _Proc. ICASSP_ , 2023 DOI: 10.1109/ICASSP49357.2023.10095931
  * [17] Lester Phillip Violeta and Tomoki Toda “An Analysis of Personalized Speech Recognition System Development for the Deaf and Hard-of-Hearing” In _Proc. APSIPA_ , 2023
  * [18] Joel Shor et al. “Personalizing ASR for Dysarthric and Accented Speech with Limited Data” In _Proc. Interspeech_ , 2019, pp. 784–788
  * [19] Jordan R. Green et al. “Automatic Speech Recognition of Disordered Speech: Personalized Models Outperforming Human Listeners on Short Phrases” In _Proc. Interspeech_ , 2021, pp. 4778–4782
  * [20] Jimmy Tobin and Katrin Tomanek “Personalized Automatic Speech Recognition Trained on Small Disordered Speech Datasets” In _Proc. ICASSP_ , 2022, pp. 6637–6641 DOI: 10.1109/ICASSP43922.2022.9747516
  * [21] Suyoun Kim, Takaaki Hori and Shinji Watanabe “Joint CTC-attention based end-to-end speech recognition using multi-task learning” In _Proc. ICASSP_ , 2017, pp. 4835–4839
  * [22] Takuma Okamoto, Yoshinori Shiga and Hisashi Kawai “Hi-Fi-CAPTAIN: High-fidelity and high-capacity conversational speech synthesis corpus developed by NICT”, 2023
  * [23] Sungwon Kim, Heeseung Kim and Sungroh Yoon “Guided-TTS 2: A diffusion model for high-quality adaptive text-to-speech with untranscribed data” In _arXiv preprint arXiv:2205.15370_ , 2022
  * [24] Jonathan Ho and Tim Salimans “Classifier-free diffusion guidance” In _Proc. NeurIPS_ , 2021
  * [25] Jungil Kong, Jaehyeon Kim and Jaekyoung Bae “HiFi-GAN: Generative adversarial networks for efficient and high fidelity speech synthesis” In _Proc. NeurIPS_ , 2020, pp. 17022–17033
  * [26] Shintaro Ando and Hiromasa Fujihara “Construction of a Large-Scale Japanese ASR Corpus on TV Recordings” In _Proc. ICASSP_ , 2021, pp. 6948–6952 DOI: 10.1109/ICASSP39728.2021.9413425
  * [27] Ryosuke Sonobe, Shinnosuke Takamichi and Hiroshi Saruwatari “JSUT corpus: free large-scale Japanese speech corpus for end-to-end speech synthesis” arXiv, 2017 DOI: 10.48550/ARXIV.1711.00354
  * [28] Shinnosuke Takamichi et al. “JVS corpus: free Japanese multi-speaker voice corpus” In _arXiv preprint arXiv:1908.06248_ , 2019
  * [29] J. C. and K. “CSTR VCTK corpus: English multi-speaker corpus for CSTR voice cloning toolkit” In _University of Edinburgh. The Centre for Speech Technology Research (CSTR)_ , 2012 URL: %5Curl%7Bhttp://dx.doi.org/10.7488/ds/1994%7D
  * [30] Anmol Gulati et al. “Conformer: Convolution-augmented Transformer for Speech Recognition” In _Proc. Interspeech_ , 2020, pp. 5036–5040 DOI: 10.21437/Interspeech.2020-3015
  * [31] Jinglin Liu et al. “Diffsinger: Singing voice synthesis via shallow diffusion mechanism” In _arXiv preprint arXiv:2105.02446_ 2, 2021
  * [32] Yuzi Yan et al. “AdaSpeech 2: Adaptive Text to Speech with Untranscribed Data”, 2021 arXiv:2104.09715 [cs.SD]