Gridless Multisnapshot Variational Line Spectral Estimation from Coarsely Quantized Samples

TL;DR

This paper introduces a low-complexity, gridless variational line spectral estimation method combined with expectation propagation for estimating frequencies from coarsely quantized samples, improving performance in array signal processing.

Contribution

It proposes the MVALSE-EP algorithm that effectively estimates frequencies from low-resolution data, incorporating a CRB benchmark and analyzing parameter effects.

Findings

Snapshots improve frequency estimation accuracy.

The proposed method outperforms existing techniques in coarsely quantized scenarios.

Numerical experiments validate the effectiveness of MVALSE-EP, including real data.

Abstract

Due to the increasing demand for low power and higher sampling rates, low resolution quantization for data acquisition has drawn great attention recently. Consequently, line spectral estimation (LSE) with multiple measurement vectors (MMVs) from coarsely quantized samples is of vital importance in cutting edge array signal processing applications such as range estimation and DOA estimation in millimeter wave radar systems. In this paper, we combine the low complexity gridless variational line spectral estimation (VALSE) and expectation propagation (EP) and propose an MVALSE-EP algorithm to estimate the frequencies from coarsely quantized samples. In addition, the Cram\'{e}r Rao bound (CRB) is derived as a benchmark performance of the proposed algorithm, and insights are provided to reveal the effects of system parameters on estimation performance. It is shown that snapshots benefits the frequency estimation, especially in coarsely quantized scenarios. Numerical experiments are conducted to demonstrate the effectiveness of MVALSE-EP, including real data set.

Tables

Table 1. TABLE I: The range estimation performance of the MVALSE-EP for experiment 1.

B \T	T=1 (range (m), Amp. (dB))	T=16 (range (m), Amp. (dB))
1 bit	(2.32,)	(2.32,)
12 bit	(2.32,46.9),(0.08,23.6), (2.46, 26.4)	(2.32, 46.6), (0.09, 23.3),(2.44, 27.3)

Table 2. TABLE II: The range estimation performance of the MVALSE-EP for experiment 2.

B \T	T=1 (range (m), Amp. (dB))	T=16 (range (m), Amp. (dB))
1 bit	(2.65,0),(2.21,-6.21)	(2.66,0),(2.21,-5.99)
12 bit	(2.61,50.0),(2.31,46.4), (0.08, 23.7),(1.31, 22.6)	(2.61, 50.0), (2.31, 46.3),(0.11, 24.0)

Table 3. TABLE III: The range estimation performance of the MVALSE-EP for experiment 3.

B \T	T=1 (range (m), Amp. (dB))	T=16 (range (m), Amp. (dB))
1 bit	(2.29,0)	(2.29,0),(3.97,-6.68)
12 bit	(2.30,48.6),(3.97,43.3), (0.12, 22.9),(5.44, 23.0)	(2.30, 48.6), (3.97, 43.2),(0.12, 23.0),(5.41, 22.7)

Table 4. TABLE IV: The range estimation performance of the MVALSE-EP for experiment 4.

B \T	T=1 (range (m), Amp. (dB))	T=16 (range (m), Amp. (dB))
1 bit	(3.42,0),(2.35,-4.9),(0.13,-5.7)	(3.42,0),(2.40,-6.58)
12 bit	(3.46,27.4),(2.32,21.2), (3.77, 21.7),(0.12, 17.8)	(3.46, 27.2), (3.78, 21.7),(2.36, 20.4),(0.12, 18.4)

Table 5. TABLE V: The range estimation performance of the MVALSE-EP for experiment 5.

B \T	T=1 (range (m), Amp. (dB))	T=16 (range (m), Amp. (dB))
1 bit	(3.55,0),(2.75,-1.7)	(3.55,0),(2.74,-1.7)
12 bit	(3.56,30.7),(2.77,28.1), (0.08, 21.5),(4.02, 18.8)	(3.56, 30.5), (2.76, 28.8),(0.08, 22.0),(4.00, 16.9)

Equations

$$\mathbf{r}(t)=\mathbf{A}(\mathbf{\omega })\mathbf{x}(t)+\mathbf{w}(t),t=1,2,\cdots ,T,$$

$$\mathbf{A}=[\mathbf{a}({\omega }_1),\mathbf{a}({\omega }_2),\cdots ,\mathbf{a}({\omega }_k),\cdots ,\mathbf{a}({\omega }_K)]\in {\mathbb{C}}^{N\times K}$$

$$\mathbf{y}(t)=\mathcal{Q}\left\{\Re \left\{\mathbf{r}(t)\right\}\right\}+\mathrm{j}\mathcal{Q}\left\{\Im \left\{\mathbf{r}(t)\right\}\right\},t=1,2,\cdots ,T,$$

$$\mathcal{Q}(a)=o_l.$$

$${\mathrm{SNR}}_k=10{\log }_{10}\frac{\sum _{t=1}^T|x_k(t){|}^2/T}{{\sigma }^2},$$

$$Q(a)=o_l,\mathrm{if}a\in [{\tau }_l,{\tau }_{l+1}).$$

$$\mathbf{I}(\mathbf{\kappa })=\frac{2}{{\sigma }^2}{\left[\frac{\partial \mathbf{\mu }(\mathbf{\kappa })}{\partial {\mathbf{\kappa }}^{\mathrm{T}}}\right]}^{\mathrm{T}}\mathbf{\Lambda }\left[\frac{\partial \mathbf{\mu }(\mathbf{\kappa })}{\partial {\mathbf{\kappa }}^{\mathrm{T}}}\right],$$

$$\displaystyle\frac{\partial{\bm{\mu}}({\bm{\kappa}})}{\partial{\bm{\kappa}}^{\rm T}}=\left[\begin{array}[]{cccc}\frac{\partial[{\bm{\mu}}({\bm{\kappa}})]}{\partial\kappa_{1}}&\frac{\partial[{\bm{\mu}}({\bm{\kappa}})]}{\partial\kappa_{2}}&\cdots&\frac{\partial[{\bm{\mu}}({\bm{\kappa}})]}{\partial\kappa_{P}}\end{array}\right]\in{\mathbb{R}}^{N\times P},$$

$${\Lambda }_{i,i}=h({\mu }_i(\mathbf{\kappa }),{\sigma }^2),$$

$$h(x,{\sigma }^2)=\sum _{d=0}^{|\mathcal{D}|-1}\frac{{\left[\varphi \left(\frac{{\tau }_{d+1}-x}{\sigma /\sqrt{2}}\right)-\varphi \left(\frac{{\tau }_d-x}{\sigma /\sqrt{2}}\right)\right]}^2}{\Phi \left(\frac{{\tau }_{d+1}-x}{\sigma /\sqrt{2}}\right)-\Phi \left(\frac{{\tau }_d-x}{\sigma /\sqrt{2}}\right)}.$$

$$h(x,{\sigma }^2)=\frac{{\varphi }^2\left(\frac{x}{\sigma /\sqrt{2}}\right)}{\Phi \left(\frac{x}{\sigma /\sqrt{2}}\right)\Phi \left(\frac{-x}{\sigma /\sqrt{2}}\right)}=\frac{1}{2\pi }\frac{{\mathrm{e}}^{-\frac{2x^2}{{\sigma }^2}}}{\Phi \left(\frac{x}{\sigma /\sqrt{2}}\right)\Phi \left(\frac{-x}{\sigma /\sqrt{2}}\right)}.$$

