# Sparse Digital Cancellation of Receiver Nonlinear Distortion in Carrier
Aggregation Systems

Ahmad Gomaa, Charles Chien, Ming Lei, ChihYuan Lin and Chun-Ying Ma
Mediatek Incorporation

###### Abstract

In carrier-aggregation systems, digital baseband cancelation of self-
interference generated by receiver nonlinearity requires the estimation of
several reference signals contributions. As the nonlinearity order and
frequency selectivity of the chip response increase, the number of reference
signals significantly increases rendering the estimation of their
contributions more complex. We propose a sparsity-based approach for the
selection of the reference signals to match the distortion interference using
a few reference signals. Simulation results show significant performance
improvement over prior art with the same complexity.

##  I Introduction

In carrier aggregation (CA) systems, signals are transmitted and/or received
over multiple frequency bands simultaneously. CA is adopted in the long-term
evolution (LTE) standard [1]. Due to imperfect chip and board isolations, one
or more uplink (UL) signals can then leak into the low-noise amplifiers (LNAs)
of the receive chains of downlink (DL) signals. Practical LNAs exhibit non-
linearity behaviour [2] generating harmonic distortion (HD) and inter-
modulation distortion (IMD) of the leaked UL signals that can lie at the DL
band. This problem causes significant self-interference that degrade the
performance of the victim band [1].

In [3], digital baseband interference cancellers are used, where the UL
baseband signals were utilized to regenerate the distortion and subtract it
from the received signal. Before reaching the LNA input, UL signals are shaped
by the chip frequency response. Hence, the LNA-generated distortion is a
nonlinear function of the weighted summation of several delayed versions of
the original UL signals. The summation weights represent the chip channel
response at the corresponding lags. To avoid using nonlinear cancellation
filters, the distortion signal is decomposed into the linear combination of
the UL signals raised to different powers and delayed by different lags,
called reference signals. The reference signals weights are functions of the
original chip response and can be estimated using any linear estimation
technique [3]. In [4, 5] and [6], digital cancellation of nonlinearity
distortion generated at the transmitter and receiver, respectively. In [6],
the IMD was between one modulated signal and an unmodulated tone, unlike this
paper and [3] where both inter-modulated signals are modulated. Unlike
modulated signals, unmodulated tones do not experience frequency-selective
channel.

The reference signals in [3] are static and pre-defined regardless of the chip
channel response. Instead, we propose a sparsity-based solution to dynamically
select the reference signals, where all candidate reference signals are
included in a dictionary matrix. Then, based on their auto-correlation and
cross-correlation with the observed signal, a subset of them is selected to
represent the distortion signal. The number of selected reference signals is
flexibly set based on the design constraints on complexity and power. Unlike
[3], our approach enables different subset selection for different chip
responses. The rest of the paper is organized as follows. The system model and
problem formulation are described in Sections II and III, respectively. Our
proposed sparsity-based approach is presented in Section IV. Simulation
results are provided in Section V, and the paper is concluded in Section VI.

Notations: Lower and upper case bold letters denote vectors and matrices,
respectively, and $\mathbf{0}$ denotes the all-zero vector. Also, $(\,)^{*}$,
$(\,)^{T}$ and $(\,)^{H}$ denote the complex conjugate, transpose and
conjugate transpose operations, respectively. The notation $\left|\;\right|$
denotes the absolute value.

##  II System Model

We describe the two main nonlinearity distortion models known in the
literature, namely, HD and IMD.

###  II-A Harmonic Distortion Model

