Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
xx.xxxx

This work was supported in part by the Japan Society for the Promotion of
Science (JSPS) KAKENHI under Grant JP21K17827 and JP21H01592.

Corresponding author: Hiroki Kuroda (e-mail: kuroda@vos.nagaokaut.ac.jp).

# Convex Estimation of Sparse-Smooth Power Spectral Densities from Mixtures of
Realizations with Application to Weather Radar

HIROKI KURODA1  DAICHI KITAHARA2  EIICHI YOSHIKAWA3, 4  HIROSHI KIKUCHI5  and
TOMOO USHIO2  Department of Information and Management Systems Engineering,
Nagaoka University of Technology, Niigata 940-2188, Japan  Division of
Electrical, Electronic and Infocommunications Engineering, Osaka University,
Osaka 565-0871, Japan  Aeronautical Technology Directorate, Japan Aerospace
Exploration Agency, Tokyo 181-0015, Japan  Department of Electrical and
Computer Engineering, Colorado State University, Fort Collins 80521 CO, USA
Center for Space Science and Radio Engineering, The University of Electro-
Communications, Tokyo 182-8585, Japan

###### Abstract

In this paper, we propose a convex optimization-based estimation of sparse and
smooth power spectral densities (PSDs) of complex-valued random processes from
mixtures of realizations. While the PSDs are related to the magnitude of the
frequency components of the realizations, it has been a major challenge to
exploit the smoothness of the PSDs because penalizing the difference of the
magnitude of the frequency components results in a nonconvex optimization
problem that is difficult to solve. To address this challenge, we design the
proposed model that jointly estimates the complex-valued frequency components
and the nonnegative PSDs, which are respectively regularized to be sparse and
sparse-smooth. By penalizing the difference of the nonnegative variable that
estimates the PSDs, the proposed model can enhance the smoothness of the PSDs
via convex optimization. Numerical experiments on the phased array weather
radar, an advanced weather radar system, demonstrate that the proposed model
achieves superior estimation accuracy to the existing sparse estimation models
combined with or without post-smoothing.

###### Index Terms:

Power spectral density estimation, random process, sparsity, smoothness,
regularization, convex optimization, weather radar.

##  I Introduction

Power spectral density (PSD) of a random process describes how power of the
random process is distributed over frequency. Estimation of the PSD from
realizations of a random process is a fundamental problem in science and
engineering [1, 2, 3]. For weather radar applications, the PSD estimation is
essential for analysis of weather phenomena because the PSD contains
information of the precipitation intensity and the Doppler velocity
distribution [4, 5, 6, 7, 8]. For example, the parabolic Doppler weather radar
[4] transmits a pencil beam and observes backscattered signals in a narrow
range of elevation angles, which can be regarded as realizations of a single
random process whose PSD reflects the weather condition in the narrow range.

We consider the estimation of PSDs from mixtures of realizations of random
processes, which is much more challenging than the classical case of single
random process. Our primary interest is in the phased array weather radar
(PAWR) [9, 10, 11, 12, 13, 14, 15, 16], which is developed to detect hazardous
weather phenomena. The Doppler weather radar is not capable of detection of
hazardous weather because of its mechanical vertical scan to observe
backscattered signals in multiple elevation angles, which requires a long
observation time. To shorten the observation time, the PAWR transmits a fan
beam and observes backscattered signals in a wide range of elevation angles.
The backscattered signals observed by the PAWR can be modeled as mixtures of
realizations of random processes whose PSDs reflect the weather conditions in
finely-divided ranges. Thus, to obtain the weather condition in a fine
resolution, the PAWR needs digital signal processing which recovers the PSDs
in the finely-divided ranges.

Since the estimation of PSDs from mixtures of realizations is a challenging
problem, the major existing methods employ a two-step approach that first
estimates the frequency components of the realizations and then estimates the
PSDs. For the frequency component estimation, sparsity-aware methods have
achieved significant improvements on estimation accuracy over the classical
linear ones in many fields [17, 18, 19, 20, 21, 22, 23, 16]. Sparsity in the
spatial domain is exploited in [17, 18, 19] by assuming that signals arrive
from a few angles. Unfortunately, this assumption is not suitable for the PAWR
because targets such as clouds and raindrops exist at many angles [16]. In
[20, 21, 22, 23], isolated sparse frequency components, called _line spectra_
, are estimated based on the $\ell_{1}$ regularization. It is demonstrated in
[16] that a block-sparse regularization model using the mixed
$\ell_{2}/\ell_{1}$ norm [24, 25, 26, 27, 28, 29, 30] is more effective for
weather radar applications because the frequency components are clustered in a
few blocks due to the narrow-bandwidth of the PSDs. After the frequency
component estimation, the _periodogram_ , i.e., the squared magnitude of the
frequency components, is usually employed to estimate the PSD because of its
asymptotic unbiasedness. However, the periodogram has a drawback of large
variance and erratic oscillation [1, 2, 3, 4]. For the classical case of
single random process, smoothing techniques, e.g., those shown in [1, 2], are
often used to reduce the variance and the erratic oscillation. While the
existing smoothing techniques can be used as a post-processing, such a two-
stage approach would be sub-optimal because the smoothness is not exploited in
the frequency component estimation. Since the PSD is estimated by the
periodogram that is the squared magnitude of the frequency components, one may
add the penalty for the difference between the magnitude of the frequency
components in the frequency component estimation. However, due to the
nonconvexity of this type of penalty (see, e.g., [31]), it is hard for this
approach to obtain an optimal solution, and the performance dependency on the
initial estimate and the optimization algorithm is difficult to elude.

Another line of works derive an approximated observation model between the
realizations and the PSD for the cases of, e.g., single random process [32,
33, 34, 35] and spatially independent random processes [36, 37]. Since the
approximated observation model is written in terms of the (nonnegative) PSD,
the smoothness of the PSD could be exploited via convex optimization. However,
this approach takes the magnitude of the observation model to derive the
approximated observation models, which implies that half information of the
original observation model is lost as the phase information is discarded. In
particular, this approach is not applicable to the PAWR because signals from
different angles cannot be distinguished by the magnitude information (see (6)
in Example 1).

In this paper, we propose a convex optimization-based method that
simultaneously estimates block-sparse frequency components and block-sparse
and smooth PSDs from mixtures of realizations.111If a target is sparse and
smooth, the target is also block-sparse since nonzero components are clustered
due to the smoothness. To design the proposed method, we first apply the
optimally structured block-sparse model of [38] for the frequency component
estimation. Then, we newly leverage the latent variable of the designed model,
which is originally introduced to optimize the block structure, for the PSD
estimation. More precisely, we demonstrate that the latent variable is in fact
related to the square root of the PSDs, which enables to exploit the
smoothness of the PSDs via convex optimization. The main contributions are
summarized as follows.

  * •

We present, for the first time in the literature, a convex optimization-based
method that can exploit the smoothness of the PSDs for the estimation from
mixtures of realizations.

  * •

We show that many smoothness priors designed for real-valued signals, e.g.,
the high-order total variation [39, 40] and the total generalized variation
[41], can be directly incorporated into the proposed framework to enhance the
smoothness of the PSDs of complex-valued random processes thanks to the
nonnegative latent variable.

  * •

We conduct thorough numerical simulations on the PAWR, which demonstrate that
the proposed method achieves superior estimation accuracy to the existing
sparse estimation models with or without post-smoothing.

The rest of this paper is organized as follows. In Section II, we formulate
the estimation of PSDs from mixtures of realizations of random processes, and
clarify its relation to weather radar applications. In Section III, we design
the proposed convex optimization model that simultaneously estimates block-
sparse frequency components and block-sparse and smooth PSDs. Section IV
presents numerical experiments on the PAWR, followed by conclusion in Section
V.

A preliminary short version of this paper was presented at a conference [42].

Notations: $\mathbb{N}$, $\mathbb{R}$, $\mathbb{R}_{+}$, $\mathbb{R}_{++}$,
and $\mathbb{C}$ respectively denote the sets of nonnegative integers, real
numbers, nonnegative real numbers, positive real numbers, and complex numbers.
We use $\imath\in\mathbb{C}$ to denote the imaginary unit, i.e.,
$\imath=\sqrt{-1}$. For every $x\in\mathbb{C}$, $x^{*}$ denotes the complex
conjugate of $x$, and $|x|:=\sqrt{x^{*}x}$ denotes the absolute value of $x$.
For matrices or vectors, we denote the transpose and the Hermitian transpose
respectively by $(\cdot)^{\top}$ and $(\cdot)^{\mathrm{H}}$. The identity
matrix of order $N$ is denoted by $\mathbf{I}_{N}\in\mathbb{R}^{N\times N}$.
We denote the diagonal matrix with components of $\mathbf{w}\in\mathbb{C}^{N}$
on the main diagonal by $\mathrm{diag}(\mathbf{w})\in\mathbb{C}^{N\times N}$.
The cardinality of a set $\mathcal{A}$ is denoted by $|\mathcal{A}|$. The
$\ell_{2}$ (Euclidean) norm, the $\ell_{1}$ norm, and the $\ell_{0}$ pseudo-
norm of $\mathbf{x}=(x_{1},\ldots,x_{N})^{\top}\in\mathbb{C}^{N}$ are
respectively denoted by $\|\mathbf{x}\|:=\sqrt{\sum_{n=1}^{N}|x_{n}|^{2}}$,
$\|\mathbf{x}\|_{1}:=\sum_{n=1}^{N}|x_{n}|$, and
$\|\mathbf{x}\|_{0}:=\left|\\{n\in\\{1,\ldots,N\\}\,|\,x_{n}\neq 0\\}\right|$.
The expectation operator is denoted by $E[\cdot]$.

##  II Problem Formulation

We consider the estimation problem of power spectral densities (PSDs) of $N$
random processes from noisy mixtures of their realizations. We denote the
$n$-th complex-valued discrete-time random process $(n=1,\ldots,N)$ by

$\displaystyle X_{n}^{\star}[\ell]\in\mathbb{C}\qquad(\ell=0,\pm 1,\pm
2,\ldots).$ (1)

To define the PSD of $X_{n}^{\star}[\ell]$, we assume that
$X_{n}^{\star}[\ell]$ is zero-mean and wide-sense stationary, which imply that
$E[X_{n}^{\star}[\ell]]=0$ for any $\ell$, and the auto-correlation

$\displaystyle E[X_{n}^{\star}[m+\ell](X_{n}^{\star}[m])^{*}]=:R_{n}[\ell]$

does not depend on $m$ for any $\ell$. In addition, suppose
$\sum_{\ell=-\infty}^{\infty}|R_{n}[\ell]|<\infty$. Then, the PSD of
$X_{n}^{\star}[\ell]$ is given by

$\displaystyle
S^{\star}_{n}(f):=\sum_{\ell=-\infty}^{\infty}R_{n}[\ell]e^{-\imath 2\pi
f\ell}\quad\left(f\in\left[-\frac{1}{2},\frac{1}{2}\right)\right).$ (2)

We denote $L$ consecutive realizations of $X_{n}^{\star}[\ell]$ by

$\displaystyle\bar{x}_{j,n}[\ell]\in\mathbb{C}\qquad(\ell=1,\ldots,L),$ (3)

where $j\in\\{1,\ldots,J\\}$ is the index of trials. Note that
$\bar{x}_{j,n}[\ell]$ for $j=1,\ldots,J$ are assumed to be realizations of the
common random process $X_{n}^{\star}[\ell]$ (see Remark 1 for validity of this
assumption). We define the observation model by

