1

Newtonized Orthogonal Matching Pursuit for Line
Spectrum Estimation with Multiple Measurement
Vectors
Lin Han, Jiang Zhu, Rick S. Blum and Zhiwei Xu

arXiv:1802.01266v1 [cs.IT] 5 Feb 2018

Abstract
A Newtonized orthogonal matching pursuit (NOMP) algorithm is proposed to estimate continuous frequencies and amplitudes
of a mixture of sinusoids with multiple measurement vectors (MMVs). The proposed algorithm includes two key steps: Detecting
a new sinusoid on an oversampled discrete Fourier transform (DFT) grid and refining the parameters of already detected sinusoids
to avoid the problem of basis mismatch. We provide a stopping criterion based on the overestimating probability of the model
order. In addition, the convergence of the proposed algorithm is also proved. Finally, numerical results are conducted to investigate
the effectiveness of the proposed algorithm when compared against the state-of-the-art algorithms in terms of frequency estimation
accuracy and run time.

Index terms— Orthogonal matching pursuit, frequency estimation, line spectrum, Newton refinement, multiple measurement
vectors
I. I NTRODUCTION
One of the classical problems in digital communication and radar processing applications is to estimate continuous-valued
frequencies of sinusoids in additive white Gaussian noise (AWGN) environments from a small number of measurements [1, 2].
On the one hand, several classical subspace methods have been proposed to perform the frequency estimation, such as MUSIC
and ESPRIT [3, 4], which exploit the autocorrelation matrix’s low-rank structure to estimate the underlying signal subspace. As
the signal-to-noise ratio (SNR) decreases, the performance of these two algorithms for estimating closely spaced frequencies
will degrade [5]. On the other hand, a variety of methods based on sparse representation and compressed sensing (CS) have
also been proposed to estimate frequencies for multiple sinusoids [6, 7]. Basically, the estimation problem can be transformed
to that of seeking a sparse approximation of the multiple sinusoids by referring to an infinite-dimensional dictionary. In fact,
of all the frequencies lie on the discrete Fourier transform (DFT) grid, it can be shown that the signal can exactly be recovered
by utilising convex optimization from randomly selected samples with high probability [8].
However, there exists a major grid mismatch problem induced by the fact that the measurements are sparsely represented
under a finite discrete dictionary, which badly deteriorates the performance of various reconstruction algorithms. In fact, one has
to make a reasonable tradeoff between the oversampling rate and the computational cost when implementing sparse methods.
This unavoidable grid mismatch problem is studied in [9–11] in detail. Moreover, sparse reconstruction methods usually entail
one or more parameters which in fact are not necessarily known, such as the number of the sinusoids, the regularization
parameters, the variance of the noise and so on. Recently, the semiparametric iterative covariance-based estimation (SPICE)
algorithm [12–14] has been proposed to alleviate the drawbacks of the discretization operation to a great extent, which uses
the covariance fitting criterion from a statistical perspective and no user-parameters are required.
A. Related work
An exact discretization-free method, called the sparse and parametric approach (SPA) for uniform and sparse linear arrays,
is proposed in [15]. This algorithm utilises the advantage of multiple measurement vectors (MMVs) to perform line spectrum
estimation by solving a semidefinite programming problem. In [16], two approaches are developed to solve the problem of
line spectrum denoising and estimation, which estimate an ensemble of spectrally-sparse signals composed of the same set
of continuous-valued frequencies from MMVs. The iterative reweighed approach (IRA) is proposed to deal with both single
measurement vector (SMV) and MMV settings [17], where all the frequencies are updated in parallel with good initial frequency
estimates.
The most related work is the Newtonized orthogonal matching pursuit (NOMP) algorithm proposed in [18] considered for
SMV settings. The authors provided a constant false alarm rate (CFAR) based on the stopping criterion, a detailed convergence
analysis result and a comprehensive comparison with state-of-the-art algorithms. The work presented here has focused on the
formulation of the NOMP algorithm with MMVs to estimate frequencies over the continuum, which is especially important
in array signal processing fields. We not only extend the NOMP to MMVs scenarios, but also characterize convergence and
provide a stopping criterion based on the given probability of overestimating the model order.

2

B. Contributions
The main contributions are summarized as follows:
• We extend the NOMP algorithm to deal with line spectrum estimation in the MMVs setting. Similar to the NOMP [18],
our proposed algorithm avoids the basis mismatch problem by using a Newton refinement step as feedback to improve
the estimation of already detected sinusoids. Besides this, the model order K is not needed in advance during estimation
of frequencies of the sinusoids in the mixture, and a stopping criterion based on probability is provided.
• We analyze the convergence of the NOMP algorithm in MMVs settings by using the theory of dual norms. Specifically,
we provide an upper bound on the number of iterations required by the NOMP method with MMVs to meet the provided
stopping criterion and also obtain a bound on the convergence rate.
• By implementing numerical simulations to compare the estimation performance of various algorithms against the CramérRao bound (CRB) in different scenarios, we show that the algorithm can achieve a near-optimal performance in terms of
estimation accuracy.
• We analyze the computation time in four scenarios conducted in the simulation section. It can be shown that the proposed
algorithm is faster than SPA, IRA and basis pursuit de-noising (BPDN), but slower than the classical MUSIC method.
Outline: In Section II, we set up our problem model. We propose our algorithm in Section III. In Section IV, we present the
stopping criterion based on the probability of overestimating the model order, along with an analytical expression of the miss
probability of detecting the sinusoids. We present the analysis of convergence in Section V. In Section VI, we conduct numerical
experiments to compare the estimation accuracy of the NOMP with MMVs against the state-of-the-art methods. Section VII
concludes the paper. Notation: We use (·)H , (·)T and (·)∗ to denote the conjugate transpose, transpose and conjugate operator
respectively. CN denotes the complex Gaussian distribution. The Frobenius norm, the real part of the complex number a and
the trace operator are denoted by k·k2F , R{a} and tr{·} respectively. bac denotes the greatest integer that is less than or equal
to a. ⊗ is the Kronecker product operator.
II. P ROBLEM S ETUP
In an MMV model, we consider a line spectrum estimation scenario with N sensors and T snapshots collecting measurements
of K distinct frequency components. The measurements at the array output can be expressed as
Y = ΦAX + Z,

(1)

iT
1 h jωl
1, e , · · · , ej(N −1)ωl ,
a(ωl ) , √
N

(2)

where the lth column of A is defined as

where Y ∈ CM ×T is the noisy measurement collected by all N sensors and T snapshots, Φ ∈ CM ×N is the known compression
matrix and A = [a(ω1 ), · · · , a(ωK )] ∈ CN ×K . Each ωl is continuous-valued in [0, 2π). Furthermore, zij is independent and
identically distributed (i.i.d.) Gaussian random variable and follows zij ∼ CN (0, σ 2 ). X = [x1 , · · · , xT ] ∈ CK×T contains all
the sidnusoid amplitudes xij for each snapshot. The uniform linear array (ULA) scenario for DOA can also be formulated as
model (1) [15].
Before moving forward we make the following assumptions on the statistics of the sources and noises:
1) Sources are uncorrelated.
2) Sources and noises are independent.
3) Noises at different sensors are uncorrelated.
III. A LGORITHM
We first look into the estimation problem of a single sinusoid in the compressed scenarios, and then generalize the results
to a mixture of sinusoids.
A. Single frequency
In this scenario, the model (1) simplifies to
Y = ΦaxT + Z,
(3)

T √
T
where Φ ∈ CM ×N , x = [x1 , · · · , xT ] ∈ CT ×1 , and a = 1, ejω , · · · , ej(N −1)ω / N .
The Maximum Likelihood (ML) estimate of the amplitudes x and frequency ω can be calculated by minimizing the residual
power kY − ΦaxT k2F , which is equal to maximizing the function
S(x, ω) = 2tr{R(YH aΦ xT )} − kaΦ xT k2F ,
where aΦ = Φa.