In practical transceivers, the UL signal up-converted to carrier frequency
$f_{\text{tx}}$ leaks into the receiver LNA. Due to its inevitable
nonlinearity, the LNA generates the $Q$-th order harmonic of the leakage UL
signal at frequency $Q\times f_{\text{tx}}$. In FDD systems where the receiver
oscillator frequency $f_{\text{rx}}=Qf_{\text{tx}}$, the HD of the UL leakage
signal will interfere with the desired downlink (DL) signal de-sensitizing the
whole receiver chain. In CA systems, $f_{\text{tx}}$ and $f_{\text{rx}}$ are
the UL and DL carrier frequencies of two different aggregated bands, for
example, bands 12 (710MHz) and 4 (2.13 GHz), respectively, in LTE [1]. We
denote the baseband time-domain (TD) UL signal at the digital-to-analog
converter (DAC) input by $s(n)$. The UL signal leaks to the LNA input through
a finite impulse response (FIR) channel representing the chip response
$\left\\{h(k)\right\\}_{k=0}^{L-1}$ where $L$ is the channel length. The TD
baseband equivalent of the leakage signal at the LNA input is given by:

$x(n)=\sum_{k=0}^{L-1}h(k)s(n-k)$ (1)

We write the $Q$-th order HD signal at the analog-to-digital converter (ADC)
output as follows:

$p(n)=c_{0}x^{Q}(n)$ (2)

where $c_{0}$ is related to the $Q^{\text{th}}$ order input-referred intercept
point (IIP) of the LNA [2].

###  II-B Inter-Modulation Distortion Model

In carrier aggregation, the two UL signals of different frequencies,
$f_{\text{tx,1}}$ and $f_{\text{tx,2}}$ leak to the receiver LNA. As a result
of the LNA nonlinearity, these two leakage signals inter-modulate together
creating an IMD signal sitting at a new frequency
$pf_{\text{tx,1}}+qf_{\text{tx,2}}$, where $p$ and $q$ are integer nonzero
numbers. If this new frequency equals to the frequency of the downlink signal
$f_{\text{rx}}$, then the system is said to suffer from IMD. The IMD order is
given by $Q_{I}=|p|+|q|$. In LTE, bands 3 (1750MHz) and 20 (850MHz) can cause
$3^{\text{rd}}$-order IMD ($p=2,q=-1$) to band 7 (2660MHz) [1]. We write the
$Q_{I}$-th order IMD signal at the ADC output for $p,q>0$ as follows:

$\displaystyle p(n)=c_{0}x_{1}^{|p|}(n)x_{2}^{|q|}(n)$ (3)

where $x_{1}(n)$ and $x_{2}(n)$ are the TD baseband equivalent of the two UL
leakage signals seen at the LNA input, given by:

$\displaystyle
x_{1}(n)=\sum_{k=0}^{L_{1}-1}h_{1}(k)s_{1}(n-k),\;x_{2}(n)=\sum_{k=0}^{L_{2}-1}h_{2}(k)s_{2}(n-k)$
(4)

where $s_{1}(n)$ and $s_{2}(n)$ are the complex baseband TD signals at the
inputs of the DACs associated with the two UL chains. Furthermore, $h_{1}(k)$
and $h_{2}(k)$ are the FIR channels of lengths $L_{1}$ and $L_{2}$,
respectively, representing the chip responses between the two UL chains and
the LNA input. The IMD expression in (3) can be also written for $p<0$ and/or
$q<0$ with the following modification. If $p$ or $q$ is negative, we replace
the corresponding signal $x_{1}(n)$ or $x_{2}(n)$, respectively, in Equation
(3) by its complex conjugate.

##  III Problem Formulation and prior Art

###  III-A Problem Formulation

The distortion signal $p(n)$ in (2) or (3) for HD or IMD schemes,
respectively, interferes with the desired DL signal $y(n)$ of power $P_{s}$
yielding:

$r(n)=p(n)+y(n)+z(n)$ (5)