$\displaystyle\mathbf{y}_{j}:=\sum_{n=1}^{N}\mathbf{A}_{n}\bar{\mathbf{x}}_{j,n}+\boldsymbol{\varepsilon}_{j}\in\mathbb{C}^{d}\quad(j=1,\ldots,J),$
(4)

where the realizations in (3) are collectively denoted by

$\displaystyle\bar{\mathbf{x}}_{j,n}$
$\displaystyle:=(\bar{x}_{j,n}[1],\bar{x}_{j,n}[2],\ldots,\bar{x}_{j,n}[L])^{\top}\in\mathbb{C}^{L},$
(5)

$\mathbf{A}_{n}\in\mathbb{C}^{d\times L}$ is the known matrix that models the
observation process for the $n$-th source, and
$\boldsymbol{\varepsilon}_{j}\in\mathbb{C}^{d}$ is the (unknown) observation
noise. Our goal is to estimate the PSDs $S^{\star}_{n}(f)\quad(n=1,\ldots,N)$
from $\mathbf{y}_{j}\quad(j=1,\ldots,J)$ and
$\mathbf{A}_{n}\quad(n=1,\ldots,N)$ in (4). Note that the classical PSD
estimation problem for a single random process [1, 2, 3], e.g., for the
Doppler weather radar [4], is a special instance of (4) for $N=1$, $d=L$, and
$\mathbf{A}_{n}=\mathbf{I}_{L}$. The generalized observation model (4) is
introduced to handle the PSD estimation for the phased array weather radar
(PAWR) [9, 10, 11, 12, 13, 14, 15, 16], which is our primary interest.

######  Example 1 (Phased array weather radar (PAWR)).

For the PAWR, $X_{n}^{\star}[\ell]$ corresponds to the sum of backscattered
signals in the angular interval
$\left[\theta_{n}-\frac{\Delta\theta}{2},\theta_{n}+\frac{\Delta\theta}{2}\right]$,
where $\theta_{n}\quad(n=1,\ldots,N)$ are the equally spaced angles with the
spacing $\Delta\theta$. By using an $M$-element uniform linear array, the PAWR
observes noisy mixtures of realizations by

$\displaystyle\mathbf{y}_{j}[\ell]:=\sum_{n=1}^{N}\mathbf{a}(\theta_{n})\bar{x}_{j,n}[\ell]+\boldsymbol{\varepsilon}_{j}[\ell]\in\mathbb{C}^{M}\quad(\ell=1,\ldots,L)$
(6)

for each $j\in\\{1,\ldots,J\\}$, where
$\mathbf{a}(\theta_{n})\in\mathbb{C}^{M}$ is the known steering vector for the
angle $\theta_{n}$, and $\boldsymbol{\varepsilon}_{j}[\ell]\in\mathbb{C}^{M}$
is the white Gaussian noise. More precisely, $\mathbf{a}(\theta)$ is defined
by

$\displaystyle\mathbf{a}(\theta):=\left(1,e^{-\imath\frac{2\pi\Delta\sin\theta}{\lambda_{\mathrm{cw}}}},\ldots,e^{-\imath\frac{2(M-1)\pi\Delta\sin\theta}{\lambda_{\mathrm{cw}}}}\right)^{\top}\in\mathbb{C}^{M},$

where $\Delta$ is the inter-element spacing of the uniform linear array, and
$\lambda_{\mathrm{cw}}$ is the carrier wavelength. The observation model (6)
for the PAWR can be written in the form of (4) by setting

$\displaystyle\mathbf{y}_{j}$
$\displaystyle=(\mathbf{y}_{j}[1]^{\top},\mathbf{y}_{j}[2]^{\top},\ldots,\mathbf{y}_{j}[L]^{\top})^{\top}\in\mathbb{C}^{ML},$
$\displaystyle\boldsymbol{\varepsilon}_{j}$
$\displaystyle=(\boldsymbol{\varepsilon}_{j}[1]^{\top},\boldsymbol{\varepsilon}_{j}[2]^{\top},\ldots,\boldsymbol{\varepsilon}_{j}[L]^{\top})^{\top}\in\mathbb{C}^{ML},$

and $\mathbf{A}_{n}\in\mathbb{C}^{ML\times L}$ to the block-diagonal matrix
that contains $L$ copies of $\mathbf{a}(\theta_{n})$ on the diagonal blocks,
i.e.,

$\displaystyle\mathbf{A}_{n}=$
$\displaystyle\begin{pmatrix}\mathbf{a}(\theta_{n})&&&\\\
&\mathbf{a}(\theta_{n})&&\\\ &&\ddots&\\\
&&&\mathbf{a}(\theta_{n})\end{pmatrix}.$

The estimation of the PSDs $S^{\star}_{n}(f)$ of
$X_{n}^{\star}[\ell]\quad(n=1,\ldots,N)$ is essential for the PAWR because
$S^{\star}_{n}(f)$ contains information about the weather condition, i.e., the
precipitation intensity and the Doppler velocity distribution222For weather
radar applications, even when aliasing occurs, the continuous-time PSD can be
recovered from the discrete-time PSD in (2) unless the variation of Doppler
velocity is extremely large (see, e.g., [4, Chapter 5] and [16, Section
II-A])., in the narrow angular interval
$\left[\theta_{n}-\frac{\Delta\theta}{2},\theta_{n}+\frac{\Delta\theta}{2}\right]$.

######  Remark 1 (Tradeoff between $L$ and $J$).

To derive the observation model in (4) where
$\bar{x}_{j,n}[\ell]\quad(j=1,\ldots,J)$ are realizations of the common random
process $X_{n}^{\star}[\ell]$, similarly to the case of a single random
process $(N=1)$ [1, 2, 3, 4], we split the whole observations into $J$
subsets. For weather radar applications, $JL$ corresponds to the observation
time, and thus is set small enough so that the statistics of the targets such
as clouds and raindrops are (approximately) unchanged. Since $L$ defines the
frequency resolution (see (11)), $L$ is desired to be large enough. Thus, $J$
is typically set very small, and $J=1$ is of particular interest in weather
radar applications.

We rewrite (4) to an observation model in terms of frequency components of the
time-domain realizations $\bar{\mathbf{x}}_{j,n}$ in (5) because they are
related to the PSD in (2) more directly than $\bar{\mathbf{x}}_{j,n}$. More
precisely, we represent $\bar{\mathbf{x}}_{j,n}$ as

$\displaystyle\bar{\mathbf{x}}_{j,n}=\mathbf{G}\bar{\mathbf{u}}_{j,n}$ (7)

for $j=1,\ldots,J$ and $n=1,\ldots,N$, where

$\displaystyle\bar{\mathbf{u}}_{j,n}:=(\bar{u}_{j,n}[1],\bar{u}_{j,n}[2],\ldots,\bar{u}_{j,n}[L])^{\top}\in\mathbb{C}^{L}$
(8)

is used as the frequency components, and $\mathbf{G}\in\mathbb{C}^{L\times L}$
is a suitable synthesis matrix. Substituting the representation (7) to the
observation model (4), we have

$\displaystyle\mathbf{y}_{j}=\sum_{n=1}^{N}\mathbf{A}_{n}\mathbf{G}\bar{\mathbf{u}}_{j,n}+\boldsymbol{\varepsilon}_{j}\in\mathbb{C}^{d}\quad(j=1,\ldots,J),$
(9)

which is used as the observation model for the frequency components
$\bar{\mathbf{u}}_{j,n}$. The representation (7) covers popular frequency
analysis methods used in weather radar applications, e.g., the discrete
Fourier transform (DFT) and the windowed DFT.

######  Example 2 (DFT).

Let $\bar{\mathbf{u}}^{\mathrm{(DFT)}}_{j,n}$ be the normalized DFT
coefficients of $\bar{\mathbf{x}}_{j,n}$. Define the normalized DFT matrix
$\mathbf{F}\in\mathbb{C}^{L\times L}$ by

$\displaystyle\mathbf{F}:=\frac{1}{\sqrt{L}}\begin{pmatrix}1&e^{-\imath 2\pi
f_{1}}&\cdots&e^{-\imath 2\pi f_{1}(L-1)}\\\ 1&e^{-\imath 2\pi
f_{2}}&\cdots&e^{-\imath 2\pi f_{2}(L-1)}\\\ \vdots&\vdots&\ddots&\vdots\\\
1&e^{-\imath 2\pi f_{L}}&\cdots&e^{-\imath 2\pi f_{L}(L-1)}\end{pmatrix},$
(10)

where

$\displaystyle f_{k}:=\frac{k-1-L/2}{L}\quad(k=1,\ldots,L)$ (11)

are uniform sampling points in $[-1/2,1/2)$ used as a frequency grid, and $L$
is assumed to be even for simplicity. Then,
$\bar{\mathbf{u}}_{j,n}^{\mathrm{(DFT)}}$ is given by

$\displaystyle\bar{\mathbf{u}}_{j,n}^{\mathrm{(DFT)}}=\mathbf{F}\bar{\mathbf{x}}_{j,n}.$

Due to the unitarity of $\mathbf{F}$, we have

$\displaystyle\bar{\mathbf{x}}_{j,n}=\mathbf{F}^{\mathrm{H}}\bar{\mathbf{u}}_{j,n}^{\mathrm{(DFT)}},$

which corresponds to (7) with $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$.

######  Example 3 (Windowed DFT).

To mitigate the frequency sidelobes, a window function
$\mathbf{w}\in\mathbb{R}_{++}^{L}$ is applied to $\bar{\mathbf{x}}_{j,n}$
before the DFT in some cases [3, 43]. The windowed DFT coefficients are given
by

$\displaystyle\bar{\mathbf{u}}_{j,n}^{\mathrm{(WDFT)}}=\mathbf{F}\mathbf{W}\bar{\mathbf{x}}_{j,n},$

where $\mathbf{F}$ is the DFT matrix in (10), and
$\mathbf{W}:=\mathrm{diag}(\mathbf{w})\in\mathbb{R}_{++}^{L\times L}$. Since
$(\mathbf{F}\mathbf{W})^{-1}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$, we have

$\displaystyle\bar{\mathbf{x}}_{j,n}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}\bar{\mathbf{u}}_{j,n}^{\mathrm{(WDFT)}},$

which corresponds to (7) with
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$.

### Major Challenge in PSD Estimation from Mixtures of Realizations

The square of the magnitude of the frequency components in (8), i.e.,

$\displaystyle|\bar{u}_{j,n}[k]|^{2}\quad(k=1,\ldots,L),$ (12)

is called the _periodogram_ 333Strictly speaking, (12) is called the
periodogram when $\mathbf{G}$ in (7) is the inverse DFT
$\mathbf{F}^{\mathrm{H}}$ in Example 2, and (12) with
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ is called, e.g., the
_windowed periodogram_ [1] or the _modified periodogram_ [3]. and widely used
as an estimate of the PSD $S^{\star}_{n}(f)$ on the frequency grid defined in
(11). It should be noted that the periodogram has to be estimated from the
mixtures of realizations in (9) for our problem. The periodogram with the DFT
shown in Example 2 is an asymptotically unbiased estimator of the PSD under a
mild condition [1, 2]. More precisely, since
$\bar{u}_{j,n}^{\mathrm{(DFT)}}[k]$ can be regarded as a realization of the
random variable

