A linear programming approach to the tracking of partials

TL;DR

This paper introduces a linear programming method for tracking sinusoidal chirps that outperforms classical algorithms in noisy conditions and has lower complexity than existing approaches.

Contribution

It presents a novel LP formulation for sinusoidal partial tracking that improves efficiency and accuracy over traditional greedy and Viterbi-based methods.

Findings

LP approach outperforms classical algorithm in noisy environments

Complexity of LP method is less than previous approaches

LP method effectively tracks sinusoidal chirps in high noise levels

Abstract

A new approach to the tracking of sinusoidal chirps using linear programming is proposed. It is demonstrated that the classical algorithm of McAulay and Quatieri is greedy and exhibits exponential complexity for long searches, while approaches based on the Viterbi algorithm exhibit factorial complexity. A linear programming (LP) formulation to find the best $L$ paths in a lattice is described and its complexity is shown to be less than previous approaches. Finally it is demonstrated that the new LP formulation outperforms the classical algorithm in the tracking of sinusoidal chirps in high levels of noise.

Figures

Tables

Table 1. Table 1: Parameters of q 𝑞 q th chirp. ν 0 subscript 𝜈 0 \nu_{0} and ν 1 subscript 𝜈 1 \nu_{1} are the initial and final frequency of the chirp in Hz.

$q$	${\varphi }_q$	${\omega }_q$	${\psi }_q$	${\nu }_0$	${\nu }_1$
0	0	0.20	2.45 $\times {10}^{-6}$	500	600
1	0	0.39	4.91 $\times {10}^{-6}$	1000	1200
2	0	0.59	7.36 $\times {10}^{-6}$	1500	1800

Equations

$$\rho =\left\{(i,j):\mathcal{D}({\theta }_i,{\theta }_j)\le {\Delta }_{\text{LP}},0\le i<M,0\le j<M,i\ne j\right\}$$

$${\mathbf{c}}_{\rho }=\left\{D({\theta }_i,{\theta }_j)\forall (i,j)\in \rho \right\}$$