$$\displaystyle{\bm{\kappa}}=\left[\begin{array}[]{c}\Re\{{\mathbf{x}}(1)\}\\ \Im\{{\mathbf{x}}(1)\}\\ \Re\{{\mathbf{x}}(2)\}\\ \Im\{{\mathbf{x}}(2)\}\\ \vdots\\ \vdots\\ \Re\{{\mathbf{x}}(T)\}\\ \Im\{{\mathbf{x}}(T)\}\\ {\bm{\omega}}\end{array}\right],{\bm{\mu}}({\bm{\kappa}})=\left[\begin{array}[]{c}\Re\{{\mathbf{A}}{\mathbf{x}}(1)\}\\ \Im\{{\mathbf{A}}{\mathbf{x}}(1)\}\\ \Re\{{\mathbf{A}}{\mathbf{x}}(2)\}\\ \Im\{{\mathbf{A}}{\mathbf{x}}(2)\}\\ \vdots\\ \vdots\\ \Re\{{\mathbf{A}}{\mathbf{x}}(T)\}\\ \Im\{{\mathbf{A}}{\mathbf{x}}(T)\}\end{array}\right].$$

$$\bar{\mathbf{A}}$$

$${\mathbf{D}}_{\omega }=\left[\frac{\mathrm{d}\mathbf{a}({\omega }_1)}{\mathrm{d}{\omega }_1},\frac{\mathrm{d}\mathbf{a}({\omega }_2)}{\mathrm{d}{\omega }_2},\cdots ,\frac{\mathrm{d}\mathbf{a}({\omega }_K)}{\mathrm{d}{\omega }_K}\right],$$

$$\mathbf{X}(t)\triangleq \mathrm{diag}(\mathbf{x}(t)),$$

$$\displaystyle\bar{\mathbf{D}}(t)=\left[\begin{array}[]{cc}\frac{\partial{\Re\{{\mathbf{Ax}}(t)\}}}{\partial{\bm{\omega}}}\\ \frac{\partial{\Im\{{\mathbf{Ax}}(t)\}}}{\partial{\bm{\omega}}}\end{array}\right]=\left[\begin{array}[]{cc}\Re\{{\mathbf{D}}_{\omega}{\mathbf{X}}(t)\}\\ \Im\{{\mathbf{D}}_{\omega}{\mathbf{X}}(t)\}\end{array}\right],$$

$$\displaystyle{\bm{\lambda}}(t)=\left[\begin{array}[]{c}h(\Re\{{\mathbf{A}}{\mathbf{x}}(t)\},{\sigma}^{2})\\ h(\Im\{{\mathbf{A}}{\mathbf{x}}(t)\},{\sigma}^{2})\end{array}\right],$$

$$\mathbf{\Lambda }(t)=\mathrm{diag}(\mathbf{\lambda }(t)),$$

$$\mathbf{\Lambda }=\left(\begin{array}{cccc}\mathbf{\Lambda }(1)& & & \\ & \mathbf{\Lambda }(2)& & \\ & & \ddots & \\ & & & \mathbf{\Lambda }(T)\end{array}\right),$$

$$\bar{\mathbf{H}}(t)={\bar{\mathbf{A}}}^{\mathrm{T}}\mathbf{\Lambda }(t)\bar{\mathbf{A}},$$

$$\bar{\mathbf{\Delta }}(t)={\bar{\mathbf{A}}}^{\mathrm{T}}\mathbf{\Lambda }(t)\bar{\mathbf{D}}(t).$$

$$\mathrm{CRB}(\mathbf{\omega })={\left(\frac{2}{{\sigma }^2}\sum _{t=1}^T\left({\bar{\mathbf{D}}}^{\mathrm{T}}(t)\mathbf{\Lambda }(t)\bar{\mathbf{D}}(t)-{\bar{\mathbf{\Delta }}}^{\mathrm{T}}(t){\bar{\mathbf{H}}}^{-1}(t)\bar{\mathbf{\Delta }}(t)\right)\right)}^{-1}.$$

$$\mathrm{CRB}(\mathbf{\theta })=\mathrm{diag}\left({\left(\frac{\lambda }{2\pi d\cos (\mathbf{\theta })}\right)}^2\right)\mathrm{CRB}(\mathbf{\omega }).$$

$$\displaystyle{\bm{\kappa}}=\left[\begin{array}[]{c}\psi(1)\\ \psi(2)\\ \vdots\\ \psi(T-1)\\ \psi(T)\\ \omega\end{array}\right],{\bm{\mu}}({\bm{\kappa}})=\left[\begin{array}[]{c}g\cos(\omega{\mathbf{n}}+\psi(1))\\ g\sin(\omega{\mathbf{n}}+\psi(1))\\ g\cos(\omega{\mathbf{n}}+\psi(2))\\ g\sin(\omega{\mathbf{n}}+\psi(2))\\ \vdots\\ \vdots\\ g\cos(\omega{\mathbf{n}}+\psi(T))\\ g\sin(\omega{\mathbf{n}}+\psi(T))\\ \end{array}\right].$$

$$\frac{\partial \mathbf{\mu }(\mathbf{\kappa })}{\partial {\mathbf{\kappa }}^{\mathrm{T}}}=\left(\begin{array}{ccccc}\bar{\mathbf{a}}(1)& & & & \bar{\mathbf{d}}(1)\\ & \bar{\mathbf{a}}(2)& & & \bar{\mathbf{d}}(2)\\ & & \ddots & & ⋮\\ & & & \bar{\mathbf{a}}(T)& \bar{\mathbf{d}}(T)\end{array}\right)\in {\mathbb{R}}^{2NT\times (T+1)}$$

$$\displaystyle\bar{\mathbf{a}}(t)=\left[\begin{array}[]{cc}-g\sin({\mathbf{n}}\omega+\psi(t))\\ g\cos({\mathbf{n}}\omega+\psi(t))\end{array}\right],$$

$$\displaystyle\bar{\mathbf{d}}(t)=\left[\begin{array}[]{cc}\frac{\partial{\Re\{{\mathbf{ax}}(t)\}}}{\partial{\bm{\omega}}}\\ \frac{\partial{\Im\{{\mathbf{ax}}(t)\}}}{\partial{\bm{\omega}}}\end{array}\right]=g\left[\begin{array}[]{cc}-{\mathbf{n}}\odot\sin({\mathbf{n}}\omega+\psi(t))\\ {\mathbf{n}}\odot\cos({\mathbf{n}}\omega+\psi(t))\end{array}\right],$$

$$\displaystyle{\bm{\lambda}}(t)=\left[\begin{array}[]{cc}h(g\cos({\mathbf{n}}\omega+\psi(t)),\sigma^{2})\\ h(g\sin({\mathbf{n}}\omega+\psi(t)),\sigma^{2})\end{array}\right]$$

$${\mathrm{CRB}}^{-1}(\omega )=\frac{2}{{\sigma }^2}\sum _{t=1}^T\left({\bar{\mathbf{d}}}^{\mathrm{T}}(t)\mathbf{\Lambda }(t)\bar{\mathbf{d}}(t)-\frac{{\bar{\Delta }}^2(t)}{\bar{H}(t)}\right),$$

$$\bar{\Delta }(t)={\bar{\mathbf{a}}}^{\mathrm{T}}(t)\mathrm{diag}(\mathbf{\lambda }(t))\bar{\mathbf{d}}(t)=g^2\sum _{n=0}^{N-1}\chi \left(n,\omega ,\psi ,\frac{g}{\sigma }\right)n,$$