$\displaystyle
U_{n}^{\star}[k]:=\frac{1}{\sqrt{L}}\sum_{\ell=1}^{L}X_{n}^{\star}[\ell]e^{-\imath
2\pi f_{k}(\ell-1)}\quad(k=1,\ldots,L),$

the asymptotic unbiasedness means that

$\displaystyle\lim_{L\rightarrow\infty}E\left[|U_{n}^{\star}[k]|^{2}\right]=S^{\star}_{n}(f_{k})\quad(k=1,\ldots,L).$
(13)

A simple proof of (13) is provided in [16, Section II-A] under a mild
sufficient condition $\sum_{\ell=-\infty}^{\infty}|\ell R_{n}[\ell]|<\infty$.
The periodogram with the windowed DFT in Example 3 is also asymptotically
unbiased if the window function is properly designed [3].

Figure 1: The PSD and the periodogram computed from the realizations for a
simulation.

The drawback of the periodogram lies in its large variance and often-observed
erratic oscillation [1, 2, 3, 4] (see also Fig. 1 for an example from the
experiments in Section IV). Indeed, this drawback is theoretically shown for
the periodogram with the DFT when $X_{n}^{\star}[\ell]$ is a linear
combination of i.i.d. Gaussian random variables and $L\rightarrow\infty$ [1,
2]. Namely, for this case, the variance is as large as the square of the PSD,
i.e.,

$\displaystyle\lim_{L\rightarrow\infty}E\left[\left(|U_{n}^{\star}[k]|^{2}-S^{\star}_{n}(f_{k})\right)^{2}\right]=(S^{\star}_{n}(f_{k}))^{2}$

for $k=1,\ldots,L$, and $U_{n}^{\star}[k]$ and $U_{n}^{\star}[k^{\prime}]$
$(k\neq k^{\prime})$ are uncorrelated when $L\rightarrow\infty$. These facts
validate the often-observed erratic oscillation of the periodogram.

Although several approaches are developed to reduce the variance of the
periodogram, they are not suitable for our problem where the PSDs have to be
estimated from the mixtures of realizations in (9). A simple way is to exploit
the situation that $\bar{u}_{j,n}[k]\quad(j=1,\ldots,J)$ are realizations of
the common random variable, i.e., to use the ensemble average of the
periodograms:

$\displaystyle\frac{1}{J}\sum_{j=1}^{J}|\bar{u}_{j,n}[k]|^{2}\quad(k=1,\ldots,L).$
(14)

Unfortunately, since $J$ is very small for weather radar applications (see
Remark 1), the ensemble average cannot sufficiently reduce the variance and
the erratic oscillation. Since the PSD is usually smooth in weather radar
applications [4], another promising approach is to exploit the smoothness of
the PSD. However, existing smoothing techniques, e.g., those shown in [1, 2],
are not directly applicable to our problem because these techniques suppose
that the frequency components $\bar{\mathbf{u}}_{j,n}$ of the realizations are
known. Using smoothing techniques as a post-processing would be sub-optimal
because the smoothness of the PSD is not considered in the estimation of the
frequency components. Thus, it remains a major challenge to exploit the
smoothness of the PSDs when they have to be estimated from the mixtures of
realizations.

##  III Proposed Approach

To exploit the sparsity and the smoothness for the PSD estimation from the
observed mixtures of realizations in (9), we design a convex model that
jointly estimates the frequency components $\bar{u}_{j,n}[k]$ and the PSDs
$S^{\star}_{n}(f_{k})$. In Section III-A, we first apply the optimally
structured block-sparse model of [38] for the estimation of the frequency
components. Then, in Section III-B, we leverage its latent variable, which is
originally introduced for the block structure optimization, for the estimation
of sparse and smooth PSDs.

###  III-A Block-Sparse Estimation of Frequency Components

We design a block-sparse penalty for the frequency components
$\bar{\mathbf{u}}_{j,n}$ by applying the optimally structured block-sparse
model [38] with the knowledge of the PAWR [16]. For simplicity, we begin by
designing a penalty for each $n\in\\{1,\ldots,N\\}$. As demonstrated in [16],
the PSD $S^{\star}_{n}(f)$ is usually narrow-band for the PAWR, which implies
that $\bar{\mathbf{u}}_{j,n}$ is block-sparse for each source
$n\in\\{1,\ldots,N\\}$ and trial $j\in\\{1,\ldots,J\\}$ due to the relation
(13).444Although we present (13) for $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ in
Example 2 for simplicity, $\bar{\mathbf{u}}_{j,n}$ is also block-sparse when
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ in Example 3 because the
window function in Example 3 is designed to reduce the heights of the
sidelobes and slightly increases the width of the mainlobe. Moreover, since
$\bar{\mathbf{u}}_{j,n}\quad(j=1,\ldots,J)$ are realizations of the common
random variable, suitable block partitions for
$\bar{\mathbf{u}}_{j,n}\quad(j=1,\ldots,J)$ are the same. Thus, using the
mixed $\ell_{2}/\ell_{1}$ norm that is suitable for the block-sparsity, we
introduce a penalty for $\mathbf{u}_{n}:=(\mathbf{u}_{j,n})_{j=1}^{J}$ as

$\displaystyle\left\|\mathbf{u}_{n}\right\|_{2,1}^{(\mathcal{B}_{m,n})_{m=1}^{h_{n}}}:=$
$\displaystyle\sum_{m=1}^{h_{n}}\sqrt{J|\mathcal{B}_{m,n}|}\left\|((u_{j,n}[k])_{j=1}^{J})_{k\in\mathcal{B}_{m,n}}\right\|$
$\displaystyle=$
$\displaystyle\sum_{m=1}^{h_{n}}\sqrt{J|\mathcal{B}_{m,n}|}\sqrt{\sum_{j=1}^{J}\sum_{k\in\mathcal{B}_{m,n}}|u_{j,n}[k]|^{2}},$

where $\mathcal{B}_{m,n}\subset\\{1,\ldots,L\\}\quad(m=1,\ldots,h_{n})$ is a
block partition in the frequency domain of the $n$-th source. By suppressing
the mixed $\ell_{2}/\ell_{1}$ norm, the block-sparsity is promoted because the
components $((u_{j,n}[k])_{j=1}^{J})_{k\in\mathcal{B}_{m,n}}$ in the same
block are forced to be zeros together. The problem in [16] is that an
appropriate block partition is unknown a priori because it depends on the
unknown Doppler velocity distribution. To solve the problem of unknown block
partition, following the approach of [38], we minimize the mixed
$\ell_{2}/\ell_{1}$ norm over the partition of at most $H_{n}$ blocks, i.e.,

$\displaystyle\psi_{H_{n}}(\mathbf{u}_{n}):=\min_{h_{n}\in\\{1,\ldots,H_{n}\\}}\left[\min_{\left(\mathcal{B}_{m,n}\right)_{m=1}^{h_{n}}\in\mathcal{P}_{h_{n}}}\|\mathbf{u}_{n}\|_{2,1}^{(\mathcal{B}_{m,n})_{m=1}^{h_{n}}}\right].$
(15)

The constraint set $\mathcal{P}_{h_{n}}$ consists of all $h_{n}$ block
partitions of $\\{1,\ldots,L\\}$, i.e.,

$\displaystyle\left(\mathcal{B}_{m,n}\right)_{m=1}^{h_{n}}\in\mathcal{P}_{h_{n}}$
$\displaystyle\Leftrightarrow$ $\displaystyle\left\\{\begin{aligned}
&\bigcup_{m=1}^{h_{n}}\mathcal{B}_{m,n}=\\{1,\ldots,L\\},\\\
&\mathcal{B}_{m,n}\neq\varnothing\quad(m=1,\ldots,h_{n}),\\\
&\mathcal{B}_{m,n}\cap\mathcal{B}_{m^{\prime},n}=\varnothing\quad(m\neq
m^{\prime}),\\\
&\mathcal{B}_{m,n}=\mathfrak{M}_{L}\left(\\{\ell\in\mathbb{N}\mid
a_{m,n}\leq\ell\leq b_{m,n}\\}\right)\\\ &\mbox{ for some
}a_{m,n},b_{m,n}\in\mathbb{N}\quad(m=1,\ldots,h_{n}),\end{aligned}\right.$

where

$\displaystyle\mathfrak{M}_{L}(\mathcal{I}):=\left\\{\left.\ell-L\left\lfloor\frac{\ell-1}{L}\right\rfloor\in\\{1,\ldots,L\\}\,\right|\,\ell\in\mathcal{I}\right\\},$

and $\lfloor\cdot\rfloor$ is the floor function. For instance, when
$\mathcal{I}=\\{L-1,L,L+1\\}$, $\mathfrak{M}_{L}(\mathcal{I})=\\{L-1,L,1\\}$.
Differently from the standard design in [38], this design makes
$\mathcal{P}_{h_{n}}$ includes blocks connected by the first and the last
entries, and is suitable for weather radar applications where aliasing may
occur.555Since aliasing is not a serious matter in weather radar applications,
the anti-aliasing filter is usually not employed (see also Footnote 2).

Although $\psi_{H_{n}}(\mathbf{u}_{n})$ is difficult to use directly due to
the combinatorial optimization in (15), we can construct a tight convex
relaxation of $\psi_{H_{n}}(\mathbf{u}_{n})$ as follows. Let
$\phi\colon\mathbb{C}^{J}\times\mathbb{R}_{+}\rightarrow\mathbb{R}_{+}\cup{\\{\infty\\}}$
be a lower semicontinuous convex function defined by

$\displaystyle\phi(\mathbf{v},\sigma):=\begin{dcases}\frac{\|\mathbf{v}\|^{2}}{2\sigma}+\frac{J}{2}\sigma,&\mbox{{\rm
if} }\sigma>0;\\\ 0,&\mbox{{\rm if} }\mathbf{v}=\boldsymbol{0}\mbox{ {\rm and}
}\sigma=0;\\\ \infty,&\mbox{{\rm otherwise}}.\end{dcases}$ (16)

Then, similarly to [38, Section II], $\psi_{H_{n}}(\mathbf{u}_{n})$ can be
rewritten as

$\displaystyle\psi_{H_{n}}\left(\mathbf{u}_{n}\right)=\min_{\begin{subarray}{c}\boldsymbol{\sigma}_{n}\in\mathbb{R}_{+}^{L}\vspace{1.5pt}\\\
\|\mathbf{D}\boldsymbol{\sigma}_{n}\|_{0}\leq
H_{n}-1\end{subarray}}\sum_{k=1}^{L}\phi\left((u_{j,n}[k])_{j=1}^{J},\sigma_{n}[k]\right),$
(17)

where $\mathbf{D}\in\mathbb{R}^{L\times L}$ is the first-order difference
operator with the periodic boundary condition, i.e., the difference operator
on the ring graph [44]. More precisely, $\mathbf{D}$ is defined by

$\displaystyle\mathbf{D}:=\begin{bmatrix}-1&1&0&0&\cdots&0&0\\\
0&-1&1&0&\cdots&0&0\\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&\cdots&-1&1\\\ 1&0&0&0&\cdots&0&-1\end{bmatrix}\in\mathbb{R}^{L\times
L}.$ (18)

Note that (17) is a slight extension of the result in [38] to the case where
the blocks are fixed over the trials $j\in\\{1,\ldots,J\\}$. We can obtain a
tight convex relaxation of (17) by replacing the $\ell_{0}$ pseudo-norm in the
constraint with its best convex relaxation, i.e., the $\ell_{1}$ norm:

$\displaystyle\tilde{\psi}_{\alpha_{n}}(\mathbf{u}_{n}):=\min_{\begin{subarray}{c}\boldsymbol{\sigma}_{n}\in\mathbb{R}_{+}^{L}\vspace{1.5pt}\\\
\|\mathbf{D}\boldsymbol{\sigma}_{n}\|_{1}\leq\alpha_{n}\end{subarray}}\sum_{k=1}^{L}\phi\left((u_{j,n}[k])_{j=1}^{J},\sigma_{n}[k]\right),$

where $\alpha_{n}\in\mathbb{R}_{+}$ is a tuning parameter related to the
number of blocks. Although the sum
$\sum_{n=1}^{N}\tilde{\psi}_{\alpha_{n}}(\mathbf{u}_{n})$ can be used for the
penalty of $(\mathbf{u}_{n})_{n=1}^{N}$, tuning $\alpha_{n}$ for each
$n\in\\{1,\ldots,N\\}$ could be troublesome. Thus, to simplify the tuning, we
propose a convex penalty for $\mathbf{u}:=(\mathbf{u}_{n})_{n=1}^{N}$ as

$\displaystyle\Psi_{\alpha}(\mathbf{u}):=\min_{\begin{subarray}{c}(\boldsymbol{\sigma}_{n})_{n=1}^{N}\in\mathbb{R}_{+}^{NL}\vspace{1.5pt}\\\
\sum_{n=1}^{N}\|\mathbf{D}\boldsymbol{\sigma}_{n}\|_{1}\leq\alpha\end{subarray}}\sum_{n=1}^{N}\sum_{k=1}^{L}\phi\left((u_{j,n}[k])_{j=1}^{J},\sigma_{n}[k]\right),$
(19)

where the single tuning parameter $\alpha\in\mathbb{R}_{+}$ is related to the
number of total blocks for $n=1,\ldots,N$.

Using the proposed convex penalty $\Psi_{\alpha}(\mathbf{u})$ in (19), we
estimate the frequency components by the regularized least squares for the
observation model in (9), i.e.,

$\displaystyle\operatorname*{minimize}_{\mathbf{u}\in\mathbb{C}^{JNL}}\frac{1}{2}\sum_{j=1}^{J}\left\|\mathbf{y}_{j}-\sum_{n=1}^{N}\mathbf{A}_{n}\mathbf{G}\mathbf{u}_{j,n}\right\|^{2}+\lambda\Psi_{\alpha}\left(\mathbf{u}\right),$
(20)

where $\lambda>0$ is the regularization parameter that controls the importance
of the block-sparsity. Substituting the definition of
$\Psi_{\alpha}\left(\mathbf{u}\right)$ in (19) into (20), we can reduce the
optimization problem (20) to

$\left.\begin{aligned}
\operatorname*{minimize}_{\mathbf{u}\in\mathbb{C}^{JNL},\boldsymbol{\sigma}\in\mathbb{R}_{+}^{NL}}\frac{1}{2}&\sum_{j=1}^{J}\left\|\mathbf{y}_{j}-\sum_{n=1}^{N}\mathbf{A}_{n}\mathbf{G}\mathbf{u}_{j,n}\right\|^{2}\\\
+&\lambda\sum_{n=1}^{N}\sum_{k=1}^{L}\phi\left((u_{j,n}[k])_{j=1}^{J},\sigma_{n}[k]\right)&\\\
&\mbox{subject to
}\sum_{n=1}^{N}\|\mathbf{D}\boldsymbol{\sigma}_{n}\|_{1}\leq\alpha\end{aligned}\right\\},$
(21)