where $z(n)$ is the complex background additive Gaussian (AWG) noise of
single-sided power spectral density $N_{o}$. The cancellation algorithm
exploits the knowledge of the transmitted UL signal $s(n)$ (or $s_{1}(n)$ and
$s_{2}(n)$ in IMD cases) to constructs the distortion signal $p(n)$ and cancel
it from $r(n)$. The algorithm requires the estimation of the FIR channel
$h(k)$ (or $h_{1}(k)$ and $h_{2}(k)$) representing the chip response to
construct the leakage signal $x(n)$ (or $x_{1}(n)$ and $x_{2}(n)$) seen at the
LNA input. However, the observed signal $r(n)$ is not linear in the unknown
channel $h(k)$ due to the LNA nonlinearity in (2) and (3). Hence, we need to
expand the polynomials of the distortion signal, e.g.,
$\left(\sum_{k=0}^{L-1}h(k)s(n-k)\right)^{Q}$ and
$\left(\sum_{k=0}^{L_{1}-1}h_{1}(k)s_{1}(n-k)\right)^{|p|}\left(\sum_{k=0}^{L_{2}-1}h_{2}(k)s_{2}(n-k)\right)^{|q|}$,
using the following multinomial theorem:

$\left(\sum_{k=1}^{K}a_{k}\right)^{N}=\sum_{\begin{subarray}{c}t_{1}+t_{2}+\ldots+t_{N}=N,\\\
t_{1},t_{2},..,t_{N}=0,1,2,..\end{subarray}}\frac{N!}{t_{1}!t_{2}!\ldots
t_{N}!}\prod_{1\leq k\leq N}a_{k}^{t_{k}}$ (6)

For instance, the first polynomial in the right-hand side (RHS) of (2) is
expanded for $Q=3$ and $L=2$ as follows:

$\displaystyle
c_{0}x^{Q}(n)=c_{0}\left(\sum_{k=0}^{1}h(k)s(n-k)\right)^{3}=\sum_{t=0}^{3}\bar{h}_{t}s^{t}(n)s^{3-t}(n-1)$
(7)

where $\left\\{\bar{h}_{t}\right\\}_{t=0}^{3}$ are the new parameters to be
estimated. Rewriting Eqn. (5) in matrix-vector format, we get

$\mathbf{r}=\mathbf{D}\mathbf{v}+\mathbf{e}$ (8)

where $\mathbf{r}=\left[r(0),r(1),\ldots,r(P-1)\right]^{T}$, and $P$ is the
number of observed samples used in the estimation process. Furthermore, the
columns of the matrix $\mathbf{D}$ represent the distortion reference signals
obtained by the expansion of the distortion polynomials. Following the example
in (7), the $t$-th columns of $\mathbf{D}$ is filled by the samples
$\left\\{s^{t}(n)s^{3-t}(n-1)\right\\}_{n=0}^{P-1}$. The length-$\bar{L}$
vector $\mathbf{v}$ contains the new parameters to be estimated representing
the contributions of the reference signals in the columns of $\mathbf{D}$.
Finally, the vector $\mathbf{e}$ is the error vector containing the desired DL
signal $y(n)$, background noise $z(n)$, and the part of the distortion signal
not represented in the columns of $\mathbf{D}$ due to its weak contribution.
The linear least squares (LLS) solution of $\mathbf{v}$ is given by the
solution of the following minimization problem [7]:

$\hat{\mathbf{v}}=\underset{\bar{\mathbf{v}}}{\text{min}}\left\|\mathbf{r}-\mathbf{D}\bar{\mathbf{v}}\right\|^{2}=\mathbf{R}^{-1}\mathbf{q}$
(9)

where $\mathbf{R}=\mathbf{D}^{H}\mathbf{D}$ is the auto-correlation matrix of
the reference signals. Moreover, $\mathbf{q}=\mathbf{D}^{H}\mathbf{r}$
represents the cross-correlation vector between the observed signal and the
reference signals. Transforming the nonlinear estimation problem into a linear
one comes at the expense of increasing the problem dimension, especially for
high distortion levels ($Q$ and $Q_{I}$) and channel length $L$. In HD
schemes, modeling only $x^{Q}(n)$ in the reference matrix $\mathbf{D}$
increases the problem dimension from $L$ to
$\bar{L}=\frac{(Q+L-1)!}{Q!(L-1)!}$. For $Q=3$ and $L=4$, this corresponds to
increasing the problem dimension from $L=4$ to $\bar{L}=20$. The new problem
dimension (after multinomial expansion) becomes even larger for IMD schemes.
In Section IV, we show how sparsity-based techniques are used to reduce the
number of parameters to be estimated.