$$\bar{H}(t)={\bar{\mathbf{a}}}^{\mathrm{T}}(t)\mathrm{diag}(\mathbf{\lambda }(t))\bar{\mathbf{a}}(t)=g^2\sum _{n=0}^{N-1}\chi \left(n,\omega ,\psi ,\frac{g}{\sigma }\right),$$

$$\chi (n,\omega ,\psi ,\frac{g}{\sigma })={\sin }^2(n\omega +\varphi (t))h(g\cos (n\omega +\psi (t)),{\sigma }^2)+{\cos }^2(n\omega +\psi (t))h(g\sin (n\omega +\psi (t)),{\sigma }^2).$$

$${\mathrm{ECRB}}^{-1}(\omega )=\frac{2}{{\sigma }^2}T{\mathrm{E}}_{\psi }\left[{\bar{\mathbf{d}}}^{\mathrm{T}}\mathbf{\Lambda }\bar{\mathbf{d}}-\frac{{\bar{\Delta }}^2}{\bar{H}}\right].$$

$${\mathrm{ECRB}}_{\infty }(\omega )=\frac{{\sigma }^2}{|g{|}^2}\frac{6}{N(N-1)(N+1)}\frac{1}{T}.$$

$${\mathrm{E}}_{\varphi }\left[\chi (n,\omega ,\psi ,\frac{g}{\sigma })\right]={\mathrm{E}}_{\psi }\left[{\sin }^2(\psi )h(g\cos (\psi ),{\sigma }^2)+{\cos }^2(\psi )h(g\sin (\psi ),{\sigma }^2)\right]\triangleq r\left(\frac{g}{\sigma }\right)$$

Full text

keywords: line spectral estimation, multiple measurement vectors, quantization, variational Bayes, expectation propagation

I Introduction

As radar and communication systems scale up in bandwidth, the cost and power consumption of high-precision (e.g., $10-12$ bits) analog-to-digital converters (ADCs) become the limiting factor for various applications, such as cognitive radio and radar. As a remedy, low resolution quantization has been gradually drawn attention in recent years [1, 2, 3, 4, 5, 6, 7, 8]. In addition, in military radar systems incorporating very large arrays, low resolution quantization is a promising technique as it provides a reasonable estimation accuracy, and it also enables a very high sampling rate for a wideband system. Moreover, it dramatically reduces the memory and transmission requirements for the samples.

For low resolution quantization, the loss of information per sample is large, and new, low-resolution quantization adapted digital signal processing algorithms should be developed [9, 10]. From the linear signal processing point of view, the spectrum of 1-bit line spectral signal is analyzed. As shown in [17, 6], the binary data contains plentiful self-generated [17] and cross-generated harmonics [6]. For low signal to noise (SNR) scenario, the strengths of harmonics decay quickly and conventional fast Fourier transform (FFT) performs well. While for high SNR scenario, FFT will overestimate the model order. From the nonlinear signal processing point of view, the direction of arrival (DOA) estimation from 1 bit data is treated as a nonlinear parameter estimation problem. In [18], the frequency estimation from 1 bit quantized samples has been addressed. The problems in [18] adopts the deterministic signal modeling by treating the signal samples as a deterministic, yet unknown set of parameters that has to be estimated. By comparing the derived Cramèr-Rao bound (CRB) for 1-bit quantization with that of infinite quantization, it is found that 1 bit quantization gives a dramatic increase of variance at certain frequencies, and a slightly worse performance for other frequencies. In [19], a single stochastic Gaussian point source model is assumed to study the DOA estimation which is a spatial analogy to the problem of temporal line spectral estimation (LSE) with multiple measurement vectors (MMVs), and the CRB for a two-sensor array case is derived. It is shown that the estimation error has weak dependency on signal-to-noise ratio (SNR), and there exist two singular DOA angles $0^{\circ}$ and $30^{\circ}$ for which higher SNR results in better estimation performance.

In this work, we present an approach named multisnapshot VALSE-EP (MVALSE-EP) for LSE with MMVs from coarsely quantized measurements. Compared to [19], we use a deterministic signal model by treating the amplitudes of the signal as unknown deterministic parameters, instead of stochastic Gaussian model. For the stochastic signal model, the spatial signal samples are clearly dependent (as the sensors are correlated), and the covariance matrix is not diagonal. In order to apply the arcsine law for LSE, Gaussian assumption for the signal amplitude must be adopted. This implies that when the signal amplitude distribution violates the Gaussian assumption, the estimation performance of the covariance matrix based approach may degrade significantly. Besides, the number of snapshots is usually needed to be large to obtain accurate reconstruction of the covariance matrix, e.g., $10^{4}$. By using the deterministic signal model, the signal samples can be regarded as independent identically distributed (i.i.d.) variables, and consequently, performing approximate Bayesian inference becomes tractable [13, 14, 16], as we show later. Compared to [16] which proposes the VALSE-EP to tackle the LSE in the single snapshot case, we extend it to address the multisnapshot LSEs. This is meaningful especially in coarsely quantized scenarios as model may become unidentifiable in the single snapshot case in certain scenarios [21], while multisnapshots makes the estimation problem unlikely identifiable and benefits the estimation performance [22].

In this paper, the LSE estimation with MMVs from coarsely quantized samples is studied. Firstly, to provide a bench mark performance of the LSE, the CRB is derived. It is shown that under $1$ bit quantization, the CRB is inversely proportional to the number of snapshots and the cubic of the number of antennas. For lower SNR scenario, the CRB is inversely proportional to SNR, while the CRB is inversely proportional to the square root of the SNR for high SNR scenario. Secondly, the MVALSE-EP is proposed to estimate the DOAs and the complex weight. Since MVALSE-EP performs Newton step to refine the frequencies, it overcomes the model mismatch [23] issue incurred by on-grid assumptions [24]. The expectation maximization (EM) algorithm is incorporated to estimate the noise variance automatically. Finally, numerical experiments are conducted to demonstrate the effectiveness of MVALSE-EP, including real data set.

The rest of this paper is organized as below. Section II describes the system model. Section III derives CRB. The MVALSE-EP algorithm and the details of the updating expressions are presented in Section IV. Substantial numerical experiments are provided in Section V and Section VI concludes the paper.