where $\boldsymbol{\sigma}$ denotes $(\boldsymbol{\sigma}_{n})_{n=1}^{N}$.
Although the proposed regularization model (21) is a relatively difficult
convex optimization problem due to the discontinuous function $\phi$, we can
obtain a globally optimal solution of (21) by applying the proximal splitting
techniques [45, 46, 47, 48, 49] with the interpretation of $\phi$ in (16) as a
_perspective function_[50, 51, 52]. A concrete algorithm based on the
alternating direction method of multipliers (ADMM) and its derivation are
provided in Appendix A.

###  III-B Leveraging Latent Variable for PSD Estimation

We demonstrate that the solution for the latent variable $\boldsymbol{\sigma}$
of the proposed model (21) is in fact suitable for the PSD estimation. Let
$\hat{\mathbf{u}}$ and $\hat{\boldsymbol{\sigma}}$ be the solutions of (21)
respectively for the variables $\mathbf{u}$ and $\boldsymbol{\sigma}$. While
it is possible to compute the periodogram as in (14) for $\hat{\mathbf{u}}$,
$\hat{\boldsymbol{\sigma}}$ is more suitable for the estimation of smooth
PSDs. To see this, we show that $\hat{\boldsymbol{\sigma}}$ corresponds to the
square root of smoothed and averaged periodogram as follows.

  1. i)

We begin by considering the case of $\alpha\rightarrow\infty$, which is not of
our interest but easy to analyze. In this case, since $\boldsymbol{\sigma}$
minimizes

$\displaystyle\lambda\sum_{n=1}^{N}\sum_{k=1}^{L}\phi\left((u_{j,n}[k])_{j=1}^{J},\sigma_{n}[k]\right)$
(22)

in (21), the solutions $\hat{\mathbf{u}}$ and $\hat{\boldsymbol{\sigma}}$
satisfy the relation

$\displaystyle\hat{\sigma}_{n}[k]=\sqrt{\frac{1}{J}\sum_{j=1}^{J}|\hat{u}_{j,n}[k]|^{2}}$
(23)

for each $k=1,\ldots,L$ and $n=1,\ldots,N$, which can be shown based on [38,
Lemma 1]. The relation (23) means that $\hat{\boldsymbol{\sigma}}$ is the
square root of the averaged periodogram in (14) computed with
$\hat{\mathbf{u}}$ when $\alpha\rightarrow\infty$.

  2. ii)

Next, we consider the case of our interest, where $\alpha$ is set to a finite
value. In this case, the constraint

$\displaystyle\sum_{n=1}^{N}\|\mathbf{D}\boldsymbol{\sigma}_{n}\|_{1}\leq\alpha$

in (21) penalizes the smoothness of
$\boldsymbol{\sigma}_{n}\quad(n=1\,\ldots,N)$ since $\mathbf{D}$ is the
difference operator. Meanwhile, the other part (22) of the proposed model
forces $\boldsymbol{\sigma}$ to be the averaged periodogram in (23). Thus, by
the combination of these terms, roughly speaking, $\hat{\boldsymbol{\sigma}}$
is smoothed around the square root of the averaged periodogram in (23).

Besides the smoothness, $\hat{\boldsymbol{\sigma}}$ is block-sparse because
$\hat{\mathbf{u}}_{j,n}$ $(j=1,\ldots,J)$ are regularized to have a common
block-sparse support. Since the PSDs are smooth and block-sparse for the PAWR,
the square of components of $\hat{\boldsymbol{\sigma}}$, i.e.,

$\displaystyle\hat{S}_{n}(f_{k})=(\hat{\sigma}_{n}[k])^{2}\quad(k=1,\ldots,L),$
(24)

is expected to be a better estimate of the PSDs than the periodogram computed
with $\hat{\mathbf{u}}$.

Intuitively, the proposed model is expected to accurately estimate smooth and
block-sparse PSDs by the following mechanism. Since $\phi(\mathbf{v},\sigma)$
is basically $\frac{\|\mathbf{v}\|^{2}}{2\sigma}+\frac{J}{2}\sigma$ (see
(16)), roughly speaking, the part (22) acts as

$\displaystyle\lambda\sum_{n=1}^{N}\sum_{k=1}^{L}\sum_{j=1}^{J}\left(\frac{|u_{j,n}[k]|^{2}}{2\sigma_{n}[k]}+\frac{\sigma_{n}[k]}{2}\right).$

Since $\sigma_{n}[k]$ estimates the square root of the PSDs,

$\displaystyle\frac{|u_{j,n}[k]|^{2}}{2\sigma_{n}[k]}$

is expected to be an effective regularization for $\mathbf{u}$ because the
expectation of the squared magnitude of the realizations $\bar{\mathbf{u}}$ is
close to the PSDs (see (13) and Fig. 1). Refining $\mathbf{u}$ leads to a
refined estimate of $\boldsymbol{\sigma}$ because $\boldsymbol{\sigma}$ is
smoothed around the value in (23). Thanks to these interactions, the proposed
model is expected to effectively estimate the frequency components and the
PSDs simultaneously.

While $\|\mathbf{D}\boldsymbol{\sigma}_{n}\|_{1}$ with the first difference
operator $\mathbf{D}$ in (18) is a good choice for the control of block
structure, more advanced smoothness priors can be incorporated to further
improve the estimation accuracy of the PSDs. Thanks to the nonnegativity of
$\boldsymbol{\sigma}_{n}$, many convex smoothness penalties designed for real-
valued signals, e.g., the high-order total variation [39, 40] that uses
$\mathbf{D}^{r}$ $(r\geq 2)$ instead of $\mathbf{D}$ and the total generalized
variation [41], can be directly applied to $\boldsymbol{\sigma}_{n}$. For
instance, the proposed model with the high-order total variation

$\left.\begin{aligned}
\operatorname*{minimize}_{\mathbf{u}\in\mathbb{C}^{JNL},\boldsymbol{\sigma}\in\mathbb{R}_{+}^{NL}}\frac{1}{2}&\sum_{j=1}^{J}\left\|\mathbf{y}_{j}-\sum_{n=1}^{N}\mathbf{A}_{n}\mathbf{G}\mathbf{u}_{j,n}\right\|^{2}\\\
+&\lambda\sum_{n=1}^{N}\sum_{k=1}^{L}\phi\left((u_{j,n}[k])_{j=1}^{J},\sigma_{n}[k]\right)&\\\
&\mbox{subject to
}\sum_{n=1}^{N}\|\mathbf{D}^{r}\boldsymbol{\sigma}_{n}\|_{1}\leq\alpha\end{aligned}\right\\}$
(25)

can be solved similarly to the case of (21) by the ADMM-based algorithm shown
in Appendix A. In contrast, when these penalties are applied to, e.g., the
magnitude of $\mathbf{u}$, their convexity is lost (see, e.g., [31]), which
implies that a globally optimal solution is difficult to obtain. Note that
applying these penalties to $\mathbf{u}$, which is complex-valued, is not a
suitable strategy because the magnitude of $\mathbf{u}$ is smooth but the
phase of $\mathbf{u}$ is not smooth in most applications.