###  III-B Prior Art and other approaches

The prior-art IMD cancellation approach in [3] models the
$3^{\text{rd}}$-order IMD reference signal as follows:

$p_{\text{prior}}(n)=\left(\sum_{k_{1}=0}^{L_{1}-1}h_{1}(k)s_{1}^{*}(n-k_{1})\right)\left(\sum_{k_{2}=0}^{L_{2}-1}h_{2}(k_{2})s_{2}^{2}(n-k_{2})\right)$
(10)

where the total number of terms in the expansion of (10) is
$J_{\text{prior}}=L_{1}\times L_{2}$. Comparing (10) with the exact
representation in Eqns. (3) and (4), we find the model in [3] is different
from the exact model. The model in [3] squares $s_{2}(n-k_{2})$, while the
squaring should go over $s_{2}$ after passing through the filter $h_{2}(k)$,
i.e., the squaring should go over
$\left(\sum_{k_{2}=0}^{L_{2}-1}h_{2}(k_{2})s_{2}(n-k_{2})\right)$. The IMD
reference signal was modeled in [3] as if $s_{2}(n)$ go through the non-
linearity before getting filtered which is not true as discussed in Section
II.

Another candidate approach would be to model the distortion as if the
nonlinearity is generated before the channel. We call it the Hammerstein-based
approach as it follows the well-known Hammerstein model [8]. For example, the
first nonlinearity term in the RHS of (3) modeling the IMD signal is
approximated by:

$p_{\text{Hammerstein}}(n)=\sum_{k=0}^{J_{\text{Hamm}}-1}h(k)s_{1}^{|p|}(n-k)s_{2}^{|q|}(n-k)$
(11)

For both the prior-art and Hammerstein-based approaches, the corresponding
reference matrix ($\mathbf{D}_{\text{prior}}$ or $\mathbf{D}_{\text{Hamm}}$)
will have $J_{\text{prior}}$ or $J_{\text{Hamm}}$ columns, respectively. Then,
the LLS-based estimation of the distortion signature is then obtained for the
prior-art approach (and similarly for the Hammerstein-based approach) as
follows:

$\hat{\mathbf{v}}_{\text{prior}}=\left(\mathbf{D}_{\text{prior}}^{H}\mathbf{D}_{\text{prior}}\right)^{-1}\mathbf{D}_{\text{prior}}^{H}\mathbf{r}$
(12)

##  IV Proposed Algorithm

We propose a novel approach for improved reference signal design while
controlling the problem dimension. Our approach comprises two main steps.
First, we use the multinomial expansion in (6) and include all the reference
signals out of this expansion in the reference matrix $\mathbf{D}$. The second
step is to obtain a $J_{\text{sparsity}}$-sparse solution of $\mathbf{v}$ in
(8) with only $J_{\text{s}}<<\bar{L}$ nonzero entries. This sparse solution
$\hat{\mathbf{v}}_{s}$ should keep $\mathbf{D}\hat{\mathbf{v}}_{s}$ close to
the observation vector $\mathbf{r}$. This requirement can be casted into the
following optimization problem:

$\hat{\mathbf{v}}_{s}=\underset{\bar{\mathbf{v}}}{\text{min}}\left\|\mathbf{r}-\mathbf{D}\bar{\mathbf{v}}\right\|^{2}\quad\text{subject
to}\quad\left\|\bar{\mathbf{v}}\right\|_{0}=J_{\text{s}}$ (13)