For a complex vector ${\mathbf{x}}\in{\mathbb{C}}^{M}$, let ${\Re}\{{\mathbf{x}}\}$ and $\Im\{\mathbf{x}\}$ denote the real and imaginary part of $\mathbf{x}$, respectively, let $|\mathbf{x}|$ and $\angle{\mathbf{x}}$ denote the componentwise amplitude and phase of $\mathbf{x}$, respectively. For the square matrix $\mathbf{A}$, let ${\rm diag}(\mathbf{A})$ return a vector with elements being the diagonal of $\mathbf{A}$. While for a vector ${\mathbf{a}}$, let ${\rm diag}(\mathbf{a})$ return a diagonal matrix with the diagonal being $\mathbf{a}$, and thus ${\rm diag}({\rm diag}(\mathbf{A}))$ returns a diagonal matrix. Let ${\rm j}$ denote the imaginary number. Let ${\mathcal{S}}\subset\{1,\cdots,N\}$ be a subset of indices and $|{\mathcal{S}}|$ denote its cardinality. For the matrix $\mathbf{J}\in\mathbb{C}^{N\times N}$, let $\mathbf{J}_{\mathcal{S}}$ denote the submatrix by choosing both the rows and columns of $\mathbf{J}$ indexed by $\mathcal{S}$. Similarly, let $\mathbf{h}_{\mathcal{S}}$ denote the subvector by choosing the elements of $\mathbf{h}$ indexed by $\mathcal{S}$. Let ${(\cdot)}^{*}_{\mathcal{S}}$, ${(\cdot)}^{\rm T}_{\mathcal{S}}$ and ${(\cdot)}^{\rm H}_{\mathcal{S}}$ be the conjugate, transpose and Hermitian transpose operator of ${(\cdot)}_{\mathcal{S}}$, respectively. For the matrix $\mathbf{A}$, let $|\mathbf{A}|$ denote the elementwise absolute value of $\mathbf{A}$. Let $\mathbf{I}_{L}$ denote the identity matrix of dimension $L$. $``\sim i"$ denotes the indices ${\mathcal{S}}$ excluding $i$. Let ${\mathcal{CN}}({\mathbf{x}};{\bm{\mu}},{\bm{\Sigma}})$ denote the complex normal (CN) distribution of ${\mathbf{x}}$ with mean ${\bm{\mu}}$ and covariance ${\bm{\Sigma}}$. For a frequency $\theta$, let ${\mathcal{V}M}(\theta;\mu,\kappa)$ be the von Mises distribution with $\mu$ and $\kappa$ being the mean and concentration parameters. For a random vector $\mathbf{x}$ with probability density function (PDF) $p({\mathbf{x}})$, let ${\rm Proj}[p({\mathbf{x}})]$ denote the projection of $p({\mathbf{x}})$ onto Gaussian PDF with diagonal covariance matrix, where the means and variances are matched with that of $p({\mathbf{x}})$. Let $\phi(x)=\exp(-x^{2}/2)/{\sqrt{2\pi}}$ and $\Phi(x)=\int_{-\infty}^{x}\phi(t){\rm d}t$ denote the standard normal probability density function (PDF) and cumulative distribution function (CDF), respectively.

II Problem Setup

Consider a line spectral estimation problem with MMVs. For the $t$th snapshot, the signal received before quantization can be modeled as

[TABLE]

where ${\mathbf{x}}(t)=[x_{1}(t),x_{2}(t),\cdots,x_{K}(t)]^{\rm T}$ denotes the complex coefficients of the unknown frequencies ${\bm{\omega}}=[\omega_{1},\omega_{2},\cdots,\omega_{K}]^{\rm T}\in{\mathbb{R}}^{K}$, $\omega_{i}\in[0,2\pi)$, and $K$ denotes the number of frequencies,

[TABLE]

denotes the nominal array manifold matrix, with the steering vectors ${\mathbf{a}}(\omega_{k})=[1,{\rm e}^{{\rm j}\omega_{k}},\cdots,{\rm e}^{{\rm j}(N-1)\omega_{k}}]^{\rm T}$ as its columns, ${\mathbf{w}}(t)\in{\mathbb{C}}^{M}$ denotes the additive internal (e.g., thermal) receiver noise modeled as spatially and temporally independent, identically distributed (i.i.d.) zero-mean circular CN with a covariance matrix ${\mathbf{R}}_{w}\triangleq{\rm E}[{\mathbf{w}}(t){\mathbf{w}}^{\rm H}(t)]=\sigma^{2}{\mathbf{I}}_{N}$ where $\sigma^{2}$ is unknown.

The signal ${\mathbf{r}}(t)$ is then coarsely quantized to obtain

[TABLE]

where ${\mathcal{Q}}(\cdot)$ is a quantizer which maps the continuous-valued observations into a finite number of bits. Note that for a quantizer with bit-depth $B$, the cardinality of the output of the quantizer is $|{\mathcal{D}}|=2^{B}$. Assume that the quantization intervals for the quantizer ${\mathcal{Q}}(\cdot)$ are $\{(\tau_{l},\tau_{l+1})\}_{l=0}^{|{\mathcal{D}}|-1}$, where $\tau_{0}=-\infty$, $\tau_{|{\mathcal{D}}|}=\infty$, $\bigcup_{l=0}^{|{\mathcal{D}}|-1}[\tau_{l},\tau_{l+1})={\mathbb{R}}$. Given a real number $a\in[\tau_{l},\tau_{l+1})$ input, the quantizer output is

[TABLE]

For example, one-bit quantization refers to $B=1$, $D=2$, $\tau_{0}=-\infty$, $\tau_{1}=0$ and $\tau_{2}=\infty$, ${\mathcal{Q}}(\cdot)$ reduces to the signum function, i.e., $Q(\cdot)={\rm sign}(\cdot)$.

The SNR of the $k$th target is defined as

[TABLE]

which corresponds to the sample SNR. Note that after coherent integration, the sample SNR can be improved up to $10\log_{10}N$ dB.

It is worth noting that model (1) and model (3) has various applications such as range, velocity and angle estimation in millimeter wave LFMCW radar [25]. Below we present the details.

• •

range estimation: here $\omega_{k}=2\pi\times\frac{2\mu T_{s}}{c}r_{k}$ where $\mu$, $T_{s}$, $c$ and $r_{k}$ denote the chirp rate, sampling interval, electromagnetic speed and radial distance of the $k$th target.

• •

velocity estimation: here $\omega_{k}=2\pi\times\frac{2T_{c}}{\lambda}v_{k}$ where $T_{c}$, $\lambda$ and $v_{k}$ denote the repetition interval, wavelength and radial velocity of the $k$th target.

• •

DOA estimation: here $\omega_{k}=2\pi\times\frac{d}{\lambda}\sin\theta_{k}$ where $d$, $\lambda$ and $\theta_{k}$, denote the element spacing, wavelength and direction of arrival of the $k$th target.

• •

vital signs detection [26]: The vital sign extraction from multiple independent self-injection locking doppler radars can be formulated as a line spectral estimation problem with MMVs, see [26, eq. (5)].

III Cramér Rao bound

The CRB is a lower bound of unbiased estimators and provides a benchmark against which we can compare the performance of the proposed algorithm [27, 28]. To derive the CRB, $K$ is assumed to be known, the frequencies ${\bm{\omega}}\in{\mathbb{R}}^{K}$ and weights ${\mathbf{X}}=[{\mathbf{x}}(1),\cdots,{\mathbf{x}}(T)]\in{C}^{K\times T}$ are treated as unknown deterministic parameters. As for the quantizer $Q(\cdot)$, the quantization intervals are $\{(\tau_{l},\tau_{l+1})\}_{l=0}^{|{\mathcal{D}}|-1}$, where $\tau_{0}=-\infty$, $\tau_{{\mathcal{D}}}=\infty$, $\bigcup_{l=0}^{{\mathcal{D}}-1}[\tau_{l},\tau_{l+1})={\mathbb{R}}$. Given a real number $a\in[\tau_{l},\tau_{l+1})$, the representation is

[TABLE]

Note that for a quantizer with bit-depth $B$, the cardinality of the output of the quantizer is $|{\mathcal{D}}|=2^{B}$.

The CRB equals to the inverse of the Fisher information matrix (FIM). For calculating the FIM, the following lemma can be utilized.

Lemma 1