$${\mathbf{A}}_{r_{\text{I}},\mathcal{B}(p)}^{\text{I}}=\{\begin{array}{ll}1& \text{if }{p}_1=r_{\text{I}}+N_0\\ 0& \text{otherwise}\end{array},0\le r_{\text{I}}<R_{\text{I}},p\in \rho$$

$$0\le {\mathbf{A}}^{\text{I}}\mathbf{x}\le 1$$

$${\mathbf{A}}_{r_{\text{O}},\mathcal{B}(p)}^{\text{O}}=\{\begin{array}{ll}1& \text{if }{p}_0=r_{\text{O}}\\ 0& \text{otherwise}\end{array},0\le r_{\text{O}}<R_{\text{O}},p\in \rho$$

$$0\le {\mathbf{A}}^{\text{O}}\mathbf{x}\le 1$$

$${\mathbf{A}}_{r_{\text{B}}}^{\text{B}}={\mathbf{A}}_{r_{\text{B}}}^{\text{B}}-{\mathbf{A}}_{r_{\text{B}}+N_0}^{\text{B}}\text{ for rows }0\le r_{\text{B}}<R_{\text{B}}$$

$${\mathbf{A}}^{\text{B}}\mathbf{x}=0$$

$${\mathbf{A}}_{r_{\text{C}}}^{\text{C}}=\sum _{k=a}^b{\mathbf{A}}_k^{\text{O}}$$

$${\mathbf{A}}^{\text{C}}\mathbf{x}=L1$$

$$\underset{\mathbf{x}}{\min }{\mathbf{c}}_{\rho }^T\mathbf{x}$$

$$\mathbf{G}\mathbf{x}=\left[\begin{array}{c}{\mathbf{A}}^{\text{I}}\\ {\mathbf{A}}^{\text{O}}\\ -\mathbf{I}\end{array}\right]\mathbf{x}\le \left[\begin{array}{c}1\\ 1\\ 0\end{array}\right]$$

$$\mathbf{A}\mathbf{x}=\left[\begin{array}{c}{\mathbf{A}}^{\text{B}}\\ {\mathbf{A}}_{K-2}^{\text{C}}\end{array}\right]\mathbf{x}=\left[\begin{array}{c}0\\ L\end{array}\right]$$

$$\left[\begin{array}{cc}-\mathbf{D}{\mathbf{G}}^T\mathbf{G}& {\mathbf{A}}^T\\ \mathbf{A}& 0\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right]=\left[\begin{array}{c}t{\mathbf{c}}_{\rho }+{\mathbf{A}}^T\mathbf{d}\\ 0\end{array}\right]$$

$$\begin{gathered} \frac{2}{3}(K-1)N^6+2(K-2)(K-1)N^5+ \\ 2(K-2{)}^2(K-1)N^4+\frac{2}{3}(K-2{)}^3{N}^3 \end{gathered}$$

$$s_q(n)=\exp (j({\varphi }_q+{\omega }_{q}n+\frac{1}{2}{\psi }_q{n}^2))$$

$${\mathcal{D}}_{\text{pr.}}\left({\theta }_i^k,{\theta }_j^{k+1}\right)=\left({\omega }_i^k+{\psi }_i^{k}H-{\omega }_j^{k+1}\right)$$

Full text

Index Terms—  partial tracking, linear programming, optimization, additive synthesis, atomic decomposition, regularized approximation

1 Introduction

Atomic decompositions of audio allow for the discovery of meaningful underlying structures such as musical notes [gribonval2003harmonic] or sparse representations [plumbley2010sparse]. A classical structure sought in decompositions of speech and music signals is the sum-of-sinusoids model: windowed sinusoidal atoms in the decomposition of sufficient energy and in close proximity in both time and frequency are considered as connected. The progressions of these connected atoms in time form paths or partial trajectories.

Many authors have considered the partial tracking problem, beginning with [mcaulay1986speech]. Their technique is improved upon in [lagrange2003enhanced] with the use of linear prediction to improve the plausibility of partial tracks. Rather than seeking individual paths, in [depalle1993tracking] the most plausible sequence of connections and detachments between atoms is determined via an extension of the Viterbi algorithm proposed in [wolf1989finding]. Improvements to this technique are made in [kereliuk2008improved] by incorporating the frequency slope into atom proximity evaluations.

The latter techniques seeking globally optimal sets of paths incur great computational cost due to the large number of possible solutions. For this reason, in this paper we propose a linear programming relaxation formulation of the optimal path-set problem based on an algorithm for the tracking of multiple objects in video [jiang2007linear]. It will be shown that this algorithm has favourable asymptotic complexity and performs well on the tracking of chirps in high levels of noise.

1.1 Note on notation

The atomic decompositions used in this paper consider blocks of of contiguous samples, called frames and these frames are computed every $H$ samples, $H$ being the hop-size. We will denote the $N_{k}$ sets of parameters for atoms found in the decomposition in frame $k$ as $\theta_{0}^{k},\dotsc,\theta_{N_{k}-1}^{k}$ and the $N_{k+1}$ in frame $k+1$ as $\theta_{0}^{k+1},\dotsc,\theta_{N_{k+1}-1}^{k+1}$ where $k$ and $k+1$ refer to adjacent frames.

The total number of nodes is $M=\sum_{k=0}^{K-1}N_{k}$. $\theta_{i}^{j}$ is the $i$th node of the $j$th frame, $\theta^{j}$ the set of all the nodes in the $j$th frame and $\theta_{m}$ the $m$th node out of all $M$ nodes ($0\leq m<M$).

We are interested in paths that extend across $K$ frames where each path touches only one parameter set and each parameter set is either exclusive to a single path or is not on a path.

In this paper, indexing starts at 0. If we have a vector $\boldsymbol{x}$ then $\boldsymbol{x}_{i}$ is the $i$th row or column of that vector depending on the orientation. The same notation is used for Cartesian products, e.g., if $\alpha$ and $\beta$ are sets and $A=\alpha\times\beta$ then for the pair $a\in A$ $a_{0}$ is the first item in the pair and $a_{1}$ the second.

2 A greedy method

In this section, we present the McAulay-Quatieri method of peak matching. It is conceptually simple and a set of short paths can be computed quickly, but it can be sensitive to spurious peaks and its complexity becomes unwieldly for long searches.

In [mcaulay1986speech, p. 748] the peak matching algorithm is described in a number of steps; we summarize them here in a way comparable with the linear programming formulation to be presented shortly. In that paper, the parameters of each data point are the instantaneous amplitude, phase, and frequency but here we allow for arbitrary parameter sets $\theta$. Define a distance function $\mathcal{D}\left(\theta_{i},\theta_{j}\right)$ that computes the similarity between $2$ sets of parameters. We will now consider a method that finds $L$ tuples of parameters that are closest.

We compute the cost tensor $\boldsymbol{C}=\theta^{k}\otimes_{\mathcal{D}}\mathellipsis\otimes_{\mathcal{D}}\theta^{k+K-1}$. For each $l\in\left[0\dotsc L-1\right]$, find the indices $i_{0},\dotsc,i_{K-1}$ corresponding to the shortest distance, then remove the $i_{0},\dotsc,i_{K-1}$th rows (lines of table entries) in the their respective dimensions from consideration and continue until $L$ tuples have been determined or a distance between a pair of nodes on the path exceeds some threshold $\Delta_{\text{MQ}}$. This is summarized in Algorithm 1.

This is a greedy algorithm because on every iteration the smallest cost is identified and its indices are removed from consideration. Perhaps choosing a slightly higher cost in one iteration would allow smaller costs to be chosen in successive iterations. This algorithm does not allow for that. In other terms, the algorithm does not find a set of pairs that represent a globally minimal sum of costs. Furthermore, the algorithm does not scale well: assuming equal numbers of parameter sets in all frames, the search space grows exponentially with $K$. Nevertheless, the method is simple to implement, computationally negligible when $K$ is small, and works well with a variety of audio signals such as speech [mcaulay1986speech] and music [smith1987parshl].

3 $L$ best paths through a lattice via linear programming (LP)

In this section we show how to find L paths through a lattice of $K$ frames such that the sets of nodes on each path are disjoint. The $k$th frame of the lattice contains $N_{k}$ nodes for a total of $M=\sum_{k=0}^{K-1}N_{k}$ nodes.

Similar to the McAulay-Quatieri method we define the cost $\Delta_{\text{LP}}$ as the limiting cost under which the connection between two nodes will be considered in the LP method.

The solution vector $\boldsymbol{x}$ to the linear program shall indicate the presence of a connection between a pair of nodes by having an entry equal to $1$ and otherwise have entries equal to [math]. To enumerate the set of possible connection-pairs we define

[TABLE]

The cost vector of the objective function is then

[TABLE]

and the length of $\boldsymbol{c}_{\rho}$ is $\#\rho=\#\boldsymbol{c}_{\rho}=P$, in other words, $P$ pairs of nodes. For convenience we define a bijective mapping ${\mathcal{B}:\rho\rightarrow[0,\mathellipsis,M-1]}$ giving the index in $\boldsymbol{x}$ of the pair $p\in\rho$. For the implementation considered in this paper, $\mathcal{D}(\theta_{i},\theta_{j})=\infty$ for all $i,j$ not in adjacent frames and so $P$ will be no larger than $(K-1)N^{2}$ (assuming the same number of nodes $N$ in each frame).

The total cost of the paths in the solution is then calculated through the inner product $\boldsymbol{c}_{\rho}^{T}\boldsymbol{x}$. To obtain $\boldsymbol{x}^{\ast}$ that represents $L$ disjoint paths we must place constraints on the structure of the solution. Some of the constraints presented in the following are redundant but the redundancies are kept for clarity; later we will show which constraints can be removed without changing the optimal solution $\boldsymbol{x}^{\ast}$.

All nodes in $\boldsymbol{x}^{\ast}$ will have at most one incoming connection or otherwise no connections, a constraint that can be enforced through the following linear inequality: define $\boldsymbol{A}^{\text{I}}\in\mathbb{R}^{R_{\text{I}}\times P}$ with $R_{\text{I}}=\sum_{k=1}^{K-1}N_{k}$, the number of nodes in all the frames excluding the first. We sum all the connections into the node $r_{\text{I}}+N_{0}$ represented by the respective entry in $\boldsymbol{x}$ through an inner product with the $r_{\text{I}}$th row in $\boldsymbol{A}^{\text{I}}$ and require that this sum be between [math] and $1$, i.e.,

[TABLE]

and

[TABLE]

Similarly, to constrain the number of outgoing connections into each node, we define $R_{\text{O}}=\sum_{k=0}^{K-2}N_{k}$ and $\boldsymbol{A}^{\text{O}}\in\mathbb{R}^{R_{\text{O}}\times P}$ with

[TABLE]

and

[TABLE]

To forbid breaks in the paths it is required that the number of incoming connections into a given node equal the number of outgoing connections for the $R_{\text{B}}=\sum_{k=1}^{K-2}N_{k}$ nodes potentially having both incoming and outgoing connections.

[TABLE]

and

[TABLE]

Finally we ensure that there are $L$ paths by counting the number of connections in each frame and constraining this sum to be $L$. We choose arbitrarily to count the number of outgoing connections by summing rows of $\boldsymbol{A}^{\text{O}}$ into rows of $\boldsymbol{A}^{\text{C}}\in\mathbb{R}^{(K-1)\times P}$

[TABLE]

with $a=\sum_{j=0}^{r_{\text{C}}}N_{j}$ and $b=\sum_{j=0}^{r_{\text{C}}+1}N_{j}$ and

[TABLE]

As stated above, some of these constraints are redundant and can be removed. Indeed, we have $\boldsymbol{0}\leq\boldsymbol{x}\leq\boldsymbol{1}$, therefore we will always have $\boldsymbol{A}^{\text{I}}\boldsymbol{x}\geq 0$ and $\boldsymbol{A}^{\text{O}}\boldsymbol{x}\geq 0$. Furthermore, all but the last row of (10) can be seen as constructed from linear combinations of rows of (8) and the last row of (10) so we only require $\boldsymbol{A}^{\text{C}}_{K-2}\boldsymbol{x}=L$. Finally we always have $\boldsymbol{x}\leq\boldsymbol{1}$ because of the constraint that there be a maximum of $1$ incoming and outgoing connection from each node.

The complete LP to find the $L$ best disjoint paths through a lattice described by node connections $\rho$ is then

[TABLE]

subject to

[TABLE]

[TABLE]

where $\boldsymbol{I}$ is the identity matrix. A proof that the solution $\boldsymbol{x}^{\ast}$ will have entries equal to either [math] or $1$ can be found in [parker1988discrete, p. 167].

4 Memory complexity

To simplify notation, in this section we assume there are $N$ nodes in each frame of the lattice.

Although the matrices involved in (11) are large, only a small fraction of their values are non-zero. Matrices $\boldsymbol{A}^{\text{I}},\boldsymbol{A}^{\text{O}}\in\mathbb{R}^{N(K-1)\times P}$, but each contains only $P$ non-zero entries. Furthermore $\boldsymbol{A}^{\text{B}}\in\mathbb{R}^{N(K-1)\times P}$ but contains only $2N^{2}(K-2)$ non-zero entries while $\boldsymbol{A}^{\text{C}}_{K-2}\in\mathbb{R}^{P}$ contains merely $N$. The $\boldsymbol{x}\geq\boldsymbol{0}$ constraint requires a matrix with $P$ non-zero entries. The total memory complexity including the entries in $\boldsymbol{c}_{\rho}$ and the right-hand-sides of (11) is $2N^{2}(K-2)+4P+2N(K-1)+N+1$ non-zero floating-point numbers.

5 Complexity

Here we will compare the complexity of the LP formulation of the best $L$ paths search to the greedy McAulay-Quatieri method as well a combinatorial algorithm proposed in [wolf1989finding].

Assuming the same number of nodes $N$ in each frame of the lattice, the search for the $l$th best path in the generalized McAulay-Quatieri algorithm ($0\leq l<L$) requires a search over $(N-l)^{K}$ possible paths.

The LP formulation of the $L$-best paths problem gives results equivalent to the solution to the $L$-best paths problem proposed in [wolf1989finding]. The complexity of the algorithm by Wolf in [wolf1989finding] is equivalent to the Viterbi algorithm for finding the single best path through a trellis whose $k$th frame has $\binom{N_{k}}{L}\binom{N_{k+1}}{L}L!$ connections where $N_{k}$ and $N_{k+1}$ are the number of nodes in two consecutive frames of the original lattice. Therefore, assuming a constant number $N$ of nodes in each frame, its complexity is $O((\binom{N}{L}^{2}L!)^{2}K)$.

The complexity of the algorithm presented here is polynomial in the number of variables (the size of $\boldsymbol{x}$). Assuming we use the algorithm in [karmarkar1984new] to solve the LP, our program has a complexity of $O(P^{3.5}B^{2})$ where $B$ is the number of bits used to represent each number in the input. However, this bound is conservative considering the reported complexity of modern algorithms.

For instance, the complexity of a log-barrier interior-point method is dominated by solving the system of equations

[TABLE]

some $10$s of times [boyd2004convex, p. 590]. Each iteration then takes ${\frac{2}{3}((K-1)N^{2}+(K-2)N)^{3}}$ flops (floating-point operations) to solve (12) using a standard $LU$-decomposition [golub1996matrix, p. 98]. As $\boldsymbol{D}$ is a diagonal matrix, if the number of nodes in each frame is $N$ for all frames, then $\boldsymbol{D}\boldsymbol{G}^{T}\boldsymbol{G}$ will be a block-diagonal matrix made up of $K-1$ blocks $\boldsymbol{B}_{k}\in\mathbb{R}^{N^{2}\times N^{2}}$. The system can then be solved in

[TABLE]

flops [boyd2004convex, p. 675]; this complexity is without exploiting the sparsity of $\boldsymbol{A}$ nor the structure of $\boldsymbol{B}_{k}=\boldsymbol{D}_{k}\boldsymbol{C}$ — the product of some diagonal matrix $\boldsymbol{D}_{k}$ with an unchanging symmetric matrix $\boldsymbol{C}$.

6 Partial paths on an example signal

We compare the greedy and LP based methods for peak matching on a synthetic signal. The signal is composed of $Q=3$ chirps of constant amplitude, the $q$th chirp $s$ at sample $n$ described by the equation

[TABLE]

The parameters for the $Q$ chirps are presented in Table 6.