where $\left\|\bar{\mathbf{v}}\right\|_{0}$ is the $l_{0}$-norm of the
argument vector and represents the number of nonzero entries in this vector.
However, the solution of (13) requires an intensive computation burden even
for moderate values of $\bar{L}$ and $J_{\text{s}}$ due to the huge search
space. Several techniques have been proposed in the literature to obtain
approximate solutions of (13), e.g., Orthogonal Matching Pursuit (OMP) [9],
Orthogonal Least Squares (OLS) [10], and FOCUSS [11]. We choose the OMP thanks
to its implementation simplicity and efficiency. The OMP algorithm takes as
input both the reference matrix $\mathbf{D}$, the observation vector
$\mathbf{r}$, and the required sparsity level $J_{\text{s}}$. The OMP output
is an approximate $J_{\text{s}}$-sparse solution
$\hat{\mathbf{v}}_{\text{OMP}}$ of (9) as follows:

$\hat{\mathbf{v}}_{\text{OMP}}=\text{OMP}\left(\mathbf{D},\mathbf{y},J_{\text{s}}\right)$
(14)

The OMP technique is described as follows:

Initialization: Define an empty index set $I_{0}=\phi$, set the initial
residual $\mathbf{r}_{0}=\mathbf{y}$, initialize
$\hat{\mathbf{v}}_{\text{OMP}}=\mathbf{0}$, and set $k=1$.

The $k^{\text{th}}$ iteration:

  1. 1.

Compute $\delta_{i}=\left|\mathbf{r}_{k-1}^{H}\mathbf{D}(:,i)\right|$ for all
$i\notin I_{k-1}$.

  2. 2.

Choose
$c_{k}=\underset{i}{\operatorname{arg}}\;\underset{i}{\operatorname{max}}\;\delta_{i}$.

  3. 3.

Update $I_{k}=I_{k-1}\cup c_{k}$. In this step, the indices of the nonzero
elements are augmented by $c_{k}$, the index of the $k^{\text{th}}$ nonzero
entry computed at the $k^{\text{th}}$ iteration.

  4. 4.

Compute

$\hat{\mathbf{v}}_{\text{OMP}}(I_{k})=\left(\left(\mathbf{D}(:,I_{k})\right)^{H}\mathbf{D}(:,I_{k})\right)^{-1}\left(\mathbf{D}(:,I_{k})\right)^{H}\mathbf{y}$
(15)

where $\hat{\mathbf{v}}_{\text{OMP}}(I_{k})$ holds
$\hat{\mathbf{v}}_{\text{OMP}}$ elements indexed by $I_{k}$.

  5. 5.

Compute
$\mathbf{r}_{k}=\mathbf{y}-\mathbf{D}(:,I_{k})\hat{\mathbf{v}}_{\text{OMP}}(I_{k})$
where $\mathbf{r}_{k}$ is the residual error term at the $k^{\text{th}}$
iteration.

  6. 6.

If $k=J_{\text{s}}$, exit the algorithm, else set $k=k+1$ and go to Step 1.

In words, OMP tries to find the columns (atoms) of the matrix $\mathbf{D}$
(dictionary) whose linear combination is close (matched) to $\mathbf{y}$. From
(15), we find that the maximum size of the matrix to be inverted is
$J_{\text{s}}$, controlled by the designer and can be made much smaller than
$\bar{L}$. We can actually set $J_{\text{s}}=L$ and let the OMP technique
choose the $L$ reference vectors that closely matches the observed signal. The
selected set of reference vectors are adaptive and can be different from a
channel response to another. This is clearly different from [3] where the
reference signals are pre-set regardless of the channel response. The matrix
inversion in (15) is efficiently implemented using the Cholesky decomposition
algorithm [12]. An efficient OMP implementation algorithm using adaptive
Cholesky decomposition is proposed in [13, 14], where the decomposition step
is not performed every iteration. Instead, it is observed that the matrix
$\left(\mathbf{D}(:,I_{k})\right)^{H}\mathbf{D}(:,I_{k})$ is the same as the
matrix $\left(\mathbf{D}(:,I_{k-1})\right)^{H}\mathbf{D}(:,I_{k-1})$ in the
previous iteration except for an extra augmented row and column. This
observation saves computations by using the same Cholesky decomposition of the
previous iteration with one more augmented vector obtained by forward
substitution [15].