[15, 16]** Let ${\bm{\kappa}}\in{\mathbb{R}}^{P}$ denote the set of unknown deterministic parameters. Note that in the case of quantized observations ${\mathbf{y}}={\mathcal{Q}}({\mathbf{r}})\in{\mathbb{R}}^{N}$ where ${\mathbf{r}}\sim{\mathcal{N}}({\bm{\mu}}({\bm{\kappa}}),{\sigma}^{2}{\mathbf{I}}_{N}/2)$, the FIM is given by

[TABLE]

where

[TABLE]

and ${\bm{\Lambda}}$ is a diagonal matrix with the $(i,i)$th element

[TABLE]

and $h(x,\sigma^{2})$ is

[TABLE]

For one-bit quantization where $|{\mathcal{D}}|=2$, $\tau_{0}=-\infty$, $\tau_{1}=0$, $\tau_{2}=\infty$, $h(x,\sigma^{2})$ simplifies to be

[TABLE]

* *

For unquantized system, the FIM (7) is obtained with ${\bm{\Lambda}}={\mathbf{I}}_{N}$.

In our setting, the observations are $\left[\Re\{{\mathbf{y}}(1)\};\Im\{{\mathbf{y}}(1)\};\Re\{{\mathbf{y}}(2)\};\Im\{{\mathbf{y}}(2)\};\cdots;\Re\{{\mathbf{y}}(T)\};\Im\{{\mathbf{y}}(T)\}\right]$. Note that ${\bm{\kappa}}\in{\mathbb{R}}^{2KT+K}$ and ${\bm{\mu}}({\bm{\kappa}})\in{\mathbb{R}}^{2MT}$ are

[TABLE]

Define

[TABLE]

As shown in Appendix VII-A, according to Lemma 1, the ${\rm CRB}({\bm{\omega}})$ is

[TABLE]

Note that for the DOA estimation problem where $\omega=\frac{2\pi d}{\lambda}\sin\theta$ with $d$ and $\lambda$ being the element spacing and wavelength, by using vector parameter CRB for transformations [27, eq. (3.30)], the CRB of DOAs is

[TABLE]

We now hope to provide insight into the relationship between the system parameters and the estimation performance. As a result, we provide an asymptotic analysis of the CRB of a single target under one bit quantization for both low SNR and high SNR scenarios. It is assumed that the amplitude of the frequency $g=|x(t)|$ is known and is the same for all the snapshots. Besides, the phase $\psi(t)=\angle x(t)$ of the frequency at the $t$th snapshot is unknown. We assume that $\psi(t)$ is uniformly drawn from $[0,2\pi)$. Define ${\mathbf{n}}=[0,1,\cdots,N-1]^{\rm T}$. Similar to the derivation of the CRB, here the unknown deterministic parameters ${\bm{\kappa}}$ and ${\bm{\mu}}({\bm{\kappa}})$ are

[TABLE]

Then

[TABLE]

where $\bar{\mathbf{a}}(t)$ and $\bar{\mathbf{d}}(t)$ are defined as

[TABLE]

respectively, $\odot$ denotes the Hardarmard product, ${\bm{\lambda}}(t)$ is

[TABLE]

with $h(\cdot,\cdot)$ being given by (12). It can be shown that the CRB for a single frequency is

[TABLE]

where $\bar{\Delta}(t)$ and $\bar{H}(t)$ are

[TABLE]

where

[TABLE]

Now suppose that $\psi(t)\sim{\mathcal{U}}(0,2\pi)$. We define the expected ${\rm CRB}^{-1}({\omega})$ as

[TABLE]

It can be seen that given that the number of snapshots is large, ${\rm ECRB}^{-1}({\omega})$ (65) provides a good approximation for ${\rm CRB}^{-1}({\omega})$ (61). Note that for unquantized measurements, one has $\chi(n,\omega,\phi,\frac{g}{\sigma})=1$. It is straightforward to show that ${\rm ECRB}_{\infty}({\omega})$ for a single frequency under unquantized measurements is

[TABLE]

Another property is that the integration

[TABLE]

holds, i,e, the integration is independent of $n\omega$.

Note that

[TABLE]

It is hard to directly calculate $\frac{2}{\sigma^{2}}{\rm E}_{\psi}\left[\frac{\bar{\Delta}^{2}}{\bar{H}}\right]$. Instead we use the Jensen inequality to obtain a bound

[TABLE]

Here

[TABLE]

Furthermore, we have

[TABLE]

Thus

[TABLE]

According to (68), (69) and (72), ${\rm ECRB}^{-1}({\omega})$ (65) can be approximated as

[TABLE]

due to $\sum\limits_{n=0}^{N-1}n^{2}=\frac{(N-1)N(2N-1)}{6}$ and $\sum\limits_{n=0}^{N-1}n=\frac{N(N-1)}{2}$.

The Chernoff bound [30] $\Phi(x)\Phi(-x)\leq\frac{1}{4}{\rm e}^{-\frac{x^{2}}{2}}$ can be used to approximate $h(x;\sigma^{2})$ as

[TABLE]

which is also a very tight approximation. Utilizing the integration [29]

[TABLE]

where $I_{v}(\cdot)$ are the modified Bessel functions of the first kind, $r\left(\frac{g}{\sigma}\right)$ can be approximated as

[TABLE]

where $\stackrel{{\scriptstyle a}}{{\leq}}$ is due to (74). Substituting (III) in (III), ${\rm ECRB}^{-1}({\omega})$ can be further approximated as

[TABLE]

Therefore

[TABLE]

i.e., ${\rm ELCRB}({\omega})$ is a lower bound of ${\rm ECRB}({\omega})$ and approximates ${\rm ECRB}({\omega})$ well.

For low SNR scenario, $I_{0}(x)\approx 1$ and $I_{1}(x)\approx 0$, ${\rm ELCRB}^{-1}({\omega})$ (III) is approximated as

[TABLE]

and

[TABLE]

For high SNR scenario, one has

[TABLE]

Substituting (81) into (III) yields

[TABLE]

According to (79) and (82), we have

[TABLE]

Because ${\rm ELCRB}({\omega})$ approximates ${\rm CRB}({\omega})$ well, it can be concluded that under one-bit quantization, the CRB is inversely proportional to the number of snapshots and the cubic of the number of antennas. For lower SNR scenario, the CRB is inversely proportional to SNR, while the CRB is inversely proportional to the square root of the SNR for high SNR scenario.

IV MVALSE-EP Algorithm

In this section, MVALSE-EP algorithm is developed to estimate the frequencies. The key step in designing the MVALSE-EP is to introduce the suitable hidden variable and expand the factor graph. Here we define

[TABLE]

and the EP is adopted to iteratively approximate the nonlinear measurement model as a pseudo linear measurement model. In addition, EM and VALSE are adopted to estimate the nuisance parameters and the frequencies. In the following text, we first present the modeling setup. Then we adopt the modularized point of view to design MVALSE-EP algorithm based on existing superresolution VALSE approach.

IV-A Modeling Setup

Since the number of targets $K$ is usually unknown, an overcomplete model in which the signal consisting of $N$ targets is assumed [11]

[TABLE]

where ${\mathbf{A}}({\bm{\omega}})=[{\mathbf{a}}({\omega}_{1}),\cdots,{\mathbf{a}}({\omega}_{N})]$ and satisfies $N>K$. To model the unknown nature of $K$, the binary hidden variables ${\mathbf{s}}(t)$ are introduced, where $s_{k}(t)=1$ means that the $k$th target is active in the $t$th snapshot, otherwise deactive (${x}_{k}(t)=0$). In addition, it is required that the number of sources are the same across the snapshots, which requires that