TABLE I: Comparison of the methods in terms of the NMAE of the PSDs, where the result is averaged over $100$ independent simulations. Values shown in parenthesis for the existing methods are the NMAEs with the post-smoothing. Settings | MMSE BF |  $\ell_{1}$ reg. |  Mixed $\ell_{2}/\ell_{1}$ reg. | Proposed
---|---|---|---|---
$L=32$, $J=1$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ (standard DFT in Example 2) |  $1.1775$ $(0.9217)$ |  $0.8771$ $(0.6802)$ |  $0.7725$ $(0.6447)$ | $\mathbf{0.6278}$
$L=32$, $J=1$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT in Example 3) |  $1.1055$ $(0.9494)$ |  $0.8236$ $(0.7239)$ |  $0.7644$ $(0.7004)$ | $\mathbf{0.6780}$
$L=128$, $J=1$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ (standard DFT in Example 2) |  $1.1808$ $(0.7396)$ |  $0.8161$ $(0.4874)$ |  $0.6843$ $(0.4906)$ | $\mathbf{0.4693}$
$L=128$, $J=1$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT in Example 3) |  $1.1767$ $(0.7975)$ |  $0.7899$ $(0.5428)$ |  $0.6970$ $(0.5244)$ | $\mathbf{0.5097}$
$L=32$, $J=2$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ (standard DFT in Example 2) |  $0.9345$ $(0.7767)$ |  $0.6917$ $(0.5582)$ |  $0.6201$ $(0.5590)$ | $\mathbf{0.5580}$
$L=32$, $J=2$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT in Example 3) |  $0.9292$ $(0.8318)$ |  $0.6596$ $(0.5865)$ |  $0.6169$ $(0.5790)$ | $\mathbf{0.5718}$
$L=128$, $J=2$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ (standard DFT in Example 2) |  $0.9815$ $(0.7057)$ |  $0.6676$ $(0.4345)$ |  $0.5700$ $(0.4623)$ | $\mathbf{0.4335}$
$L=128$, $J=2$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT in Example 3) |  $0.9763$ $(0.7418)$ |  $0.6435$ $(0.4684)$ |  $0.5723$ $(0.4771)$ | $\mathbf{0.4548}$

##  IV Simulation Results

To demonstrate the effectiveness of the proposed approach, we conduct
numerical simulations on the PSD estimation for the PAWR shown in Example 1.
We basically follow the simulation setting in [12, 16]. The elevation angles
$\theta_{1},\ldots,\theta_{N}$ are chosen uniformly from $None$ to $None$
degrees with $N=110$. We synthesize the (discrete-time) PSD $S_{n}^{\star}(f)$
by

$\displaystyle
S_{n}^{\star}(f)=\frac{1}{T}\sum_{m=-\infty}^{\infty}G_{n}^{\star}\left(\frac{f-m}{T}\right)$

for each $n=1,\ldots,N$, where $T$ is the pulse repetition time, and
$G_{n}^{\star}(f)$ is a continuous-time Gaussian-shaped PSD

$\displaystyle
G_{n}^{\star}(f)=\frac{P_{n}}{\sqrt{2\pi}\varsigma_{n}}e^{-\frac{(f-\mu_{n})^{2}}{2\varsigma_{n}^{2}}},$

which is an appropriate model when, e.g., the atmospheric turbulence is
dominant [4]. The power $P_{n}$ is set from the actual reflection intensity
measured by the PAWR at Osaka University on March 30, 2014. We define the mean
Doppler frequency $\mu_{n}$ by the certain sine curve used in [16]. The
Doppler frequency width $\varsigma_{n}$ is converted from the Doppler velocity
width chosen randomly from the uniform distribution of $[1,3]$
$\mathrm{[m/s]}$. Note that this setting is more realistic than that of [16]
where the Doppler velocity width is simply fixed to $2$ $\mathrm{[m/s]}$ at
every elevation angle. We set $X_{n}^{\star}[\ell]$ in (1) to the Gaussian
process that has the specified PSD $S_{n}^{\star}(f)$, and then generate its
realizations $\bar{x}_{j,n}[\ell]\quad(\ell=1,\ldots,L)$ based on the
probability distribution of $X_{n}^{\star}[\ell]\quad(\ell=1,\ldots,L)$
computed as in [4, 16]. The observation vector $\mathbf{y}_{j}$ is given by
(6), where $\boldsymbol{\varepsilon}_{j}$ is generated as the white Gaussian
noise of the standard deviation $\sqrt{2.5}$. The parameters of the PAWR are
set as follows: $M=128$, $\lambda_{\mathrm{cw}}=31.8\,[\mathrm{mm}]$,
$\Delta=16.5\,[\mathrm{mm}]$, and $T=0.4\,[\mathrm{ms}]$. For the synthesis
matrix $\mathbf{G}$ in the observation model (9) in terms of the frequency
components, we test both $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ and
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ respectively for the
standard DFT in Example 2 and the windowed DFT in Example 3. We use the
hamming window for the window function $\mathbf{w}$ in Example 3, which is
normalized to $\|\mathbf{w}\|=\sqrt{L}$, i.e., to the norm of the rectangular
window $(1,1,\ldots,1)^{\top}$ [43].

We compare the proposed approach that jointly estimates the frequency
components and the PSDs with the existing approach that estimates the PSD
after the estimation of the frequency components. The proposed approach
computes the estimate $\hat{S}_{n}(f_{k})$ of the PSD by (24) with the
solution $\hat{\boldsymbol{\sigma}}$ of the proposed model (25) for the
variable $\boldsymbol{\sigma}$. For frequency component estimators used in the
existing approach, we employ the mixed $\ell_{2}/\ell_{1}$ regularization
model using fixed overlapping blocks of small size [16], which is state-of-
the-art for the PAWR, and the $\ell_{1}$ regularization model as a non-
structured sparse model. For the mixed $\ell_{2}/\ell_{1}$ regularization
model, we adopt the latent group lasso [29, 30] based formulation that selects
relevant blocks from the pre-defined overlapping blocks in the mixed
$\ell_{2}/\ell_{1}$ norm because its estimation accuracy is (slightly) better
than the simple overlapping blocks formulation in [16]. While the above
nonlinear methods outperform linear ones for the frequency component
estimation in [16], we also include the minimum mean square error (MMSE)
beamformer [12], which performs best among the linear methods in [16], for
comparison. Since the MMSE beamformer is a time-domain method that estimates
$\bar{\mathbf{x}}_{j,n}$ from (4), we compute the frequency components from
the estimate of $\bar{\mathbf{x}}_{j,n}$ by using the DFT or the windowed DFT
shown in Examples 2 and 3 respectively. From the estimated frequency
components $\hat{\mathbf{u}}$, the existing approach constructs the estimate
of the PSD as the averaged periodogram

$\displaystyle\hat{S}_{n}^{\mathrm{(AP)}}(f_{k})=\frac{1}{J}\sum_{j=1}^{J}|\hat{u}_{j,n}[k]|^{2}\quad(k=1,\ldots,L).$
(26)

We also test the post-smoothing for the existing approach. Specifically, we
employ the standard smoothing technique, i.e., _Daniell method_ [1, 2]:

$\displaystyle\hat{S}_{n}^{\mathrm{(SAP)}}(f_{k})=\frac{1}{2R+1}\sum_{k^{\prime}=k-R}^{k+R}\hat{S}_{n}^{\mathrm{(AP)}}(f_{k^{\prime}})\quad(k=1,\ldots,L),$

where $2R$ neighbor frequency bins666When $k^{\prime}\notin\\{1,\ldots,L\\}$,
we instead use $k^{\prime}-L\lfloor(k^{\prime}-1)/L\rfloor$ because the anti-
aliasing filter is not employed (see also Footnotes 2 and 5). are used for the
smoothing.

Table I shows the normalized mean absolute error (NMAE)

$\displaystyle\frac{\sum_{n=1}^{N}\sum_{k=1}^{L}\left|S_{n}^{\star}(f_{k})-\hat{S}_{n}(f_{k})\right|}{\sum_{n=1}^{N}\sum_{k=1}^{L}S_{n}^{\star}(f_{k})},$

which is averaged over $100$ independent simulations. The regularization
parameters $\lambda$, $\alpha$ for the proposed model, the block size for the
mixed $\ell_{2}/\ell_{1}$ regularization model, and $R$ in the post-smoothing
are tuned to obtain the best NMAE for each model and setting. We simply set
$r=2$ in the proposed model, though the tuning of $r$ could improve the
estimation accuracy. While the case of $J=1$ is of particular interest in
weather radar applications (see Remark 1), to see the effect of $J$, we also
show the results when $J$ is increased to $2$. In Table I, the proposed model
achieves the best estimation accuracy for all the settings. The post-smoothing
improves the estimation accuracies of the existing models, though their
accuracies are not better than the proposed model.

Figure 2: Ground-truth of the PSDs and their estimates for a simulation of the
PAWR with setting $L=32$, $J=1$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$
(standard DFT). Figure 3: Ground-truth of the PSDs and their estimates for a
simulation of the PAWR with setting $L=32$, $J=1$,
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT). Figure 4:
Ground-truth of the PSDs and their estimates for a simulation of the PAWR with
setting $L=128$, $J=1$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ (standard DFT).
Figure 5: Ground-truth of the PSDs and their estimates for a simulation of the
PAWR with setting $L=128$, $J=1$,
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT). Figure 6:
Ground-truth of the PSDs and their estimates for a simulation of the PAWR with
setting $L=32$, $J=2$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ (standard DFT).
Figure 7: Ground-truth of the PSDs and their estimates for a simulation of the
PAWR with setting $L=32$, $J=2$,
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT). Figure 8:
Ground-truth of the PSDs and their estimates for a simulation of the PAWR with
setting $L=128$, $J=2$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$ (standard DFT).
Figure 9: Ground-truth of the PSDs and their estimates for a simulation of the
PAWR with setting $L=128$, $J=2$,
$\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$ (windowed DFT).

We show the ground-truth and the estimates of the PSDs for examples of
simulations: Fig. 2 for $L=32$, $J=1$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$,
Fig. 3 for $L=32$, $J=1$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$,
Fig. 4 for $L=128$, $J=1$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$, Fig. 5 for
$L=128$, $J=1$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$, Fig. 6
for $L=32$, $J=2$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$, Fig. 7 for $L=32$,
$J=2$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$, Fig. 8 for
$L=128$, $J=2$, $\mathbf{G}=\mathbf{F}^{\mathrm{H}}$, Fig. 9 for $L=128$,
$J=2$, $\mathbf{G}=\mathbf{W}^{-1}\mathbf{F}^{\mathrm{H}}$. From Figs.
2–9(b)(d)(f), we see that the estimates of the existing models exhibit erratic
oscillation as they do not exploit the smoothness of the PSDs. In Figs.
6–9(b)(d)(f) where $J$ is increased to $2$, the erratic oscillation is
slightly reduced but still clearly visible, which suggests the limitation of
the ensemble average (26) when $J$ is small (see Remark 1 for the reason why
$J$ is set small in weather radar applications). The post-smoothing reduces
the erratic oscillation to some extent as shown in Figs. 2–9(c)(e)(g).
However, the sparsity of the estimate is impaired, i.e., the number of entries
of large magnitude that are not present in the ground-truth increases, because
the sparsity is not considered in the post-smoothing. Note that the impairment
of the sparsity is more noticeable in Figs. 4, 5, 8, and 9 for $L=128$. In
contrast, in Figs. 2–9(h), the proposed approach obtains the estimates that
have both the sparsity and the smoothness.

We also see that the erratic oscillation is still visible for the MMSE
beamformer and the $\ell_{1}$ regularization model after the post-smoothing in
Figs. 2–9(c)(e), which suggests the limitation of the post-smoothing. In
particular, while the objective accuracies of the $\ell_{1}$ regularization
model is close to the proposed model when $J$ is increased to $2$ and the
standard DFT is used, the erratic oscillation is not eliminated as seen from
Figs. 6(e) and 8(e). Compared to the $\ell_{1}$ regularization model, the
erratic oscillation is considerably reduced in the mixed $\ell_{2}/\ell_{1}$
regularization model after the post-smoothing in Figs. 2–9(g). This is
hypothesized because the block-sparse model does not necessarily promote the
smoothness but the components in the same block are forced to be zeros or
nonzeros together, which would be beneficial to promote the smoothness by the
post-processing. Even compared with the second-best mixed $\ell_{2}/\ell_{1}$
regularization model, the proposed model yields more smooth estimates in
particular when $L=128$ in Figs. 4, 5, 8, and 9.