(4)

3

It’s difficult to directly optimize S(x, w) over all amplitudes and frequency. As a result, a two stage procedure is adopted:
(1) Detection stage, in which a coarse estimate of ω is found by restricting it to a discrete set, (2) Refinement stage, where
we iteratively refine the estimates of amplitudes and frequency of detected sinusoids. By using the least squares method
procedure, we can then obtain the expression of xT that maximizes S(x, ω) for any given ω, i.e., by setting ∂S(x, ω)/∂x∗ =
H
T
aH
Φ Y − aΦ aΦ x = 0 [19], we obtain
−1 H
2
xT = (aH
aΦ Y = aH
Φ aΦ )
Φ Y/kaΦ k2 .

(5)

Substituting xT in S(x, w) allows us to find the ML estimate of ω, which is the solution to the following optimization
problem:
ω̂ = argmax GY (ω),

(6)

ω

where
GY (ω) = kYH aΦ k22 /kaΦ k22 .

(7)

Detection: By restricting ω to a finite discrete set denoted by Ω , {k(2π/γN ) : k = 0, 1, · · · , (γN − 1)}, where γ is the
oversampling factor relative to the DFT grid, we can obtain a coarse estimate of ω. We treat the ωc ∈ Ω that maximizes the
2
cost function (7) as the output of this stage, and the corresponding xT vector estimate is aH
Φ (ωc )Y/kaΦ (ωc )k2 .
Refinement: Let (x̂, ω̂) denote the current estimate, then the Newton procedure for frequency refinement can be denoted as
ω̂ 0 = ω̂ − ĠY (ω̂)/G̈Y (ω̂),

(8)

where
( (
ĠY (ω̂) = tr R Y
G̈Y (ω̂) =
( (
YH
tr R
kaΦ k22
−

YH aΦ aH
ΦY
2
(kaΦ k2 )2

H



daΦ H
daH
aΦ + aΦ Φ
dω
dω



Y/kaΦ k22

−Y

H

aΦ aH
ΦY



daH
daΦ
Φ
aΦ + aH
Φ
dω
dω



))
/(kaΦ k22 )2

,

(9)

  H

 2 H !

d2 aΦ H
daΦ daH
daΦ
d aΦ
2YH
daΦ H
daH
H daΦ
Φ
Φ
a +2
+ aΦ
a + aΦ
Y
aΦ + aΦ
Y−
dω 2 Φ
dω dω
dω 2
(kaΦ k22 )2 dω Φ
dω
dω
dω
!
 2 H
 H
2 ))
H
H
H
2
d aΦ
2Y aΦ aΦ Y daΦ
da daΦ
d aΦ
daΦ
+ aH
+
aΦ + aH
aΦ + 2 Φ
.
(10)
Φ
Φ
2
2
2
3
dω
dω dω
dω
(kaΦ k2 )
dω
dω

We maximize GY (ω) by employing the update rule (8) on the condition that the function is locally concave.
B. Multiple frequency
Assume that we have already detected k sinusoids, and let P = {(xl , wl ), l = 1, · · · , L} denote the set of estimates of the
detected sinusoids. The residual measurement corresponding to this estimate can be given by
Yr (P ) = Y −

L
X

aΦ (ωl )xT
l

(11)

l=1

The method of estimating multiple frequencies proceeds by employing the single sinusoid procedure to perform Newtonized
coordinate descent on the residual energy kYr (P )k2F . One step of this coordinate descent involves adjusting all ωl . The
procedure to refine the lth sinusoid is as follows: Yr (P \{xl , ωl }) now is referred to as the measurement Y and the single
frequency update step is utilised to refine (xl , ωl ).
Refinement Acceptance Condition (RAC): This refinement step is accepted when it results in a strict improvement in
GYr (P \{xl ,ωl }) (ω), namely, GYr (ω̂ 0 ) > GYr (ω̂). By doing this, we can make sure that the adopted refinement must decrease
the overall residual energy.
In summary, firstly, we detect a frequency ω̂ over the discrete set Ω by maximizing the cost function (7). Then we use the
knowledge of the first-order and second-order derivative of the cost function to refine the estimate of ω̂. Next, we use the
information of all the other previously detected sinusoids to further improve the estimation performance of every previously
detected sinusoid one at a time. This step is crucial for the convergence and accuracy of the algorithm. Finally, we update x
by least squares methods. The whole algorithm is summarized in Algorithm 1.

4

Algorithm 1 Multi snapshot NOMP.
1: Procedure EXTRACTSPECTRUM (Y, τ ) :
2: m ← 0, P0 = {}
3: while maxω∈DFT GYr (Pm ) (ω) > τ do
4:
m←m+1
5:
IDENTIFY
ω̂ = arg maxω∈Ω GYr (Pm−1 ) (ω)
and its corresponding x vector
estimate

2
x̂T ← aH
(ω̂)Y
(P
)
/ka
r
m−1
Φ (ω̂)k2 .
Φ
0
6:
Pm
← Pm−1 ∪ {(x̂, ω̂)}
7:
SINGLE REFINEMENT: Refine (x̂, ω̂) using single frequency Newton update algorithm (Rs Newton steps) to obtain
improved estimates (x̂0 , ω̂ 0 ).
00
8:
Pm
← Pm−1 ∪ {(x̂0 , ω̂ 0 )}
00
00
00
9:
CYCLIC REFINEMENT: Refine parameters in Pm
one at a time: For each (x, ω) ∈ Pm
, we treat Yr (Pm
\{(x, ω)}) as
the measurement Y, and apply single frequency Newton update algorithm. We perform Rc rounds of cyclic refinements.
000
Let Pm
denote the new set of parameters.
00
000
10:
UPDATE all x vector estimate in Pm
by least squares: A , [aΦ (ω1 ), · · · , aΦ (ωm )], {ωl } are the frequencies in Pm
.
T
†
And [x1 , · · · , xm ] = A Y.
Let Pm denote the new set of parameters.
11: end while
12: return Pm

IV. S TOPPING CRITERION
To understand the performance of Algorithm 1, the probability of the algorithm overestimating the model order K is of
interest. An extreme scenario is that the proposed algorithm has detected K sinusoids which causes the residual to be only
AWGN in the model (1), and the stopping criterion still isn’t met. So the algorithm has to detect another sinusoid to make the
residual decrease, which corresponds to the scenario of overestimating the model order.
We use the stopping criterion to estimate the model order K. If the residual energy can be well explained by noise, up to a
target overestimating probability, then we stop. Intuitively, we choose to terminate the algorithm by comparing the magnitude
of the Fourier transform of the residual with the expected noise power. Details are given in the next section.
A. Stopping criterion
The algorithm stops when
GYr (P ) (ω) =

T
X

2

|hyri (P ), a(ω)i| < τ

(12)

i=1

for all DFT sampling frequencies {ωn , 2πn/N : n = 0, · · · , N − 1}, where
PTyri (P ) is the ith column of Yr (P ), T is the
number of snapshots and τ is the stopping threshold. Thus, we stop when i=1 kF {yri (P )}k2∞ < τ , where F {yri (P )}
denotes the discrete Fourier transform of yri (P ).
Supposedly, we have already correctly detected all sinusoids in the mixture. Under this condition, the residual is yri (P ) ≈ zi ,
where zi ∼ CN (0, σ 2 IM ). Note
PTthat F {·} corresponds to a projection operation, so the statistics of AWGN are unchanged
by it. Then by defining Rn , i=1 kaH (ωn )zi k22 , we obtain
( T
)