[TABLE]

The probability mass function of $s_{k}(t)$ is

[TABLE]

Given that $s_{k}(t)=1$, we assume that ${x}_{k}(t)\sim{\mathcal{CN}}({x}_{k}(t);0,\tau)$. Thus $(s_{k}(t),{x}_{k}(t))$ follows a Bernoulli-Gaussian distribution, that is

[TABLE]

From (87) and (88), it can be seen that the parameter $\rho$ denotes the probability of the $k$th component being active and $\tau$ is a variance parameter. The variable ${\bm{\omega}}=[\omega_{1},...,\omega_{N}]^{\rm T}$ has the prior PDF $p({\bm{\omega}})=\begin{matrix}\prod_{k=1}^{N}p(\omega_{k})\end{matrix}$. Without any knowledge of the frequency $\omega$, the uninformative prior distribution $p(\omega_{k})={1}/({2\pi})$ is used [11]. For encoding the prior distribution, please refer to [11, 12] for further details. Note that the modeling setup of (88) is different from that of [12]. In [12], only binary variables $\{s_{k}\}_{k=1}^{N}$ are introduced, and such modeling makes the algorithm overestimate the model order.

Given $\mathbf{Z}$, the PDF $p({\mathbf{Y}}|{\mathbf{Z}};\sigma_{w}^{2})=\prod\limits_{t=1}^{T}p({\mathbf{y}}(t)|{\mathbf{z}}(t);\sigma_{w}^{2})=\prod\limits_{m=1}^{N}\prod\limits_{t=1}^{T}p(y_{m}(t)|z_{m}(t);\sigma_{w}^{2})$ of $\mathbf{Y}$ can be easily calculated through (3). Let

[TABLE]

be the set of all random variables and the model parameters, respectively. According to the Bayes rule, the joint PDF $p({\mathbf{Y}},{\mathbf{Z}},{\bm{\Omega}};{\bm{\beta}})$ is

[TABLE]

Given the above joint PDF (91), the type II maximum likelihood (ML) estimation of the model parameters $\hat{\bm{\beta}}_{\rm ML}$ is

[TABLE]

Then the minimum mean squared error (MMSE) estimates of the parameters $({\mathbf{Z}},{\bm{\Omega}})$ is

[TABLE]

where the expectation is taken with respect to

[TABLE]

Directly solving the ML estimate of $\bm{\beta}$ (92) or the MMSE estimate of $({\mathbf{Z}},{\bm{\Omega}})$ (93) are both intractable. As a result, an iterative algorithm is designed in the ensuing text.

IV-B Algorithm Design

Below we describe the details of the algorithm and show that how to utilize the basic VALSE algorithm to perform the DOAs estimation. For brevity, the iteration index is omitted. The factor graph of the joint PDF (91) is shown in Fig. 1, and the algorithm is designed according to Fig. 2.

Firstly, initialize the message transmitted from the factor node $\delta({\mathbf{Z}}-{\mathbf{A}}{\mathbf{X}})$ to the variable node ${\mathbf{Z}}$ as

[TABLE]

IV-B1 Componentwise MMSE

According to EP, the message $m_{{\mathbf{Z}}\rightarrow\delta}({\mathbf{Z}})$ transmitted from the variable node ${\mathbf{Z}}$ to the factor node $\delta({\mathbf{Z}}-{\mathbf{A}}{\mathbf{X}})$ can be calculated as [13]

[TABLE]

where $\propto$ denotes identity up to a normalizing constant. First, the MMSE estimate of ${\mathbf{z}}(t)$ can be obtained, i.e.,

[TABLE]

where ${\rm E}[\cdot|q_{{\rm B}}({\mathbf{z}}(t))]$ and ${\rm Var}[\cdot|q_{{\rm B}}({\mathbf{z}}(t))]$ are the mean and variance operations taken componentwise with respect to the distribution $\propto q_{{\rm B}}({\mathbf{z}}(t))$. Here we adopt the diagonal EP and ${\rm Proj}[q_{{\rm B}}({\mathbf{z}}(t))]$ is

[TABLE]

Substituting (99) in (96), the message $m_{{\mathbf{z}}\rightarrow\delta}({\mathbf{z}}(t))$ from the variable node ${\mathbf{z}}(t)$ to the factor node $\delta({\mathbf{z}}-{\mathbf{A}}({\bm{\omega}}){\mathbf{x}}(t))$ is calculated as

[TABLE]

where ${\mathbf{z}}_{{\rm B}}^{\rm ext}(t)$ and ${\mathbf{v}}_{{\rm B}}^{\rm ext}(t)$ are [14]

[TABLE]

where $\odot$ denotes componentwise multiplication. Consequently, we have

[TABLE]

IV-B2 MVALSE module

According to (102), the message $m_{{\mathbf{Z}}\rightarrow\delta}({\mathbf{Z}})$ transmitted from the variable node ${\mathbf{Z}}$ to the factor node $\delta({\mathbf{Z}}-{\mathbf{A}}{\mathbf{X}})$ is Gaussian distributed and is independent of the snapshot $t$. Based on the definition of the factor node $\delta({\mathbf{Z}}-{\mathbf{A}}{\mathbf{X}})$, $T$ pseudo linear observation models

[TABLE]

are obtained, where $\tilde{\mathbf{w}}(t)\sim{\mathcal{CN}}({\mathbf{0}},{\rm diag}(\tilde{\bm{\sigma}}^{2}(t)))$, $\tilde{\mathbf{y}}(t)={\mathbf{z}}_{{\rm B}}^{\rm ext}(t)$ and $\tilde{\bm{\sigma}}^{2}(t)={\mathbf{v}}_{{\rm B}}^{\rm ext}(t)$. For the $t$th equation in (103), the variances of the heteroscedastic noise $\tilde{\mathbf{w}}(t)$ are different. In addition, All the snapshots share the same frequency ${\bm{\omega}}$. Then we run the MVALSE under heteroscedastic noise algorithm [31] with known noise variances. Note that we could average the noise variance with respect to snapshots or measurements. Such operation simplifies the computation especially when averaging over snapshots. The only difference is that for this modeling setup, the $k$th target is active only when

[TABLE]

where $v_{k}(t)$ and $u_{k}(t)$ are

[TABLE]

$\bm{\Sigma}^{-1}(t)={\rm diag}\left(\tilde{\bm{\sigma}}^{2}(t)\right)$, ${\mathbf{J}}(t)$ and ${\mathbf{h}}(t)$ are

[TABLE]

where ${J}_{i,j}(t)$ denotes the $(i,j)$th element of ${\mathbf{J}}(t)$ and ${\widehat{\mathbf{a}}}_{j}$ denotes the posterior mean of ${\mathbf{a}}_{j}$. The deactive case can be obtained in a similar way. Once the posterior PDF $q({\bm{\omega}}|\tilde{\mathbf{Y}})$, $q({\mathbf{X}}|\tilde{\mathbf{Y}})$ and the model parameters $\tau$ and $\rho$ are obtained, the posterior means and variances of ${\mathbf{z}}(t)={\mathbf{A}}({\bm{\omega}}){\mathbf{x}}(t)$ are also obtained as [16]

[TABLE]