The line-like artifacts in the estimates shown in Figs. 2, 4, 6, and 8 are
more or less reduced in the estimates shown in Figs. 3, 5, 7, and 9 thanks to
the window function that reduces the heights of the sidelobes. Although the
estimates with the windowed DFT are visually closer to the ground-truth than
those with the standard DFT, the objective accuracy shown in Table I is not
always improved perhaps because the window function increases the width of the
mainlobe. Since the line-like artifacts are caused due to the sidelobes of the
window function with finite $L$, the line-like artifacts are more reduced when
$L=128$. In particular, the line-like artifacts are almost completely
eliminated in the proposed estimates when $L=128$ in Figs. 5(h) and 9(h).

##  V Conclusion

We presented a convex optimization model for the estimation of sparse and
smooth PSDs of complex-valued random processes from noisy mixtures of
realizations. While the PSDs are related to the expectation of the magnitude
of the frequency components of the realizations, it has been difficult to
exploit the smoothness of the PSDs as naive penalties for the difference of
the magnitude of the frequency components induce hard nonconvex optimization
problems. To resolve this difficulty, we designed the proposed model that
jointly estimates complex-valued frequency components and nonnegative PSDs.
More precisely, we first applied the optimally structured block-sparse model
of [38] for the frequency component estimation, and then newly leveraged the
latent variable of the model, which was originally introduced to optimize the
block structure, for the PSD estimation. We demonstrate that the latent
variable is in fact related to the square root of the PSDs, which enables to
exploit the smoothness of the PSDs via convex optimization by penalizing the
difference of the nonnegative latent variable. Moreover, many smoothness
priors designed for real-valued signals can be readily incorporated into the
proposed framework for enhancement of the smoothness of the PSDs of complex-
valued random processes. Numerical experiments on the phased array weather
radar showed that the proposed approach achieved better objective accuracy and
yielded visually better estimates compared with the existing sparse estimation
models even when combined with post-smoothing.

Input: $\gamma>0$,
$\mathbf{A}:=[\mathbf{A}_{1},\ldots,\mathbf{A}_{N}]\in\mathbb{C}^{d\times
NL}$, $\mathbf{x}_{j,n}^{(0)}\in\mathbb{C}^{L}$,
$\tilde{\mathbf{u}}_{j,n}^{(0)}\in\mathbb{C}^{L}$,
$\tilde{\boldsymbol{\sigma}}_{n}^{(0)}\in\mathbb{R}_{+}^{L}$,
$\boldsymbol{\eta}_{n}^{(0)}\in\mathbb{R}^{L}$,
$\mathbf{q}_{j,n}^{(0)}\in\mathbb{C}^{L}$,
$\mathbf{r}_{j,n}^{(0)}\in\mathbb{C}^{L}$,
$\boldsymbol{\tau}_{n}^{(0)}\in\mathbb{R}^{L}$,
$\boldsymbol{\zeta}_{n}^{(0)}\in\mathbb{R}^{L}$ $(j\in\\{1,2,\ldots,J\\}\mbox{
and }n\in\\{1,2,\ldots,N\\})$.

for  _$i=0,1,2,\ldots$_ do

for  _$n\in\\{1,2,\ldots,N\\}$_ do

for  _$j\in\\{1,2,\ldots,J\\}$_ do

$\mathbf{u}_{j,n}^{(i+1)}=(\mathbf{I}_{L}+\mathbf{G}^{\mathrm{H}}\mathbf{G})^{-1}\left(\mathbf{G}^{\mathrm{H}}(\mathbf{x}_{j,n}^{(i)}-\mathbf{q}_{j,n}^{(i)})+\tilde{\mathbf{u}}_{j,n}^{(i)}-\mathbf{r}_{j,n}^{(i)}\right)$

$\boldsymbol{\sigma}_{n}^{(i+1)}=(\mathbf{I}_{L}+(\mathbf{D}^{r})^{\top}\mathbf{D}^{r})^{-1}\left((\mathbf{D}^{r})^{\top}(\boldsymbol{\eta}_{n}^{(i)}-\boldsymbol{\zeta}_{n}^{(i)})+\tilde{\boldsymbol{\sigma}}_{n}^{(i)}-\boldsymbol{\tau}_{n}^{(i)}\right)$

for  _$j\in\\{1,2,\ldots,J\\}$_ do

$\left(\mathbf{x}_{j,n}^{(i+1)}\right)_{n=1}^{N}=(\mathbf{I}_{NL}+\gamma\mathbf{A}^{\mathrm{H}}\mathbf{A})^{-1}\left(\gamma\mathbf{A}^{\mathrm{H}}\mathbf{y}_{j}+\left(\mathbf{G}\mathbf{u}_{j,n}^{(i+1)}+\mathbf{q}_{j,n}^{(i)}\right)_{n=1}^{N}\right)$

for  _$(n,k)\in\\{1,2,\ldots,N\\}\times\\{1,2,\ldots,L\\}$_ do

$\left((\tilde{u}_{j,n}^{(i+1)}[k])_{j=1}^{J},\tilde{\sigma}_{n}^{(i+1)}[k]\right)=\mathrm{prox}_{\gamma\lambda\phi}\left(\left(u_{j,n}^{(i+1)}[k]+r_{j,n}^{(i)}[k]\right)_{j=1}^{J},\sigma_{n}^{(i+1)}[k]+\tau_{n}^{(i)}[k]\right)$

$\left(\boldsymbol{\eta}_{n}^{(i+1)}\right)_{n=1}^{N}=P_{B_{1}^{\alpha}}\left(\left(\mathbf{D}^{r}\boldsymbol{\sigma}_{n}^{(i+1)}+\boldsymbol{\zeta}_{n}^{(i)}\right)_{n=1}^{N}\right)$

for  _$n\in\\{1,2,\ldots,N\\}$_ do

for  _$j\in\\{1,2,\ldots,J\\}$_ do

$\mathbf{q}_{j,n}^{(i+1)}=\mathbf{q}_{j,n}^{(i)}+\mathbf{G}\mathbf{u}_{j,n}^{(i+1)}-\mathbf{x}_{j,n}^{(i+1)}$

$\mathbf{r}_{j,n}^{(i+1)}=\mathbf{r}_{j,n}^{(i)}+\mathbf{u}_{j,n}^{(i+1)}-\tilde{\mathbf{u}}_{j,n}^{(i+1)}$

$\boldsymbol{\tau}_{n}^{(i+1)}=\boldsymbol{\tau}_{n}^{(i)}+\boldsymbol{\sigma}_{n}^{(i+1)}-\tilde{\boldsymbol{\sigma}}_{n}^{(i+1)}$

$\boldsymbol{\zeta}_{n}^{(i+1)}=\boldsymbol{\zeta}_{n}^{(i)}+\mathbf{D}^{r}\boldsymbol{\sigma}_{n}^{(i+1)}-\boldsymbol{\eta}_{n}^{(i+1)}$

Algorithm 1 Iterative solver for the proposed regularization model (25)

##  Appendix A Solver for Proposed Regularization Model

The proposed regularization model (25) can be solved by the proximal splitting
techniques [45, 46, 47, 48, 49] with the closed-form computation of the
proximity operator of $\phi$. As a concrete example, using the alternating
direction method of multipliers (ADMM) [46, 47], we provide an iterative
solver that is guaranteed to converge to an optimal solution of (25). The ADMM
solves the following convex optimization problem

$\displaystyle\operatorname*{minimize}_{v\in\mathcal{V},w\in\mathcal{W}}F(v)+G(w)\mbox{
subject to }w=\mathcal{L}v$ (27)

by the iterations

$\left\lfloor\begin{aligned}
v^{(i+1)}&\in\arg\min_{v\in\mathcal{V}}\left[\gamma
F(v)+\frac{1}{2}\|w^{(i)}-z^{(i)}-\mathcal{L}v\|^{2}\right]\\\
w^{(i+1)}&=\mathrm{prox}_{\gamma G}(\mathcal{L}v^{(i+1)}+z^{(i)})\\\
z^{(i+1)}&=z^{(i)}+\mathcal{L}v^{(i+1)}-w^{(i+1)},\end{aligned}\right.$

where we suppose that $\mathcal{V}$ and $\mathcal{W}$ are finite-dimensional
Hilbert spaces, $F$ and $G$ are proper lower semicontinuous convex functions,
$\mathcal{L}$ is a linear operator, $\mathrm{prox}_{\gamma
G}(w):=\mathrm{arg}\min_{\omega\in\mathcal{W}}\left[\gamma
G(\omega)+\frac{1}{2}\|w-\omega\|^{2}\right]$ is the proximity operator of
$\gamma G$, and $\gamma>0$.

To apply the ADMM, we rewrite (25) as

$\left.\begin{aligned}
\operatorname*{minimize}_{\mathbf{u},\boldsymbol{\sigma},\mathbf{x},\tilde{\mathbf{u}},\tilde{\boldsymbol{\sigma}},\boldsymbol{\eta}}&F(\mathbf{u},\boldsymbol{\sigma})+G(\mathbf{x},\tilde{\mathbf{u}},\tilde{\boldsymbol{\sigma}},\boldsymbol{\eta})\\\
\mbox{subject to }&\mathbf{x}_{j,n}=\mathbf{G}\mathbf{u}_{j,n}\quad(\forall
j,n)\\\ &\tilde{\mathbf{u}}_{j,n}=\mathbf{u}_{j,n}\quad(\forall j,n)\\\
&\tilde{\boldsymbol{\sigma}}_{n}=\boldsymbol{\sigma}_{n}\quad(\forall n)\\\
&\boldsymbol{\eta}_{n}=\mathbf{D}^{r}\boldsymbol{\sigma}_{n}\quad(\forall
n)\end{aligned}\right\\},$ (28)

where we let

$\displaystyle F(\mathbf{u},\boldsymbol{\sigma})$ $\displaystyle:=0,$
$\displaystyle
G(\mathbf{x},\tilde{\mathbf{u}},\tilde{\boldsymbol{\sigma}},\boldsymbol{\eta})$
$\displaystyle:=\frac{1}{2}\sum_{j=1}^{J}\left\|\mathbf{y}_{j}-\sum_{n=1}^{N}\mathbf{A}_{n}\mathbf{x}_{j,n}\right\|^{2}$
$\displaystyle+\lambda$
$\displaystyle\sum_{n=1}^{N}\sum_{k=1}^{L}\phi\left((\tilde{u}_{j,n}[k])_{j=1}^{J},\tilde{\sigma}_{n}[k]\right)+\iota_{B_{1}^{\alpha}}\left(\boldsymbol{\eta}\right),$

and $\iota_{B_{1}^{\alpha}}(\boldsymbol{\eta})$ is the indicator function of
the $\ell_{1}$ ball, i.e.,

$\displaystyle\iota_{B_{1}^{\alpha}}(\boldsymbol{\eta}):=\begin{dcases}0,&\mbox{if
}\sum_{n=1}^{N}\|\boldsymbol{\eta}_{n}\|_{1}\leq\alpha;\\\
\infty,&\mbox{otherwise}.\end{dcases}$