##  V Simulation Results and Discussion

We simulate the performance of our proposed approach and compare it with the
other approaches described in Section III-B. We use practical chip responses
provided by well-known manufacturers. In Fig. 1, we simulate the cancellation
performance for the $3^{\text{rd}}$-order IMD mechanism generated by bands 3
and 20 on band 7 as described in Section II-B. The uplink signals of bands 3
and 20 both have bandwidths of 10 MHz, and the receive bandwidth of band 7 is
also 10 MHz. The block size $P$ is set to 520 samples. The IMD level before
cancellation is -85 dBm while the IMD-to-noise power ratio (INR) is set to 0
dB. Fig. 1 shows the IMD power levels before and after cancellation over a
wide range of the DL receive signal power $P_{s}$. Using the same number of
cancellation taps $J_{\text{s}}=J_{\text{prior}}=J_{\text{Hamm}}=9$, our
sparsity-based approach shows 6 and 9 dB improvements over the prior-art
approach in [3] and the Hammerstein-based approach, respectively, described in
Section III-B. For our sparsity-based approach, we construct the columns of
the dictionary matrix $\mathbf{D}$ by setting $L_{1}=L_{2}=3$ in (3) and (4),
hence, $\mathbf{D}$ has 18 columns. For the prior-art approach in [3], we set
$L_{1}=L_{2}=3$, so $J_{\text{prior}}=L_{1}\times L_{2}=9$. In Fig. 1, we also
show the performance of our canceller with full complexity, i.e.,
$J=\bar{L}=18$, where all columns of $\mathbf{D}$ are utilized in IMD
estimation and cancellation. Our approach shows additional 6 dB interference
suppression over prior art in [3]. As expected, the residual IMD power level
gets higher than that of the original IMD at high DL signal power $P_{s}$. The
reason is the poor estimation accuracy due to low IMD to DL signal power ratio
(ISR). However, in practice the original IMD power will not be fixed to -85
dBm regardless of $P_{s}$ as in Fig. 1. Instead, it decreases as $P_{s}$
increases because the user equipment (UE) gets closer to the base station, and
the UL signal is transmitted at lower power levels.

Figure 1: Original and residual IMD power levels for different cancellation
approaches versus DL signal power.

In Figs. 2 and 3, we simulate the performance for the $3^{\text{rd}}$-order HD
mechanism in Section II-A where band 12 generates interference on band 4. The
UL and DL bandwidths of bands 12 and 4, respectively, are both set to 10MHz.
The original HD power level is set to -85 dBm with INR $=0$ dB and $P=520$
samples. Since the prior art in [3] was proposed for IMD but not for HD, we
compare our sparsity-based approach with the Hammerstein-based approach. For
our sparsity-based approach, the columns of the dictionary matrix $\mathbf{D}$
is constructed by setting $L=4$ in (1) and (2), hence $\mathbf{D}$ has 20
columns. In Fig. 2, the residual HD power level is simulated for both
approaches using $J=10$ cancellation taps over a wide range of $P_{s}$. The
performance of the full complexity ($J=20$) is also shown. In Fig. 3, we fix
$P_{s}=-95$ dBm and compare the performances of both approaches for different
number of cancellation taps $J$. The superiority of our sparsity-based
approach is clear, where increasing $J$ improves its performance as it
improves the HD estimation through including more component reference vectors
in the estimation process. However, increasing $J$ does not improve the
performance of the Hammerstein-based approach because its estimation process
does not include the contribution by most of the cross terms of the HD
expansion equation, c.f. Eqn. (7).