where $\hat{\mathbf{x}}_{\hat{S}}(t)$ and $\hat{\mathbf{C}}_{\hat{S}}(t)$ are the posterior means and covariance matrix of ${\mathbf{x}}(t)$, $\hat{\mathbf{A}}_{\hat{S}}$ is the estimate of ${\mathbf{A}}_{\hat{S}}$. The noise variance $\sigma^{2}$ can be obtained via the EM algorithm as

[TABLE]

Then we calculate the message ${m}_{\delta\rightarrow{\mathbf{z}}(t)}({\mathbf{z}}(t))$ as

[TABLE]

where ${\rm Proj}[q_{\rm A}({\mathbf{z}}(t))]$ is

[TABLE]

According to (110), ${m}_{\delta\rightarrow{\mathbf{z}}(t)}({\mathbf{z}}(t))$ is calculated to be

[TABLE]

where the extrinsic ${\mathbf{z}}_{{\rm A}}^{\rm ext}(t)$ and variance ${\mathbf{v}}_{{\rm A}}^{\rm ext}(t)$ are given by [14]

[TABLE]

and we input them to module B. The algorithm is closed and the algorithm iterates until convergence or the maximum number of iterations is reached. The MVALSE-EP algorithm is summarized as Algorithm 1. For further details, please refer to [16].

IV-C Computation Complexity

Let ${\rm Iter}_{\rm max}$ denote the number of iterations. For the proposed MVALSE-EP algorithm shown in Algorithm 1, it consists of main steps: The initialization, the componentwise MMSE operation, the mulisnapshot VALSE (MVALSE) algorithm. For initialization, the computation complexity is $O(N^{2}T)$ to obtain the posterior PDFs of the frequencies. The complexity of the MMSE module per iteration is $O(NT)$. As for the MVALSE algorithm, its complexity per iteration is dominated by calculating $\Delta_{k}$ (104), whose complexity is $O(N\hat{K}^{3}T)$ per iteration. Thus the overall computation complexity of the MVALSE-EP is $O((N^{2}+N\hat{K}^{3}\times{\rm Iter}_{\rm max})T)$.

IV-D Further Discussion

It is worth noting that the $k$th target is active only when (104) is satisfied. This is different from the MVALSE algorithm [12] where the $k$th target is active only when

[TABLE]

Because $\ln\frac{\rho}{1-\rho}<0$ as $\rho<0.5$ in general, it can be seen that with all the parameters being the same, the $k$th target is more likely to be active for the previous MVALSE algorithm [12]. We have found that if (115) is adopted as the criterion for activating the frequency, it is more likely to generate false alarms corresponding to harmonic components. Thus, the proposed criterion (104) reduces the false alarms.

It is also numerically found that MVALSE-EP tends to fit some spurious components and overestimates the number of sources for real data. As a consequence, we make the active criterion more harsh when we activate the $k$th frequency by minusing a positive bias $\lambda$ as

[TABLE]

where $\lambda>0$ is a tune parameter. For example, $\lambda=5$. Note that increasing $\lambda$ will make the model more sparse.

V Numerical Simulation

In this section, numerical experiments are conducted to evaluate the effectiveness of the proposed algorithm, through comparing with the CRBs. For the numerical simulation part, we conduct the DOA estimation and the element spacing of the linear array is half wavelength. Thus the DOA ${\bm{\theta}}$ is related to the frequency $\bm{\omega}$ via ${\bm{\omega}}=\pi\sin\left({\bm{\theta}}\right)$. For the real data case, we evaluate the range estimation for mmWave FMCW system. We evaluate the signal ${\mathbf{Z}}$ estimation error, the DOA estimation error, the correct model order estimation probability under quantized measurements in numerical simulations.

The phases of the weight coefficients are drawn i.i.d. from a uniform distribution between $[0,2\pi)$. For multi-bit quantization, a uniform quantizer is adopted and the quantization interval is restricted to $[-3\sigma_{z},3\sigma_{z}]$, where $\sigma_{z}^{2}$ is the variance of the signal $\Re\{{\mathbf{z}}_{N}\}$ or $\Im\{{\mathbf{z}}_{N}\}$. In our setting, it can be calculated that $\sigma_{z}^{2}\approx\sum\limits_{k=1}^{K}g_{k}^{2}+\sigma^{2}$. For one-bit quantization, zero is chosen as the threshold. For each target, its SNR (dB) is uniformly drawn from $[{\rm SNR}_{\rm min},{\rm SNR}_{\rm min}+\Delta]$ (dB). The amplitude of the target is determined from the SNRs and is fixed across the snapshots, the phase is uniformly drawn from $[0,2\pi)$. The noninformative prior, i.e., $p(\theta_{i})=1/(2\pi)$ is adopted for the MVALSE and MVALSE-EP algorithms. The number of maximum iterations is set as ${\rm Iter}_{\rm max}=50$. The normalized MSE (NMSE) of signal $\hat{\mathbf{Z}}$ (for unquantized and multi-bit quantized system) and the root MSE (RMSE) of the DOAs $\widehat{\bm{\theta}}$ are defined as ${{\rm NMSE}}(\hat{\mathbf{Z}})\triangleq\|\hat{\mathbf{Z}}-{\mathbf{Z}}\|_{\rm F}^{2}/\|{\mathbf{Z}}\|_{\rm F}^{2}$ and ${{\rm RMSE}}(\bm{\theta})\triangleq\sqrt{\|\widehat{\bm{\theta}}-\bm{\theta}\|_{2}^{2}/K}$, respectively. Please note that, due to magnitude ambiguity, it is impossible to recover the exact magnitude of ${x}_{k}(t)$ from one-bit measurements in the noiseless scenario. Thus for one-bit quantization, the debiased NMSE is used. The debiased NMSE of the signal defined as ${{\rm dNMSE}}({\mathbf{Z}})\triangleq\underset{\mathbf{c}}{\operatorname{min}}~{}{\|{\mathbf{Z}}^{*}-{\rm diag}({\mathbf{c}})\hat{\mathbf{Z}}\|_{\rm F}^{2}}/{\|{\mathbf{Z}}^{*}\|_{\rm F}^{2}}$ are calculated. As for the DOA error, we average only the trials in which all those algorithms estimate the correct model order. All the results are averaged over $300$ Monte Carlo (MC) trials unless stated otherwise. The empirical probability of correct model order estimation ${\rm P}(\hat{K}=K)$ is adopted as a performance metric. The RMSE of the DOA is calculated only when the model order is correctly estimated, i.e., $\hat{K}=K$.

V-A Asymptotic Property of CRB for a Single Source

This experiment is to validate the asymptotic property (83). Fig. 3 shows that the CRB under one bit quantization is inversely proportional to the cubic of the number of measurements $N$, the number of snapshots. For low SNR, the CRB under one bit quantization is inversely proportional to the SNR. While for high SNR, the CRB under one bit quantization is inversely proportional to the square root of the SNR. It can be seen that the asymptotic bound (83) provides a very good approximation for the true CRB. For the high SNR scenario shown in Fig. 3, the approximation is not accurate when $T=1$, because our bound is obtained by averaging over the phase of the amplitude of the frequency. As $T$ increases to $10$, (83) approximates the true CRB well.

V-B ${{\rm NMSE}}(\hat{\mathbf{Z}})$* versus the Number of Iteration*