X
2
kF {yri (P )}k∞ > τ = Pr
max Rn > τ = 1 − Pr
max Rn ≤ τ = 1 − Pr (R1 ≤ τ, · · · , RN ≤ τ )
Pr
i=1

n=1,··· ,N

n=1,··· ,N

(13)
Note that by defining un,i = aH (ωn )zi , we have



 H

2
2
E un1 ,i1 u∗n2 ,i2 = E aH (ωn1 )zi1 zH
i2 a(ωn2 ) = σ δi1 ,i2 E a (ωn1 )a(ωn2 ) = σ δi1 ,i2 δn1 ,n2 ,

(14)

5

0.12

Measured Poe [percent]

0.1

30dB
32dB
34dB

0.08
0.06
0.04
0.02
0
0

0.02

0.04
0.06
0.08
Nominal Poe [percent]

0.1

0.12

Fig. 1. Nominal vs measured probability of overestimating probability.

where δi1 ,i2 denotes the dirichlet function, which equals zero unless i1 = i2 holds. From (14), we can conclude that Rn is a
χ2 random variable with 2T degrees of freedom. With the degrees of freedom being even, the cumulative distribution function
(CDF) has a closed form given by


k
T
−1
X

1 Rn
 1 − e− Rσn2
, Rn > 0
k! σ 2
F2T (Rn , σ 2 ) =
(15)
k=0


0, otherwise.
Therefore, (13) can be calculated as
( T
)
X
N
Pr
kF {yri (P )}k2∞ > τ = 1 − [Pr (Rn ≤ τ )] = 1 − F2T (τ, σ 2 ) = 1 − F2T (τ /σ 2 , 1),

(16)

i=1

where F2T (Rn , 1) denotes the CDF of χ2 random variable with unit variance.
Let Poe denote a nominal overestimating probability. Thus Poe satisfies
( T
)
X
2
Pr
kF {yri (P )}k∞ > τ = Poe .

(17)

i=1
−1
By defining F2T
(·) as the inverse function of F2T (·, 1), we obtain
i
h

1
−1
τ = σ 2 F2T
(1 − Poe ) N /2.

(18)


2
Note that for a single snapshot, i.e., T = 1, F2T (Rn , σ 2 ) = 1 − e−Rn /σ and τ = −σ 2 log 1 − (1 − Poe )1/N from (18),
which is consistent with the results obtained in [18].
We conduct a numerical experiment by comparing the “measured” against “nominal” overestimating probability (17) to
substantiate the above analysis. We use Algorithm 1 to estimate frequencies in a mixture of K = 16 sinusoids of the same
fixed nominal SNR, defined as


2
ka(ωl )xT
l kF
SNRl = 10 × log10
, l = 1, · · · , K.
(19)
σ2 T
Furthermore, we generate the frequencies such that the minimal wrap-around frequency separation is ∆ωmin , where ∆ωmin =
∆ωDFT and the ∆ωDFT , 2π/N is the DFT grid separation. The parameters are set as follows: N = 64, T = 10, Rs = 1,
Rc = 3, the number of Monte Carlo (MC) trials is 300. The “measured” overestimating probability is defined as the ratio of
the overestimating events in all MC trials. From Fig. 1, it can be been that the empirical overestimating probability is close to
the nominal value at various SNRs, which means that the final residual error can be approximated as the AWGN, namely, the
frequency estimation accuracy of Algorithm 1 is good.
B. Probability of Miss
Firstly, we define the neighborhood Nωi around the true frequency ωi as Nωi , {ω : dist(ω, ωi ) ≤ 0.25 × ∆DFT }. Then
we declare a miss of ωi if none of the estimated frequencies locates in Nωi , otherwise we declare a successful detection of
ωi . Note that a miss is caused by both noise and inter-sinusoid interference, but we only discuss noise here.

6

Assuming no inter-sinusoid interference, the measurement Y can be described as Y = [a(ω)x1 +z1 ; a(ω)x2 +z2 ; · · · ; a(ω)xT +
zT ]. The signal frequency ω is not detected if

Pmiss = Pr kaH (ωd )Yk22 < τ ,
(20)
where ωd denotes the sampling frequency. Note that xi is deterministic known, and zi ∼ CN (0, σ 2 IN ). Hence, we have
( T
)
X
2
H
Pmiss = Pr
a (ωd ) [a(ω)xi + zi ] < τ .
(21)
i=1

To calculate the above probability, we introduce
vi , aH (ωd ) [a(ω)xi + zi ] = xi e

(N −1)(ω−ωd )
j
2





N (ω−ωd )
2
ω−ωd
N sin( 2 )

sin

+ aH (ωd )zi .

(22)

By defining
α,
R̃m

sin [N (ω − ωd )/2]
,
N sin [(ω − ωd )] /2
T
X
H
2
, ka (ωd )Yk2 =
|vi |2 ,

(23a)
(23b)

i=1

we can conclude that R̃m is a noncentral χ2 random variable with 2T degrees of freedom and common variance being σ 2 /2.
With degrees of freedom being an even number, the CDF of R̃m can be written in the form
!
p

√

2s
2R̃m
1 − Q
,
, R̃m > 0
T
σ
σ
(24)
F2T (R̃m ) =


0,
where QT denotes Marcum Q-function and s is defined as
v
u T
uX
s , αt
|xi |2 .

(25)

i=1

Thus it’s easy to show that
Pmiss

 v

u T
u X
p
= 1 − QT αt2
|xi |2 /σ, 2τ /σ 2  .

(26)

i=1

Supposing a frequency within a DFT grid interval follows the uniform distribution, then we obtain E[α] = 0.88, where
ω ∼ U[−π/N, π/N ]. Hence,
v


u T
u X
p
Pmiss = 1 − QT 0.88t2
|xi |2 /σ, 2τ /σ 2  .
(27)
i=1

V. C ONVERGENCE
In this section, the convergence of the proposed algorithm is studied. Firstly, upper bounds on the number of iterations
needed to reach the stopping condition are given. Then a bound on the rate of convergence of Algorithm 1 is provided, which
is a function of the “atomic norm” of original measurements Y and the oversampling factor γ. Here, we only consider the
uncompressed scenario, namely, Φ = IN .
A. Proof of convergence
It’s easy to show that the number of measurements N is a trivial upper bound of the number of iterations of Algorithm 1.
From step 10 of Algorithm 1, it can be shown that X becomes a square full-rank matrix (note that there is no frequency that
will be detected twice) after N iterations. Then the algorithm terminates because the residual is equal to zero for (N + 1)th
iteration.
In the following, we provide another upper bound on the number of iterations, which is obtained by considering how much
the residual energy will decrease when a new frequency is added to the set of estimated sinusoids.

7

Theorem 1 Let kYk2F be the original residual energy. Whena new frequency is added to the set of estimated sinusoids, the
reduction of residual energy is at least τ . Consequently, min N, bkYk2F /τ c is an upper bound on the number of iterations
of the Algorithm 1.
Pm
P ROOF Assume we have detected m sinusoids, the residual measurement is given by Yr (Pm ) = Y − l=1 a(ωl )xH
l . The
residual energy in each iteration step of the Algorithm 1 satisfies the following,
(b)

(a)

(c)