Since the constraint of (28) can be expressed in the form of (27), we can
apply the ADMM to (28), and obtain the iterative algorithm shown in Algorithm
1. For our formulation, since the minimizer of the first step of the ADMM is
unique, the convergence to an optimal solution of (28) follows from [49]. From
the equivalence between (25) and (28), $(\mathbf{u}^{(i)})_{i=1}^{\infty}$ and
$(\boldsymbol{\sigma}^{(i)})_{i=1}^{\infty}$ generated by Algorithm 1
converges to the solution of (25) for the variables $\mathbf{u}$ and
$\boldsymbol{\sigma}$ respectively.

The operators in Algorithm 1 can be computed as follows. Expressing
$\phi(\mathbf{v},\sigma)$ as the sum of the perspective function [50, 51, 52]
of $\frac{\|\mathbf{v}\|^{2}}{2}$ and the linear function $\frac{J}{2}\sigma$,
based on [52, Example 2.4], we can compute $\mathrm{prox}_{\kappa\phi}$ for
$\kappa=\gamma\lambda$ by

$\displaystyle\mathrm{prox}_{\kappa\phi}(\mathbf{v},\sigma)$
$\displaystyle=\begin{cases}(\boldsymbol{0},0),&\mbox{if
}2\kappa\sigma+\|\mathbf{v}\|^{2}\leq J\kappa^{2};\\\
(\boldsymbol{0},\sigma-\frac{\kappa J}{2}),&\mbox{if
}\mathbf{v}=\boldsymbol{0}\mbox{ and }2\sigma>J\kappa;\\\
\left(\mathbf{v}-\kappa
s\frac{\mathbf{v}}{\|\mathbf{v}\|},\sigma+\kappa\frac{s^{2}-J}{2}\right),&\mbox{otherwise},\end{cases}$
(29)

where $s>0$ is the unique positive root of

$\displaystyle
s^{3}+\left(\frac{2}{\kappa}\sigma+2-J\right)s-\frac{2}{\kappa}\|\mathbf{v}\|=0,$

and can be explicitly given based on Cardano’s formula as follows [38, 53].
Let $p=\frac{2}{\kappa}\sigma+2-J$ and
$D=-\frac{\|\mathbf{v}\|^{2}}{\kappa^{2}}-\frac{p^{3}}{27}$. Then,

$\displaystyle
s=\begin{dcases}\sqrt[3]{\frac{\|\mathbf{v}\|}{\kappa}+\sqrt{-D}}+\sqrt[3]{\frac{\|\mathbf{v}\|}{\kappa}-\sqrt{-D}},&\mbox{if
}D<0;\\\ 2\sqrt[3]{\frac{\|\mathbf{v}\|}{\kappa}},&\mbox{if }D=0;\\\
2\sqrt{-\frac{p}{3}}\cos\left(\frac{\arctan(\kappa\sqrt{D}/\|\mathbf{v}\|)}{3}\right),&\mbox{if
}D>0,\end{dcases}$

where $\sqrt[3]{\cdot}$ denotes the real cubic root. The $\ell_{1}$ ball
projection $P_{B_{1}^{\alpha}}$, which is the proximity operator of
$\iota_{B_{1}^{\alpha}}(\boldsymbol{\eta})$, can be computed in
$\mathcal{O}(NL)$ expected complexity, e.g., by the algorithm of [54].

### Implementation for weather radar applications

The matrix inversions in Algorithm 1 can be efficiently computed for the
application to weather radars as follows.

  1. a)

To efficiently compute the inversion of
$(\mathbf{I}_{L}+\mathbf{G}^{\mathrm{H}}\mathbf{G})$, we use the property

$\displaystyle\mathbf{G}\mathbf{G}^{\mathrm{H}}=\mathrm{diag}(\boldsymbol{\nu})\in\mathbb{R}_{++}^{L},$

which holds for, e.g., the DFT in Example 2 and the windowed DFT in Example 3.
More precisely, from this property, we have

$\displaystyle(\mathbf{I}_{L}+\mathbf{G}^{\mathrm{H}}\mathbf{G})^{-1}$
$\displaystyle=\mathbf{I}_{L}-\mathbf{G}^{\mathrm{H}}(\mathbf{I}_{L}+\mathbf{G}\mathbf{G}^{\mathrm{H}})^{-1}\mathbf{G}$
$\displaystyle=\mathbf{I}_{L}-\mathbf{G}^{\mathrm{H}}(\mathbf{I}_{L}+\mathrm{diag}(\boldsymbol{\nu}))^{-1}\mathbf{G}$
$\displaystyle=\mathbf{I}_{L}-\mathbf{G}^{\mathrm{H}}\mathrm{diag}\left(\left(\frac{1}{1+\nu_{\ell}}\right)_{\ell=1}^{L}\right)\mathbf{G},$

where the first equality follows from the Sherman-Morrison-Woodbury matrix
inversion lemma [55]. Note that the multiplications of $\mathbf{G}$ and
$\mathbf{G}^{\mathrm{H}}$ can be computed in $\mathcal{O}(L\log L)$ by the
fast Fourier transform (FFT) for Examples 2 and 3.

  2. b)

Since $\mathbf{D}$ in (18) is a circulant matrix, the inversion of
$(\mathbf{I}_{L}+(\mathbf{D}^{r})^{\top}\mathbf{D}^{r})$ can be computed in
$\mathcal{O}(L\log L)$ with the FFT [55].

  3. c)

For the PAWR shown in Example 1, the inversion of
$\mathbf{I}_{NL}+\mathbf{A}^{\mathrm{H}}\mathbf{A}\in\mathbb{C}^{NL\times NL}$
can be computed in $\mathcal{O}(N^{3})$, independently of the value of $L$,
because it can be translated into a block-diagonal matrix after some
permutations. We show this explicitly in another way. Notice that
$(\mathbf{x}_{j,n}^{(i+1)})_{n=1}^{N}$ in Algorithm 1 is the solution of

$\left.\begin{aligned}
\operatorname*{minimize}_{\left(\mathbf{x}_{j,n}\right)_{n=1}^{N}}\frac{\gamma}{2}&\left\|\mathbf{y}_{j}-\sum_{n=1}^{N}\mathbf{A}_{n}\mathbf{x}_{j,n}\right\|^{2}\\\
+\frac{1}{2}&\sum_{n=1}^{N}\|\mathbf{G}\mathbf{u}_{j,n}^{(i+1)}+\mathbf{q}_{j,n}^{(i)}-\mathbf{x}_{j,n}\|^{2}\end{aligned}\right\\}.$
(30)

From the definitions in Example 1, we have

$\displaystyle\left\|\mathbf{y}_{j}-\sum_{n=1}^{N}\mathbf{A}_{n}\mathbf{x}_{j,n}\right\|^{2}$
$\displaystyle=$
$\displaystyle\left\|\begin{pmatrix}\mathbf{y}_{j}[1]-\sum_{n=1}^{N}x_{j,n}[1]\mathbf{a}(\theta_{n})\\\
\vdots\\\
\mathbf{y}_{j}[L]-\sum_{n=1}^{N}x_{j,n}[L]\mathbf{a}(\theta_{n})\end{pmatrix}\right\|^{2}$
$\displaystyle=$
$\displaystyle\left\|\begin{pmatrix}\mathbf{y}_{j}[1]^{\top}-\sum_{n=1}^{N}x_{j,n}[1]\mathbf{a}(\theta_{n})^{\top}\\\
\vdots\\\
\mathbf{y}_{j}[L]^{\top}-\sum_{n=1}^{N}x_{j,n}[L]\mathbf{a}(\theta_{n})^{\top}\end{pmatrix}\right\|_{\mathrm{fro}}^{2}$
$\displaystyle=$
$\displaystyle\|\mathbf{Y}_{j}^{\top}-\mathbf{X}_{j}\mathbf{S}^{\top}\|_{\mathrm{fro}}^{2},$

where we let

$\displaystyle\mathbf{Y}_{j}:=$
$\displaystyle(\mathbf{y}_{j}[1],\ldots,\mathbf{y}_{j}[L])\in\mathbb{C}^{M\times
L},$ $\displaystyle\mathbf{X}_{j}:=$
$\displaystyle(\mathbf{x}_{j,1},\ldots,\mathbf{x}_{j,N})\in\mathbb{C}^{L\times
N},$ $\displaystyle\mathbf{S}:=$
$\displaystyle(\mathbf{a}(\theta_{1}),\ldots,\mathbf{a}(\theta_{N}))\in\mathbb{C}^{M\times
N},$

and $\|\cdot\|_{\mathrm{fro}}$ denotes the Frobenius norm of the matrix. From
this expression, it is clear that the step (30) can be solved in
$\mathcal{O}(N^{3})$ for the inversion regarding $\mathbf{S}^{\top}$. Namely,
$(\mathbf{x}_{j,n}^{(i+1)})_{n=1}^{N}$ in Algorithm 1 is obtained by

$\displaystyle\Bigl{(}$
$\displaystyle\mathbf{x}_{j,1}^{(i+1)},\ldots,\mathbf{x}_{j,N}^{(i+1)}\Bigr{)}$
$\displaystyle=\left(\gamma\mathbf{Y}_{j}^{\top}\mathbf{S}^{*}+\mathbf{G}\mathbf{U}_{j}^{(i+1)}+\mathbf{Q}_{j}^{(i)}\right)(\mathbf{I}_{N}+\gamma\mathbf{S}^{\top}\mathbf{S}^{*})^{-1},$

where $\mathbf{S}^{*}$ is the complex conjugate of $\mathbf{S}$, and

$\displaystyle\mathbf{U}_{j}^{(i+1)}:=$
$\displaystyle\left(\mathbf{u}_{j,1}^{(i+1)},\ldots,\mathbf{u}_{j,N}^{(i+1)}\right)\in\mathbb{C}^{L\times
N},$ $\displaystyle\mathbf{Q}_{j}^{(i)}:=$
$\displaystyle\left(\mathbf{q}_{j,1}^{(i)},\ldots,\mathbf{q}_{j,N}^{(i)}\right)\in\mathbb{C}^{L\times
N}.$

We remark that the inversions in Algorithm 1 are the same for all iterations,
and thus they can be computed in advance and stored in the memory.