The first experiment evaluates the ${{\rm NMSE}}(\hat{\mathbf{Z}})$ versus the number of iterations for $B=1,3,5$ and unquantized measurements. The DOAs are fixed as $[-3^{\circ},2^{\circ},75^{\circ}]$, and the SNRs are $[8,16,12]$ dB, respectively. The results are averaged over $50$ MC trials and are presented in Fig. 4. Note that MVALSE-EP converges very quickly. Meanwhile, the NMSE performance of the MVALSE-EP improves as the bit-depth increases.

V-C Performance versus SNR

The performance of the MVALSE-EP versus SNR is investigated and the results are plotted in Fig. 5. It can be seen that as ${\rm SNR}_{\rm min}$ increases, the NMSE of the signal decreases and slows down, especially for $1$ bit quantization. The model order estimation probability increases as ${\rm SNR}_{\rm min}$ increases and approaches $1$ for ${\rm SNR}_{\rm min}\geq-4$ dB. For $3$ bit quantization and no quantization, the RMSE of the DOA of the MVALSE-EP approaches to the CRB for ${\rm SNR}_{\rm min}\geq-4$ dB. While for $1$ bit quantization, there always exists a performance gap.

V-D Performance versus Number of Snapshots

The performance of the MVALSE-EP versus the number of snapshots $T$ is investigated and the results are plotted in Fig. 6. It can be seen that the NMSE of the signal decreases as the number of snapshots $T$ increases and saturates for $T\geq 50$. The successful model order estimation probability increases as the number of snapshots increases and approaches $1$. As for the RMSE of DOAs, it can be seen that the RMSE is significantly larger than their corresponding CRBs under $1$ bit, $3$ quantization and no quantization for $T\leq 2$. The reason is that when the number of snapshots is small, the number of the estimated DOAs is $3$ while the DOAs are estimated with large error, i.e., some DOAs are not detected while some false alarms are generated. As the number of snapshots increases, the probability of such events decreases and the RMSEs of DOAs approaches to the CRB. This demonstrates that increasing the number of snapshots benefits the DOA estimation.

V-E Real Data

Thus subsection uses the LFMCW AWR1642 Single-Chip to perform the range estimation. The chirp rate is set as $\kappa=29.982\times 10^{12}{\rm M}/{\rm Hz}^{2}$, the sampling frequency is $F_{s}=10$ MHz, and the maximum distance is $r_{\rm max}=\frac{cF_{s}}{2\kappa}=50$ m where $c=3\times 10^{8}$ denotes the speed of the electromagnetic wave. $\lambda$ (116) is set as 6 and 15 for 1 bit quantization and 12 bit quantization, respectively. Given that MVALSE-EP outputs the frequencies estimates $\hat{\bm{\omega}}$, the range estimates are $\hat{\mathbf{r}}=\frac{\hat{\bm{\omega}}}{2\pi}r_{\rm max}$. The number of fast time samples is $N=128$ and the number of snapshots is $16$. We first put corners as targets and test the algorithm in three experiments, as shown in Fig. 7. For the first experiment in Fig. 7, corner 1 is set just in front of the radar with the measured radial distance being 2.25m. For experiment 2, a little larger corner named corner 2 is added with the measured horizontal and vertical distance being $0.25$m and $2.5$m, and the radial distance can be calculated as $\sqrt{0.25^{2}+2.5^{2}}\approx 2.51$m. For experiment 3, the position of corner 1 is kept unchanged, and the measured horizontal and vertical distance of corner 2 is $1$m and $3.75$m, corresponding to the radial distance $3.88$m. A much larger Corner named corner 3 is added with the measured horizontal and vertical distance being $1$m and $5.25$m, and the radial distance can be calculated as $\sqrt{0.25^{2}+2.5^{2}}\approx 5.34$m.

For experiment 1, the normalized spectrum and the results are shown in Fig. 8 and Table I, respectively. It can be seen that for both 1 bit and 12 bit quantization, the highest peak occurs at $2.34$m. For 1 bit quantization, the second peak occurs at $42.96\approx 50-2.34\times 3=42.98$m, corresponding to the third harmonic and is $10$ dB lower than the highest peak. For 1 bit quantization, MVALSE-EP estimates the range of corner 1 as $2.32$ m for the single snapshot and multiple snapshots. While for 12 bit quantization, the number of detected targets is $3$, and MVALSE-EP still estimates the range of corner 1 as $2.32$ m for the single snapshot and multiple snapshots.

For experiment 2, Fig. 9 shows that FFT approach can not resolve the two targets. Table II show that MVALSE-EP detects the two corners and estimates the ranges. Besides, the amplitude of the first corner is $46$ dB, consistent with the results obtained in scenario 1. Also, the amplitude of corner 1 is about 4 dB lower than the second corner.

For experiment 3, the spectrum in Fig. 10 shows that after 1 bit quantization, the peak corresponding to corner 3 disappears. For the reconstruction results shown in Table III, MVALSE-EP only detects the strongest target corresponding to corner 1 with a single snapshot. For multiple snapshots, corner 2 is also detected and its amplitude is $6.68$ dB lower than corner 1. For high resolution data, MVALSE detects the three corners, and also the leakage component whose amplitude is comparable to corner $3$. This demonstrates that under 1 bit quantization, algorithms can not detect the weak target in the presence of the strong target, where the weak target (about $23$ dB) is $25$ dB lower than the strongest target (about $48.6$ dB).

For the fourth experiment shown in Fig. 11, the ranges of people 1 and people 2 are 3.42 m and 2.27 m, respectively. The reconstruction results are shown in Table IV. It can be seen that even under 1 bit quantization, MVALSE-EP detects two people and outputs the range estimates as $3.42$m and $2.35$m for single snapshot case, $3.42$m and $2.40$m for multiple snapshot case. Besides, the amplitude of people 1 and people 2 are about 27 dB and 21 dB, respectively.

For the last experiment, the ranges of the bicycle and people 1 are 3.41m and 2.77m, respectively. Results are shown in Table V. It can be seen that the bicycle and people 1 are detected under 1 bit and 12 bit quantization. Besides, the amplitudes of the bicycle and people 1 are about 31 dB and 28 dB, respectively. Note that the amplitudes of people 1 is consistent with the results in experiment 4.

VI Conclusion

In this paper, the multisnapshot line spectral estimation problem from coarsely quantized measurements is studied. The CRB is derived and the effects of the system parameters such as SNR, number of measurements and snapshots on the estimation accuracy are revealed in a single frequency scenario. Then, an MVALSE-EP algorithm which automatically estimates the frequencies, the model order and noise variance is proposed. Substantial numerical experiments including real data set are conducted to show the effectiveness of the MVALSE-EP.

VII Appendix

VII-A CRB for the general case

Note that

[TABLE]

where $\bar{\mathbf{A}}$ and $\bar{\mathbf{D}}(t)$ are defined in (31c) and (31j), respectively.

Substituting (117) in (7) yields

[TABLE]

where $\bar{\mathbf{H}}(t)$ and $\bar{\bm{\Delta}}(t)$ are defined in (31p) and (31q), respectively, $\bar{\bm{\Gamma}}$ is

[TABLE]

Note that for the inverse of a block matrix, one has

[TABLE]

where ${\mathbf{S}}={\mathbf{C}}-{\mathbf{B}}^{\rm T}{\mathbf{A}}^{-1}{\mathbf{B}}$. Therefore, ${\rm CRB}({\bm{\omega}})$ is

[TABLE]

and performing further simplification yields (32).