0
000 2
)k2F + GYr (Pm−1 ) (ω̂) ≥ kYr (Pm
)kF + GYr (Pm−1 ) (ω̂) ≥ kYr (Pm )k2F + GYr (Pm−1 ) (ω̂)
kYr (Pm−1 )k2F = kYr (Pm

(28)

(d)

≥ kYr (Pm )k2F + τ.

(29)

where ω̂ denotes the detected frequency. Equality in (a) holds because of the step 5 in the Algorithm 1 where we project

H
0
Yr (Pm−1 )orthogonal to the subspace spanned by a(ω̂) to get Yr (Pm
). Basically, we get Yr (Pm−1 ) − a(ω̂)xH
a(ω̂) = 0
m
H 2
H
2
kY
(P
)
−
a(ω̂)x
k
.
Then
we
take
the
Frobenius
norm
of both
and xH
=
a(ω̂)
Y
(P
)
/ka(ω̂)k
by
solving
min
r
m−1
r
m−1
m F
m
2
xm

0
sides of Yr (Pm
) = Yr (Pm−1 ) − a(ω̂)xH
m to obtain equality (a). Inequality in (b) follows from RAC, which is implemented
whenever the single refinement step is executed and (c) is a direct consequence of the step 10 of the Algorithm 1, which can
only cause a decrease in the residual energy. And the stopping criterion leads to (d).
From inequality (29), we can conclude that the reduction of the residual energy caused by the detection of a new sinusoid
frequency is always larger than τ . Thus we can get another bound on the number of iterations of Algorithm 1, which is
bkYk2F /τ c.

B. Rate of convergence
In the uncompressed and noiseless scenario, the noiseless measurement Yo can be written as
Yo =

K
X

a(ωk )xH
k =

k=1

K
X

ck a(ωk )φH
k.

(30)

k=1

T ×1
where ck = kxk k2 > 0, ωk ∈ [0, 2π), φk = c−1
with kφk k2 = 1. Then the continuous dictionary or the set of
k xk ∈ C
atoms is given as

A , A(ω, φ) = a(ω)φH : ω ∈ [0, 2π), kφk2 = 1 .
(31)

It is easy to see that Yo is a linear combination of many atoms in A. Here, we define the atomic `0 (pseudo-) norm of
Y ∈ CN ×T as the smallest number of atoms that can express it:
(
)
K
X
kYkA,0 = inf Y =
ck A(ωk , φk ), ck ≥ 0 .
(32)
K

k=1

Because the atomic `0 norm is non-convex, we utilize convex relaxation to relax the atomic `0 norm to the atomic norm, which
is defined as the gauge function of conv(A), where conv(A) is the convex hull of A:
(
)
X
X
kYkA , inf{t > 0 : Y ∈ tconv(A)} = inf
ck : Y =
ck A(ωk , φk ), ck ≥ 0 .
(33)
k

k

By defining hY, Ai = tr(A Y) and hY, AiR = R(hY, Ai), the dual norm of kYkA can be written as
H

kYk∗A ,
=

sup
kAkA ≤1

hY, AiR
hY, a(ω)φH iR =

sup
ω∈[0,2π),kφk2 =1

|hφ, YH a(ω)i| =

sup
ω∈[0,2π),kφk2 =1

sup

kYH a(ω)k2 =

ω∈[0,2π)

sup

p
GY (ω), (34)

ω∈[0,2π)

where the GY (ω) is (7).
Theorem 2 Maximizing GY (ω) (for the dictionary of unit norm sinusoids) over [0, 2π) is consistent with that over the
oversampled grid Ω with oversampling factor γ. Namely, we have

−1/2
p
p
p
4πT
max GY (ω) ≤ sup
GY (ω) ≤ 1 −
max GY (ω)
(35)
ω∈Ω
ω∈Ω
γ
ω∈[0,2π)
2

P ROOF According to (34), (kYk∗A ) can be expressed as
2
(kYk∗A )

=

sup
ω∈[0,2π)

kY

H

a(ω)k22

=

sup

T
X

ω∈[0,2π) t=1

N
1 X ∗ j(n−1)ω
√
φn,t e
N n=1

2

=

sup

T
X

ω∈[0,2π) t=1

2

|Wt (ω)| ,

(36)

8

PN
where φn,t is the (n, t)th entry of Y, and Wt (ω) , √1N n=1 φ∗n,t ej(n−1)ω . For ω1 , ω2 ∈ [0, 2π), according to Bernstein’s
theorem [20] and the result in [21] (Appendix C), we can obtain
|Wt (ω1 )| − |Wt (ω2 )| ≤ 2N |ω1 − ω2 | sup

|Wt (ω)|.

(37)

ω∈[0,2π)

Then we have
T
X

2

|Wt (ω1 )| −

t=1

T
X

!

T
X
|Wt (ω2 )| ≤
(|Wt (ω1 )| + |Wt (ω2 )|) 2N |ω1 − ω2 | sup
2

(a)

t=1

|Wt (ω)|

ω∈[0,2π)

t=1
T
X

(b)

≤ 4N |ω1 − ω2 |T sup

2

|Wt (ω)|2 = 4N |ω1 − ω2 |T (kYk∗A ) ,

(38)

ω∈[0,2π) t=1

where inequality (a) follows from plugging in (37), inequality (b) is from the fact that
PT

sup

|Wt (ω)|2 ≤ T sup

PT

|Wt (ω)|2 , and the last equality follows by (36).
o
n
−1)
2π
, we have
, · · · , 2π(γN
Let ω2 take any value of the grid points 0, γN
γN
t=1

ω∈[0,2π)

t=1

ω∈[0,2π)

2

(kYk∗A ) =
≤

T
X

sup

2

|Wt (ω)| =

ω∈[0,2π) t=1
T 
X

sup

sup

T 
X

2

2

|Wt (ω)| − |Wt (e
ω )| + |Wt (e
ω )|

2



ω∈[0,2π) t=1
2

2

|Wt (ω)| − |Wt (e
ω )|



ω∈[0,2π) t=1

+

T
X

max

d=0,··· ,γN −1

T
X



2πd
γN
t=1


T
X
2π
2πd
∗ 2
≤ 4N
Wt
T (kYkA ) +
max
d=0,··· ,γN −1
2γN
γN
t=1


T
2
X
4πT
2πd
2
=
Wt
(kYk∗A ) +
max
d=0,··· ,γN −1
γ
γN
t=1
2

≤ 4N |ω − ω
e |T (kYk∗A ) +

max

d=0,··· ,γN −1


Wt

t=1
2

2πd
γN



2



Wt

2

(39)

2

Since the maximum on the grid is a lower bound for (kYk∗A ) , we have


 2 !1/2
−1/2
T
X
2πd
4πT
max
Wt
≤ kYk∗A ≤ 1 −
d=0,··· ,γN −1
γN
γ
t=1

T
X

max

d=0,··· ,γN −1


Wt

t=1

2πd
γN

 2 !1/2
(40)

Thus,
max
ω∈Ω

p

GY (ω) ≤

sup

p


GY (ω) ≤

ω∈[0,2π)

4πT
1−
γ

−1/2
max
ω∈Ω

p

GY (ω).

(41)

To prove Theorem 3, the following lemma in [22] is introduced
Lemma 1 [22] Assume {an }n≥0 is a decreasing sequence of nonnegative numbers such that a0 ≤ U and

an−1 
an ≤ an−1 1 −
, ∀n > 0,
U
then we have an ≤

U
n+1

(42)

for all n ≥ 0.

Theorem 3 For all Y satisfying kYkA < ∞, the residual energy of Algorithm 1 at the mth iteration satisfies

−1/2
4πT
−1/2
kYr (Pm )kF ≤ (m + 1)
1−
kYkA .
γ

(43)

P ROOF From (28), we have
kYr (Pm )k2F ≤ kYr (Pm−1 )k2F − GYr (Pm−1 ) (ω̂).

(44)

9

Yr (Pm−1 ) is a direct consequence of projecting Y orthogonal to the subspace spanned by Pm−1 , therefore
q
(a)
kYr (Pm−1 )k2F = hYr (Pm−1 ), YiR ≤ kYkA kYr (Pm−1 )k∗A = kYkA sup
GYr (Pm−1 ) (ω)
ω∈[0,2π)



(b)

≤ kYkA 1 −

4πT
γ

−1/2

q
max GYr (Pm−1 ) (ω),

(45)

ω∈Ω

where inequality (a) follows from the following
kYr (Pm−1 )k∗A =


sup
kYkA ≤1

hY, Yr (Pm−1 )iR ≥

Y
, Yr (Pm−1 )
kYkA


(46)
R

and equality (b) is obtained from Theorem 2. From the step 5 of the Algorithm 1, we have
q
ω̂ = arg max GYr (Pm−1 ) (ω),

(47)

ω∈Ω


By defining η , kYkA 1 −

4πT
γ

−1/2

and combing with (47), we have
q
kYr (Pm−1 )k2F ≤ η GYr (Pm−1 ) (ω̂).

(48)

Combining (44) and (48), yields

kYr (Pm )k2F ≤ kYr (Pm−1 )k2F 1 − η −2 kYr (Pm−1 )k2F .

(49)

By using Lemma 1 and the fact that
(a)

kYr (P0 )k2F = kYk2F ≤ kYk2A ≤ η 2 ,

(50)

where inequality (a) follows from (51)
)
)
(K K
(K K
K
K X
X
XX
XX
H
|xH
|xH
≤
xk aH (ωk )a(ωl )xH
= tr
kYk2F = tr{YH Y} = tr
l xk |
l xk ||a (ωk )a(ωl )|
l
≤

k=1 l=1

k=1 l=1

k=1 l=1
K
K X
X

kxk k2 kxl k2 = kYk2A

(51)

k=1 l=1

we have
kYr (Pm )k2F ≤

η2
.
m+1

(52)

In other words,
kYr (Pm )kF ≤ (m + 1)

−1/2


−1/2
4πT
1−
kYkA .
γ

(53)

This proves the Theorem 3.
−1/2

−1

Note that the component in (53) is (1 − 4πT /γ)
instead of (1 − 2π/γ) in [18]. Given T = 1, our conclusion isn’t
−1/2
−1
consistent with that in [18]. In fact, due to (1 − 4πT /γ)
≥ (1 − 2π/γ) , the bound on the rate of convergence is tighter
in the single snapshot scenario.
C. Empirical rate of convergence
The relative residual energy of the ith iteration (averaged over 300 runs) versus the number of iterations in a noiseless scenario
is plotted to show the effects of the refinement steps on the convergence of Algorithm 1, defined as 20log10 (kYr (Pm )kF /kYkA ).
Here, we compare Algorithm 1 in the scenarios of T = 1 amd T = 10 to the following variants of OMP to show the
improvements brought by the refinement steps.
NOMP without cyclic refinements in MMV setting: This is an algorithm that has a nearly comparable performance with OMP
[23] over the continuum of atoms. We use Algorithm 1 to implement this method by setting the number of cyclic refinement
steps to 0. This algorithm lies in the class of forward greedy methods because it doesn’t have a feedback mechanism. Hence,
our analysis is also applicable to this method.
Discretized OMP: If we skip the single and cyclic refinement steps, then we can obtain a standard OMP applied to the
oversampled grid Ω. Note that this algorithm can be viewed as a special case of Algorithm 1, so the convergence analysis

10

0
NOMP (T = 1)
NOMP (T = 10)
DOMP
NOMP−

−20

Relative residual power (dB)

−40
−60
−80
−100
−120
−140
−160
−180
−200
0

10

20

30
Iteration count

40

50

60

Fig. 2. Convergence rates at noiseless case.

results are also valid. The parameters are set as follows: K = 16, N = 64, T = 10, Rs = 1, Rc = 1, SNR = 25 dB, ∆ωmin =
2.5∆ωDFT . From Fig. 2, it can be shown that a small oversampling factor (γ = 4) with single refinement step (NOMP-)
has a slightly faster convergence than a large oversampling factor (γ = 20) with no refinements (DOMP). Furthermore, we
can see a large gap between Algorithm 1 and the other two algorithms, which means that the cyclic refinement steps lead
to an extremely fast convergence for Algorithm 1. In Fig. 2, we can see that the residual energy drops sharply at the 16th
iteration, which is equal to the number of sinusoids in the mixture K. Furthermore, it can be seen that increasing the number
of snapshots has a slight effect on the convergence rate.
VI. S IMULATION
In this section, we conduct numerical simulations to compare the performance of the proposed Algorithm 1 against other
methods in terms of estimation accuracy in various scenarios.
Benchmarks: All algorithms are compared against the DFT method implemented by coarsely picking out the top K peaks,
and the oversampling rate for DFT is set to 4.
For the classical subspace method, we choose the MUSIC method to estimate frequencies of sinusoids in the mixture.
Specifically, we make one modification based on the version of the MATLAB rootmusic algorithm. In detail, a sliding window
of size W is adopted to obtain an estimate of the covariance matrix of the measurement data. Different choices of W have a
considerable impact on the estimation accuracy. In particular, if W is set too large, the estimation performance of the signal
subspace and the autocorrelation matrix will degrade. Whereas a small W will result in a less accurate frequency estimation
performance. Empirically, we found that W = 44 results in the best estimation performance in our all four scenarios. It’s worth
mentioning that, instead of implementing an empirical method to estimate K, we just directly feed the true K to the solver.
Newtonized BPDN: For the sparse method inspired by sparse representation, we employ the SPGL1 toolbox [24] to solve
the l2,1 minimization problem, also known as the MMV version of BPDN, we denote this method simply as BPDN in the
following, which is defined as
min kXk1,2 subject to kAX − YkF ≤ τ,

(54)

where kXk1,2 denotes the sum of the two-norms of the rows of X and τ is the measure of the noise level. Newtonized BPDN
can be viewed as an extension of the BPDN
√ method. By applying this toolbox to model (1) in the scenario of Φ = IN , where
the oversampling rate is set to 4 and τ = N T σ 2 , we then obtain the optimized X. By sorting the two-norms of every row
of the estimated X in a descending order, we choose the K frequencies that correspond to the top K two-norms of every row
of the estimated X as the estimated frequencies. To avoid the frequency splitting phenomenon [25, 26], we impose an extra
procedure to the Newtonized BPDN method to cope with two special cases. On the one hand, when the interval between two
adjacent estimated frequencies sorted in an ascending order is smaller or equal to 2π/(γN ), we eliminate the latter frequency.
At the same time we add the frequency corresponding to the top K + 1th two-norms of the row of the estimated X. On the
other hand, we also adopt the afore-mentioned step to improve the estimation performance when the first and last element
of the sorted frequencies are equal to 0 and the last sampling frequency, respectively. We recycle the whole procedure until
there is no occurrence of those two incidents. Then we apply the cyclic refinement step of the Algorithm 1 to the estimated
frequencies, and we set Rc = 1, Rs = 1 in this case.
We adopt the SPA with a uniform linear array (ULA) method defined in [15] to implement sparse convex optimization. We
made one modification to avoid the frequency splitting phenomenon. Specifically, according to the [15], when optimized u∗
is obtained by solving a semidefinite programming problem, we then perform eigenvalue decomposition on T (u∗ ) and treat

11

the smallest eigenvalue as the noise variance denoted by σ̂ 2 . Here T (u∗ ) is a Toeplitz matrix constructed by u∗ . Finally, we
treat T (u∗ ) − diag(σ̂ 2 ) and the true model order K as the input parameters of the MATLAB rootmusic algorithm.
For the IRA algorithm [17], we generalize the single snapshot method to the scenario with MMVs. Considering the huge
computation cost, we set the oversampling rate to 2. Other system parameters are kept unchanged.
For the NOMP algorithm with MMVs, the corresponding partial parameters in various settings are summarized in Table II.
We don’t require explicit estimates of model order and set the overestimating probability to 0.01. All possible situations can be
divided into three cases: 1) The estimated model order is smaller than K. 2) We successfully detected K frequencies. In this
case, we directly calculate the NMSE.
3) The last case is that we overestimate the model order. If we encounter this situation,


we firstly sort the diag XH X and choose the K frequencies corresponding to the top K elements of the diag XH X as
the estimated frequencies. Then we take the same procedure to avoid frequency splitting as used in the Newtonized BPDN
algorithm.
Simulation Set-up: Here, we concentrate on estimating a mixture of K = 16 complex sinusoids with length N = 64 in the
uncompressed scenario. The number of snapshots is T = 10 and the noise variance is equal to 1. Each scenario, as defined
by ∆ωmin and SNR (19) is implemented over 300 MC trials. The original frequencies of the sinusoids in the mixture are
sampled from [0, 2π) and meet the corresponding minimum separation criterion. The detailed settings are shown in Table I.
In the following, we aim at obtaining the frequency estimation accuracy performance of the Algorithm 1 as compared to the
other methods in terms of error distribution and mean squared error.
TABLE I
S ETTINGS OF DIFFERENT SCENARIOS .
Scenarios
1
2
3
4

∆ωmin /∆ωDFT
2.5
1
2.5
1

SNR (dB)
SNRnom
SNRnom
Uniform[15, 35]
Uniform[15, 35]

TABLE II
PARAMETERS OF A LGORITHM 1 IN DIFFERENT SCENARIOS .
Algorithm 1
Rc
Rs

Scenario 1
1
1

Scenario 2
3
1

Scenario 3
1
1

Scenario 4
3
1

TABLE III
T IME [ SEC ] AVERAGED OVER 300 RUNS OF EACH ALGORITHM . PARAMETERS OF A LGORITHM 1 ARE SET BY II
Time [sec]
Scenario 1
Scenario 2
Scenario 3
Scenario 4

MUSIC
0.011
0.011
0.012
0.011

M-BPDN
0.069
0.066
0.060
0.064

NOMP
0.163
0.282
0.168
0.279

SPA
37.919
40.082
36.455
39.381

IRA
184.516
250.091
278.857
229.402

A. Error Distribution
Firstly, we plot the Complementary Cumulative Distribution Function (CCDF) of the squared frequency estimation error for
all algorithms to study the error distribution. In scenario 1 and 2, SNRnom is set as 25 dB for all sinusoids in the mixture.
As for scenario 3 and 4, the SNR corresponding to each sinusoid is selected independently and uniformly from [15, 35] dB.
It should be mentioned that all the algorithms are compared against the DFT method, which is implemented by coarse peak
peaking.
From Fig. 3, we can see that the IRA and the NOMP algorithm have a comparable performance, and both of them have the
highest estimation accuracy compared to the SPA and the M-BPDN algorithm. All those four methods outperform the MUSIC
algorithm in Scenario 1. The observation that the M-BPDN algorithm has a near-optimal performance shows that BPDN has
a very good estimation of the true frequencies. While in Scenario 2, the performance of the M-BPDN algorithm degrades
because the frequencies are becoming closer to each other. The MUSIC algorithm also suffers from the small separation
between frequencies.
The error distribution in Scenarios 3 and 4 are demonstrated in Fig. 4. In this scenario, we can still see that IRA and NOMP
algorithms have a very similar error distribution, but the gap between those two algorithms and other methods has increased.
In those two scenarios, SPA outperforms M-BPDN considerably compared to Scenario 1. When frequencies are closer, it can

12

be seen that the span of squared frequency estimation error has enlarged. The performance of all algorithms degrades in the
scenario where frequencies are harder to identify, but IRA and NOMP algorithms still achieve a much better performance.
From Table III, we can conclude that the MUSIC algorithm is the fastest algorithm among all these algorithms. The NOMP
algorithm has a better performance but a slower speed compared to the M-BPDN algorithm. The computation costs of SPA
and IRA are large, where the latter method is the slowest but almost always achieves the best estimation accuracy.
0

10

CCDF

DFT (γ = 4)
NOMP (Rc = 1)
SPA
M−BPDN (γ = 4)
MUSIC
IRA (γ = 2)

−1

10

0.6
0.4
0.2
−22 −20 −18 −16
−2
10
−50
−40
−30
−20
−10
0
Squared frequency estimation error in dB (relative to DFT)

(a) Scenario 1
0

10

CCDF

DFT (γ = 4)
NOMP (Rc = 3)
SPA
M−BPDN (γ = 4)
MUSIC
IRA (γ = 2)

−1

10

0.6
0.4
0.2
−2
−24−22−20−18−16−14
10
−50
−40
−30
−20
−10
0
Squared frequency estimation error in dB (relative to DFT)

(b) Scenario 2

Fig. 3. CCDF of the frequency MSE for Scenario 1 and 2.

B. Mean squared error
In this part, we aim to compare the estimation performance of all algorithms by using NMSE and CRB in all four scenarios.
Our goal is to obtain the estimation performance of a single sinusoid (ignoring the other sinusoids). Considering the huge
computation cost of IRA algorithm, we decide not to consider it in this section.
In scenario 1 and 2, SNRnom takes values from 5 dB to 35 dB. We regard
 the mean value
 of CRB of all K sinusoids
as the CRB corresponding to a single sinusoid. The NMSE is defined as E kωest − ωtrue k22 /∆2DFT . For scenarios 3 and
4, we keep the SNR of one of the sinusoids at a given value taken from [5, 35] dB, while the other (K − 1) SNRs are
drawn independently and uniformly from [15, 35] dB. In these two scenarios, the NMSE is |ωest − ωtrue |2 /∆2DFT , where
ωtrue denotes the sinusoid according to the fixed SNR. We choose the CRB that corresponds to this sinusoid to plot along
with NMSE.
Fig. 5 shows the MSE of frequency estimation in Scenario 1 and 2. We can see that NOMP achieves the CRB when SNR
is equal to 9 dB and it benefits from increasing the SNR in Scenario 1. SPA almost achieves the CRB at a smaller SNR value
compared to the NOMP method and likewise achieves a better estimation performance as SNR increases. On the other hand,
the M-BPDN algorithm achieves performance floors when SNR is bigger than 29 dB. Although the MUSIC algorithm follows
the bound when SNR is bigger than 15 dB, it doesn’t achieve the CRB. In scenario 2, we can see that the NOMP and SPA

13

0

10

CCDF

DFT (γ = 4)
NOMP (Rc = 1)
SPA
M−BPDN (γ = 4)
MUSIC
IRA (γ = 2)

−1

10

0.6
0.4
0.2
−45 −40 −35 −30
−2
10
−60
−50
−40
−30
−20
−10
0
Squared frequency estimation error in dB (relative to DFT)

(a) Scenario 3
0

10

DFT (γ = 4)
NOMP (Rc = 3)

CCDF

SPA
M−BPDN (γ = 4)
MUSIC
IRA (γ = 2)
−1

10

0.6
0.4
0.2
−2

−40

10
−60

−30

−20

−50
−40
−30
−20
−10
Squared frequency estimation error in dB (relative to DFT)

0

(b) Scenario 4

Fig. 4. CCDF of the frequency MSE for Scenario 3 and 4.

algorithm almost equal to CRB as SNR increases. The NOMP algorithm slightly deviates from CRB when SNR is larger than
25 dB. The M-BPDN method reaches performance floors at a smaller SNR value compared to Scenario 1. MUSIC benefits
from increasing the SNR and the gap between MUSIC and CRB is larger than Scenario 1.
Fig. 6 (a) corresponds to Scenario 3. It shows that NOMP achieves CRB earlier than any other method and doesn’t reach
performance floors. SPA benefits from increasing the SNR and becomes closer to CRB. The M-BPDN algorithm reaches
performance floors at a large SNR values and has less accurate performance compared to MUSIC. The same observations can
also be applicable to Scenario 4, where the separation between frequencies is small.
VII. C ONCLUSION
The NOMP algorithm is extended to deal with MMVs settings. The algorithm uses the information of already detected
frequencies to refine the current frequency and has a stopping criterion based on a given overestimating probability. The
convergence results are provided. Compared to the classical MUSIC algorithm, more recent SPA and M-BPDN algorithm, the
NOMP algorithm has a superior performance in terms of frequency estimation accuracy and asymptotically approaches the
CRB. As for the run time, the NOMP algorithm is much faster than the SPA and IRA method.
VIII. APPENDIX: E STIMATION THEORETIC BOUNDS
Given the observation vector y which determines the estimator θ̂(y), a widely used criterion to characterize the performance
of the estimator is the mean square error (MSE), given by
mse(θ̂(y)) = Ey;θ [(θ̂(y) − θ)H (θ̂(y) − θ)].

(55)

Here, Ey;θ [·] denotes the expectation taken with respect to the probability density function (pdf) p(y; θ) of the measurement
vector y parameterized by θ. A commonly employed approach seeks unbiased estimators satisfying Ey;θ [θ̂(y)] = θ which

14

10
CRB
NOMP
SPA
MUSIC
M−BPDN

Normalized MSE

0
−10
−20
−30
−40
−50
−60
5

10

15

20
25
SNRnom [dB]

30

35

(a) Scenario 1
10
CRB
NOMP
SPA
MUSIC
M−BPDN

Normalized MSE

0
−10
−20
−30
−40
−50
−60
5

10

15

20
25
SNRnom [dB]

30

35

(b) Scenario 2

Fig. 5. Normalized frequency MSE for Scenarios 1 and 2.
TABLE IV
T HE PARTIAL DERIVATIVES OF s(x) WITH RESPECT TO x. N OTE THAT ek IS A STANDARD BASIS VECTOR WITH A k TH COMPONENT OF ONE AND ALL THE
OTHER COMPONENTS BEING ZERO

k&l\∂(·)/∂(·)

∂s
∂xk

∂s
∂xl

K + 1 ≤ k ≤ K + KT, K ≤ l ≤ K + KT
K + 1 ≤ k ≤ K + KT, K + KT + 1 ≤ l ≤ K + 2KT
K + KT + 1 ≤ k ≤ K + 2KT, K ≤ l ≤ K + KT
K + KT + 1 ≤ k ≤ K + 2KT, K + KT + 1 ≤ l ≤ K + 2KT

(IT ⊗ A) ek−K
(IT ⊗ A) ek−K
0N T
0N T

(IT ⊗ A) el−K
0N T
(IT ⊗ A) el−K
0N T

∂s
∂x∗
l

∂s
∂x∗
k

0N T
(IT ⊗ A) el−K−KT
0N T
(IT ⊗ A) el−K−KT

0N T
0N T
(IT ⊗ A) ek−K−KT
(IT ⊗ A) ek−K−KT

minimize MSE. Fortunately, the performance of an unbiased estimator can be characterised by the CRB. In the AWGN scenario,
to evaluate the CRB for model (1) in the case of Φ = I, the following lemma is utilized.
Lemma 2 [27] (General Complex Gaussian CRLB) Assume the observations y follow the complex Gaussian distribution
y ∼ CN (s(θ), Γ(θ)), then the (k, l)th element of the complex FIM J(θ) is

H

 
H







∂s
∂s
∂s
∂s
∂Γ
∂Γ
−1
−1
−1
−1
[J(θ)]kl =
Γ (θ)
+
Γ (θ)
+ tr
Γ (θ)
Γ (θ) ,
(56)
∂θk
∂θl
∂θ∗l
∂θ∗k
∂θ∗k
∂θl
where θ = [θ T , θ H ]T ∈ C2p .
By vectorizing both sides of model (1), we can obtain
y = (IT ⊗ A) x + z,

(57)

15

40
CRB
NOMP
SPA
MUSIC
M−BPDN

Normalized MSE

20

0

−20

−40

−60
5

10

15

20
25
SNRnom [dB]

30

35

(a) Scenario 3
40
CRB
NOMP
SPA
MUSIC
M−BPDN

Normalized MSE

20

0

−20

−40

−60
5

10

15

20
25
SNRnom [dB]

30

35

(b) Scenario 4

Fig. 6. Normalized frequency MSE for Scenarios 3 and 4.

T


T
T
where y = vec (Y) , x = vec (X) , z = vec (Z). In this case, θ = ω T , xT , x = xT , xH and ω = [ω1 , ω2 , · · · , ωK ] .
We can conclude that s(θ) = (IT ⊗ A) x , Γ(θ) = σ 2 IN T . We divide the calculation of FIM into four parts,


J11 J12
J=
.
(58)
J21 J22
In the following, we calculate the submatrices of J separately.
A. J11
The (k, l)th element of J is
(

H

 
H


 H 
 )
∂s
∂s
∂s
∂s
2
∂A
∂A
Jk,l =
Γ−1
+
Γ−1
= 2R
IT ⊗
x
IT ⊗
x
∂ωk
∂ωl
∂ωl
∂ωk
σ
∂ωk
∂ωl
( "
( "
# )

H # 
 )

H
∂A
∂A
∂A
∂A
(a) 2
(b) 2
H
H
= 2 R x IT ⊗
IT ⊗
x = 2 R x IT ⊗
x ,
σ
∂ωk
∂ωl
σ
∂ωk
∂ωl
H

(59)

for 1 ≤ k, l ≤ K, where equality (a) follows from (A ⊗ B) = AH ⊗BH , and (b) is a direct consequence of (A ⊗ B) (C ⊗ D) =
(AC) ⊗ (BD). According to



tr (ABCD) = vecT (B) C ⊗ AT vec DT ,
(60)

16

we can further write Jk,l as
Jk,l

2
= 2R
σ

( "
( "
#)
#)
T 
∗
H
∂A
∂A
∂A
∂A
(a) 2
∗ T
H
tr
= 2 R tr
X X
XX
,
∂ωl
∂ωk
σ
∂ωl
∂ωk

(61)

where (a) follows from R(a + bj) = R(a − bj). By defining
bl =

∂a(ωl )
,
∂ωl

(62)

Ψ , XXH ∈ CK×K ,

(63)

we obtain


∂A
∂ωl

H



H 

∂A
∂a(ωl )
∂a(ωk )

H
= 00
00
00
0 0 = [0 0 bl 0 0] [0 0 bk 0 0] = 
∂ωk
∂ωl
∂ωk

0
..
.
0

···
bH
l bk
···


0
.. 
,
. 
0 K×K

(64)

Thus
Jk,l =

2
R{bH
l bk Ψkl }, 1 ≤ k, l ≤ K.
σ2

(65)

B. J22
The (k, l)th element of J is
J22 =

1
σ2



IT ⊗ AH A
0KT ×KT



0KT ×KT 
IT ⊗ A H A


=


1
I2T ⊗ AH A .
2
σ

(66)

for K + 1 ≤ k, l ≤ K + 2KT .
C. J12
The (k, l)th element of J is


 H
 1
∂A

 
H


H

 
I
⊗
x (IT ⊗ A) el−K , (a)

T
∂s
∂s
∂s
∂s
σ2
∂ωk
−1
+
Γ
Jk,l =
Γ−1
=



∂ωk
∂xl
∂x∗l
∂ωk
∂A
1

H

IT ⊗
x, (b)
 2 [(IT ⊗ A) el−K−KT ]
σ
∂ωk


∗



1
∂A
T

T
1

X∗ [reshape(el−K , K, T )]
, (a)
∗ ∗


 σ 2 tr A
 2 aT
∂ωk
i1 bk Xkj1 , (a)
σ
(
)

T
=
=


∂A
1

 1 bT a∗ Xj k , (b).
∗
∗
T

tr
A
reshape((e
)
,
K,
T
)X
,
(b).
l−K−KT
 2
σ 2 k i2 2
σ
∂ωk
(67)
where (a) is 1 ≤ k ≤ K, K + 1 ≤ l ≤ K + KT , (b) is 1 ≤ k ≤ K, K + KT + 1 ≤ l ≤ K + 2KT , i1 = b(l − K − 1)/T c + 1,
j1 = l − K − T (b(l − K − 1)/T c), i2 = l − K − KT − K(b(l − K − KT − 1)/Kc), j2 = b(l − K − KT − 1)/Kc + 1.
D. J21
The (k, l)th element of J is



∂A
1

H

[(I
⊗
A)
e
]
I
⊗
x, (c)

k−K
T
 σ2 T
∂ωl
∂s
∂s
∂s
−1 ∂s
−1
Jk,l =
Γ
+
Γ
=

 H

∂xk
ωl
∂ωl
∂x∗k
1
∂A


IT ⊗
x (IT ⊗ A) ek−K−KT , (d),
 2
σ
∂ωl
(
)

T

1
∂A

1

∗


tr
A∗ reshape((ek−K )∗ , K, T )XT , (c)
 σ2
 2 bT
l ai3 Xj3 l , (c)
∂ωl
σ
=
=


∗




 1 aT b∗ X∗ , (d)
1
∂A
T

T
∗
 tr A
X
[reshape(e
,
K,
T
)]
,
(d).
k−K−KT
σ 2 i4 l lj4
σ2
∂ω


H







H





l

(68)
where (c) is 1 ≤ l ≤ K, K + 1 ≤ k ≤ K + KT , (d) is 1 ≤ l ≤ K, K + KT + 1 ≤ k ≤ K + 2KT , i3 = k − K − K(b(k −
K − 1)/Kc), j3 = b(k − K − 1)/Kc + 1, i4 = b(k − K − KT − 1)/T c + 1, j4 = k − K − KT − T (b(k − K − KT − 1)/T c).

17

R EFERENCES
[1] H. V. Poor, An introduction to signal detection and estimation, Springer-Verlag New York, 1994.
[2] S. M. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory, 1993, Englewood Cliffs, NJ:
Prentice Hall.
[3] R. Schmidt, “Multiple emitter location and signal parameter estimation,” IEEE Trans. Antennas Propag., vol. 34, no. 3,
pp. 276-280, 1986.
[4] R. Roy and T. Kailath, “ESPRIT - estimation of signal parameters via rotational invariance techniques,” IEEE Trans.
Acoust., Speech, Signal Process., vol. 37, no. 7, pp. 984-995, 1989.
[5] H. Teutsch, Modal Array Signal Processing: Principles and Applications of Acoustic Wavefield Decomposition, ser. Lecture
Notes in Control and Information Sciences. Berlin, Germany: Springer, 2007
[6] D. Malioutov, M. Cetin ,and A. Willsky, “A sparse signal reconstruction perspective for source localization with sensor
arrays,” IEEE Trans. Signal Process., vol. 53, no. 8, pp. 3010-2022, 2005.
[7] M. Hyder and K. Mahata, “Direction-of-arrival estimation using a mixed `2,0 norm approximation,” IEEE Trans. Signal
Process., vol. 58, no. 9, pp. 4646-4655, 2010.
[8] E. Candes and J. Romberg, “Sparsity and incoherence in compressive sampling, ”Inverse problems, vol. 23, no. 3, pp. 969,
2007.
[9] Y. Chi, L. L. Scharf, A. Pezeshki, and R. Calderbank, “Sensitivity to basis mismatch in compressed sensing,” IEEE Trans.
Signal Process., vol. 59, no. 5, pp. 2182-2195, 2011.
[10] L. L. Scharf, E. K. Chong, A. Pezeshki, and J. R. Luo, “Sensitivity considerations in compressed sensing,” in Proceedings
of the Asilomar Conference on Signals, Systems and Computers. IEEE, pp. 744-748, 2011.
[11] P. Pakrooh, L. L. Scharf, A. Pezeshki, and Y. Chi, “Analysis of fisher information and the Cramer-Rao bound for nonlinear
parameter estimation after compressed sensing,” IEEE Trans. Signal Process., vol. 63, no. 23, pp. 6423-6428, Dec. 2015.
[12] P. Stoica, P. Babu, and J. Li, “New method of sparse parameter estimation in separable models and its use for spectral
analysis of irregularly sampled data,” IEEE Trans. Signal Process., vol. 59, no. 1, pp. 35-47, 2011.
[13] P. Stoica, P. Babu, and J. Li, “SPICE: A sparse covariance-based estimation method for array processing,” IEEE Trans.
Signal Process., vol. 59, no. 2, pp. 629-638, 2011.
[14] P. Stoica and P. Babu, “SPICE and LIKES: Two hyperparameter-free methods for sparse-parameter estimation”, Signal
Process., vol. 92, no. 7, pp. 1580-1590, 2012.
[15] Z. Yang, L. Xie, and C. Zhang, “A discretization-free sparse and parametric approach for linear array signal processing,”
IEEE Trans. Signal Process., vol. 62, no. 19, pp. 4959-4973, 2014.
[16] Y. Li and Y. Chi, “Off-the-grid line spectrum denoising and estimation with multiple measurement vectors,” IEEE Trans.
Signal Process., vol. 64, no. 5, pp. 1257-1269, 2015.
[17] J. Fang, F. Wang, Y. Shen, H. Li ,and R. S. Blum, “Super-resolution compressed sensing for line spectral estimation: an
iterative reweighted approach,” IEEE Trans. Signal Process., vol. 64, no. 18, pp. 4649-4662, 2016.
[18] B. Mamandipoor, D. Ramasamy, and U. Madhow, “Newtonized orthogonal matching pursuit: Frequency estimation over
the continuum,” IEEE Trans. Signal Process., vol. 64, no. 19, pp. 5066-5081, 2016.
[19] A. Hjφrungnes, Complex-Valued Matrix Derivatives: With Applications in Signal Processing and Communications,
Cambridge University Press, 2011.
[20] A. Schaeffer, “Inequalities of a. markoff and s. bernstein for polynomials and related functions,” Bull. Amer. Math. Soc,
vol. 47, pp. 565-579, 1941.
[21] B. N. Bhaskar, G. Tang, and B. Recht, “Atomic norm denoising with applications to line spectral estimation,” IEEE Trans.
Signal Process., vol. 61, no. 23, pp. 5987-5999, 2013.
[22] A. R. Barron, A. Cohen, W. Dahmen, and R. A. Devore, “Approximation and learning by greedy algorithms,” Ann. Statist.,
pp. 64-94, 2008.
[23] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE
Trans. Inf. Theory, vol. 53, no. 12, pp. 4655-4666, 2007.
[24] E. V. D. Berg and M. P. Friedlander, “Sparse optimization with least-squares constraints,” SIAM J. Optim, vol 21, no.4,
pp. 1201-1229, 2011.
[25] Z. Yang and L. Xie, “On gridless sparse methods for line spectral estimation from complete and incomplete data,” IEEE
Trans. Signal Process., vol. 63, no. 12, pp. 3139-3153, 2015.
[26] Z. Yang and L. Xie, “Exact joint sparse frequency recovery via optimization methods,” IEEE Trans. Signal Process., vol.
64, no. 19, pp. 5145-5157, Oct 1, 2016.
[27] V. Nagesha and S. M. Kay, “Cramer-Rao lower bounds for complex parameters,” available at www.ele.uri.edu/faculty/
kay/New%20web/downloadable%20files/Nageha complex%20CRLB.pdf.