## References

  * [1] P. Stoica and R. L. Moses, _Spectral Analysis of Signals_. Prentice-Hall, 2005.
  * [2] P. J. Brockwell and R. A. Davis, _Time Series: Theory and Methods_. Springer, 2009.
  * [3] A. V. Oppenheim and R. W. Schafer, _Discrete-Time Signal Processing_ , 3rd ed. Prentice-Hall, 2009.
  * [4] V. N. Bringi and V. Chandrasekar, _Polarimetric Doppler Weather Radar: Principles and Applications_. Cambridge University Press, 2001.
  * [5] D. S. Zrnic and R. J. Doviak, “Velocity spectra of vortices scanned with a pulse-Doppler radar,” _J. Appl. Meteorol._ , vol. 14, no. 8, pp. 1531–1539, 1975.
  * [6] L. H. Janssen and G. A. Van Der Spek, “The shape of Doppler spectra from precipitation,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. AES-21, no. 2, pp. 208–219, 1985.
  * [7] Y. Wang, T.-Y. Yu, M. Yeary, A. Shapiro, S. Nemati, M. Foster, D. L. Andra, and M. Jain, “Tornado detection using a neuro–fuzzy system to integrate shear and spectral signatures,” _J. Atmos. Ocean. Technol._ , vol. 25, no. 7, pp. 1136–1148, 2008.
  * [8] T.-Y. Yu, R. R. Rondinel, and R. D. Palmer, “Investigation of non-Gaussian Doppler spectra observed by weather radar in a tornadic supercell,” _J. Atmos. Ocean. Technol._ , vol. 26, no. 3, pp. 444–461, 2009.
  * [9] D. S. Zrnic, J. F. Kimpel, D. E. Forsyth, A. Shapiro, G. Crain, R. Ferek, J. Heimmer, W. Benner, F. T. J. McNellis, and R. J. Vogt, “Agile-beam phased array radar for weather observations,” _Bull. Am. Meteorol. Soc._ , vol. 88, no. 11, pp. 1753–1766, 2007.
  * [10] D. McLaughlin, D. Pepyne, V. Chandrasekar, B. Philips, J. Kurose, M. Zink, K. Droegemeier, S. Cruz-Pol, F. Junyent, J. Brotzge, D. Westbrook, N. Bharadwaj, Y. Wang, E. Lyons, K. Hondl, Y. Liu, E. Knapp, M. Xue, A. Hopf, K. Kloesel, A. DeFonzo, P. Kollias, K. Brewster, R. Contreras, B. Dolan, T. Djaferis, E. Insanic, S. Frasier, and F. Carr, “Short-wavelength technology and the potential for distributed networks of small radar systems,” _Bull. Am. Meteorol. Soc._ , vol. 90, no. 12, pp. 1797–1818, 2009\.
  * [11] B. Isom, R. Palmer, R. Kelley, J. Meier, D. Bodine, M. Yeary, B.-L. Cheong, Y. Zhang, T.-Y. Yu, and M. I. Biggerstaff, “The atmospheric imaging radar: Simultaneous volumetric observations using a phased array weather radar,” _J. Atmos. Ocean. Technol._ , vol. 30, no. 4, pp. 655–675, Apr. 2013.
  * [12] E. Yoshikawa, T. Ushio, Z. Kawasaki, S. Yoshida, T. Morimoto, F. Mizutani, and M. Wada, “MMSE beam forming on fast-scanning phased array weather radar,” _IEEE Trans. Geosci. Remote Sens._ , vol. 51, no. 5, pp. 3077–3088, May 2013.
  * [13] D. S. Zrnic, V. M. Melnikov, R. J. Doviak, and R. Palmer, “Scanning strategy for the multifunction phased-array radar to satisfy aviation and meteorological needs,” _IEEE Geosci. Remote Sens. Lett._ , vol. 12, no. 6, pp. 1204–1208, 2015.
  * [14] F. Mizutani, T. Ushio, E. Yoshikawa, S. Shimamura, H. Kikuchi, M. Wada, S. Satoh, and T. Iguchi, “Fast-scanning phased-array weather radar with angular imaging technique,” _IEEE Trans. Geosci. Remote Sens._ , vol. 56, no. 5, pp. 2664–2673, May 2018.
  * [15] E. Yoshikawa, T. Ushio, and H. Kikuchi, “A study of comb beam transmission on phased array weather radars,” _IEEE Trans. Geosci. Remote Sens._ , vol. 59, no. 8, pp. 6346–6356, 2021.
  * [16] D. Kitahara, H. Kuroda, A. Hirabayashi, E. Yoshikawa, H. Kikuchi, and T. Ushio, “Nonlinear beamforming based on group-sparsities of periodograms for phased array weather radar,” _IEEE Trans. Geosci. Remote Sens._ , 2022.
  * [17] 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–3022, 2005.
  * [18] 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.
  * [19] L. Wang, L. Zhao, G. Bi, C. Wan, L. Zhang, and H. Zhang, “Novel wideband DOA estimation based on sparse Bayesian learning with Dirichlet process priors,” _IEEE Trans. Signal Process._ , vol. 64, no. 2, pp. 275–289, 2016\.
  * [20] 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.
  * [21] C. R. Rojas, D. Katselis, and H. Hjalmarsson, “A note on the SPICE method,” _IEEE Trans. Signal Process._ , vol. 61, no. 18, pp. 4545–4551, 2013.
  * [22] 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.
  * [23] X. Shang, J. Li, and P. Stoica, “Weighted SPICE algorithms for range-Doppler imaging using one-bit automotive radar,” _IEEE J. Sel. Topics Signal Process._ , vol. 15, no. 4, pp. 1041–1054, 2021.
  * [24] M. Yuan and Y. Lin, “Model selection and estimation in regression with grouped variables,” _J. Roy. Statist. Soc. Ser. B_ , vol. 68, no. 1, pp. 49–67, 2006\.
  * [25] M. Stojnic, F. Parvaresh, and B. Hassibi, “On the reconstruction of block-sparse signals with an optimal number of measurements,” _IEEE Trans. Signal Process._ , vol. 57, no. 8, pp. 3075–3085, Aug. 2009.
  * [26] Y. C. Eldar, P. Kuppinger, and H. Bolcskei, “Block-sparse signals: Uncertainty relations and efficient recovery,” _IEEE Trans. Signal Process._ , vol. 58, no. 6, pp. 3042–3054, Jun. 2010.
  * [27] X. Lv, G. Bi, and C. Wan, “The group lasso for stable recovery of block-sparse signal representations,” _IEEE Trans. Signal Process._ , vol. 59, no. 4, pp. 1371–1382, Apr. 2011.
  * [28] E. Elhamifar and R. Vidal, “Block-sparse recovery via convex optimization,” _IEEE Trans. Signal Process._ , vol. 60, no. 8, pp. 4094–4107, Aug. 2012\.
  * [29] L. Jacob, G. Obozinski, and J.-P. Vert, “Group lasso with overlap and graph lasso,” in _Proc. Int. Conf. Mach. Learn. (ICML)_ , Jun. 2009, pp. 433–440.
  * [30] G. Obozinski, L. Jacob, and J.-P. Vert, “Group lasso with overlaps: The latent group lasso approach,” preprint arXiv:1110.0413, 2011, 60 pages.
  * [31] P. Ciuciu, J. Idier, and J.-F. Giovannelli, “Regularized estimation of mixed spectra using a circular Gibbs-Markov model,” _IEEE Trans. Signal Process._ , vol. 49, no. 10, pp. 2202–2213, 2001.
  * [32] D. D. Ariananda and G. Leus, “Compressive wideband power spectrum estimation,” _IEEE Trans. Signal Process._ , vol. 60, no. 9, pp. 4775–4789, 2012.
  * [33] L. Wang, T. Gerkmann, and S. Doclo, “Noise power spectral density estimation using maxNSR blocking matrix,” _IEEE/ACM Trans. Audio, Speech, and Language Process._ , vol. 23, no. 9, pp. 1493–1508, 2015.
  * [34] Y. Lu and P. C. Loizou, “A geometric approach to spectral subtraction,” _Speech Communication_ , vol. 50, no. 6, pp. 453–466, 2008.
  * [35] Y. Zhao, X. Zhao, and B. Wang, “A speech enhancement method based on sparse reconstruction of power spectral density,” _Computers & Electrical Engineering_, vol. 40, no. 4, pp. 1080–1089, 2014.
  * [36] J. A. Bazerque and G. B. Giannakis, “Distributed spectrum sensing for cognitive radio networks by exploiting sparsity,” _IEEE Trans. Signal Process._ , vol. 58, no. 3, pp. 1847–1862, 2010.
  * [37] F. Li and Z. Xu, “Sparse Bayesian hierarchical prior modeling based cooperative spectrum sensing in wideband cognitive radio networks,” _IEEE Signal Process. Lett._ , vol. 21, no. 5, pp. 586–590, 2014.
  * [38] H. Kuroda and D. Kitahara, “Block-sparse recovery with optimal block partition,” _IEEE Trans. Signal Process._ , vol. 70, pp. 1506–1520, 2022\.
  * [39] T. Chan, A. Marquina, and P. Mulet, “High-order total variation-based image restoration,” _SIAM J. Sci. Comput._ , vol. 22, no. 2, pp. 503–516, 2000\.
  * [40] G. D. Maso, I. Fonseca, G. Leoni, and M. Morini, “A higher order model for image restoration: The one-dimensional case,” _SIAM J. Math. Anal._ , vol. 40, no. 6, pp. 2351–2391, 2009.
  * [41] K. Bredies, K. Kunisch, and T. Pock, “Total generalized variation,” _SIAM J. Imaging Sci._ , vol. 3, no. 3, pp. 492–526, 2010.
  * [42] H. Kuroda, D. Kitahara, E. Yoshikawa, H. Kikuchi, and T. Ushio, “Sparsity-smoothness-aware power spectral density estimation with application to phased array weather radar,” in _Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP)_ , 2023, pp. 1–5.
  * [43] V. Chandrasekar, R. M. Beauchamp, and R. Bechini, _Introduction to Dual Polarization Weather Radar: Fundamentals, Applications, and Networks_. Cambridge University Press, 2023.
  * [44] H. Kuroda and D. Kitahara, “Graph-structured sparse regularization via convex optimization,” in _Proc. IEEE ICASSP_ , 2022, pp. 5538–5542.
  * [45] A. Chambolle and T. Pock, “A first-order primal-dual algorithm for convex problems with applications to imaging,” _J. Math. Imaging Vis._ , vol. 40, no. 1, pp. 120–145, 2011.
  * [46] D. Gabay and B. Mercier, “A dual algorithm for the solution of nonlinear variational problems via finite element approximation,” _Comput. Math. Appl._ , vol. 2, no. 1, pp. 17–40, 1976.
  * [47] J. Eckstein and D. P. Bertsekas, “On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators,” _Math. Programm._ , vol. 55, no. 1, pp. 293–318, Apr. 1992.
  * [48] P. L. Combettes and J.-C. Pesquet, “Proximal splitting methods in signal processing,” in _Fixed-Point Algorithms for Inverse Problems in Science and Engineering_ , H. H. Bauschke, R. S. Burachik, P. L. Combettes, V. Elser, D. R. Luke, and H. Wolkowicz, Eds. Springer, 2011, pp. 185–212.
  * [49] L. Condat, D. Kitahara, A. Contreras, and A. Hirabayashi, “Proximal splitting algorithms for convex optimization: A tour of recent advances, with new twists,” _SIAM Rev._ , vol. 65, no. 2, pp. 375–435, 2023.
  * [50] P. L. Combettes, “Perspective functions: Properties, constructions, and examples,” _Set-Valued and Variational Analysis_ , vol. 26, no. 2, pp. 247–264, 2018.
  * [51] P. L. Combettes and C. L. Müller, “Perspective functions: Proximal calculus and applications in high-dimensional statistics,” _J. Math. Anal. Appl._ , vol. 457, no. 2, pp. 1283–1306, 2018.
  * [52] ——, “Perspective maximum likelihood-type estimation via proximal decomposition,” _Electron. J. Statist._ , vol. 14, no. 1, pp. 207–238, 2020\.
  * [53] H. H. Bauschke, M. K. Lal, and X. Wang, “Real roots of real cubics and optimization,” preprint arXiv:2302.10731, 2023.
  * [54] L. Condat, “Fast projection onto the simplex and the $\ell_{1}$ ball,” _Math. Programm._ , vol. 158, no. 1, pp. 575–585, 2016.
  * [55] G. H. Golub and C. F. Van Loan, _Matrix Computations_ , 4th ed. Johns Hopkins University Press, 2013.