Figure 2: Original and residual HD power levels for different cancellation
approaches versus DL signal power. Figure 3: Original and residual HD power
levels for different cancellation approaches versus number of cancellation
taps.

##  VI Conclusion

We proposed a sparse linear filter to cancel the distortion signal generated
by the HD and IMD of uplink signals leaking into the receiver LNA. The
reference signals selection process is dynamic and adapts itself for different
channel responses of the chip through correlating the observed vector with the
dictionary reference signals. Our sparsity-based approach provides additional
interference suppression over the prior-art approach using the same number of
filter taps.

## References

  * [1] “User equipment (UE) radio transmission and reception (Release 12),” _3GPP TS 36.101_ , 2015.
  * [2] B. Razavi, _RF Microelectronics_. NJ, USA: Prentice Hall, 1998.
  * [3] H.-T. Dabag, H. Gheidi, S. Farsi, P. Gudem, and P. Asbeck, “All-digital cancellation technique to mitigate receiver desensitization in uplink carrier aggregation in cellular handsets,” _IEEE Trans. on Microwave Theory and Techniques_ , vol. 61, no. 12, pp. 4754–4765, 2013.
  * [4] M. Kahrizi, J. Komaili, J. Vasa, and D. Agahi, “Adaptive filtering using LMS for digital TX IM2 cancellation in WCDMA receiver,” in _IEEE Radio and Wireless Symposium_ , 2008, pp. 519–522.
  * [5] M. Omer, R. Rimini, P. Heidmann, and J. Kenney, “A compensation scheme to allow full duplex operation in the presence of highly nonlinear microwave components for 4 G systems,” in _IEEE Microwave Symposium Digest (MTT)_ , 2011, pp. 1–4.
  * [6] ——, “All digital compensation scheme for spur induced transmit self-jamming in multi-receiver RF frond-ends,” in _IEEE Microwave Symposium Digest (MTT)_ , 2012, pp. 1–3.
  * [7] S. M. Kay, _Fundamentals Of Statistical Signal Processing: Estimation Theory_. New Jersey: Prentice Hall, 1993\.
  * [8] E. Bai and M. Fu, “A blind approach to Hammerstein model identification,” _IEEE Trans. Sig. Process._ , vol. 50, no. 7, pp. 1610–1619, 2002.
  * [9] G. M. Davis, S. G. Mallat, and Z. Zhang, “Adaptive time-frequency decompositions,” _SPIE J. Optical Eng._ , vol. 33, no. 7, pp. 2183–2191.
  * [10] S. Chen, S. A. Billings, and W. Luo, “Orthogonal least squares methods and their application to non-linear system identification,” _International Journal of control_ , vol. 50, no. 5, pp. 1873–1896, 1989.
  * [11] I. Gorodnitsky and B. Rao, “Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm,” _IEEE Transactions on Signal Processing_ , vol. 45, no. 3, pp. 600–616, 1997.
  * [12] G. H. Golub and C. F. V. Loan, _Matrix Computations-Third Edition_. Baltimore, MD: Johns Hopkins Univ. Press, 1996\.
  * [13] S. Singhal and B. Atal, “Amplitude optimization and pitch prediction in multipulse coders,” _IEEE Transactions on Acoustics, Speech and Signal Processing_ , vol. 37, no. 3, pp. 317–327, 1989.
  * [14] S. Cotter, J. Adler, B. Rao, and K. Kreutz-Delgado, “Forward sequential algorithms for best basis selection,” _IEE Proceedings-Vision, Image and Signal Processing_ , vol. 146, no. 5, pp. 235–244, 1999.
  * [15] R. Rubinstein, M. Zibulevsky, and M. Elad, “Efficient implementation of the K-SVD algorithm using batch orthogonal matching pursuit,” _CS Technion report_ , vol. 40, no. 8, pp. 1–15, 2008.