The problem of multidimensional (MD) frequency estimation arises in many applications. For the -D frequency estimation in wireless communication systems, the two dimensions correspond to the delay and angle of each path. Physical arguments and a growing body of experimental evidence suggest that utilizing the structures benefits and improves the channel reconstruction performance . Another application is for the -D frequency estimation in radar systems, one can interpret one dimension as time delays, one dimension as Doppler shifts, and the other as spatial frequencies. Since each pair of time delay, Doppler shift and spatial frequency specifies the localization, speed and direction of arrival (DOA) of a scatter, respectively, estimating these parameters is of great importance for target localization and tracking . Traditional methods including -D Discrete Fourier Transform (DFT), subspace based approaches such as D MUSIC , D ESPRIT , Matrix Enhancement Matrix Pencil (MEMP)  and the multidimensional folding (MDF) techniques [6, 7].
Recently, compressed sensing based approaches have been proposed for -D frequency estimation, especially for -D frequency estimation. Provided that the frequencies lie on the DFT grid, the signal can be recovered via efficient minimization . However, due to spectral leakage of the off-grid frequencies, the signal is not exactly sparse and conventional compressed sensing (CS) algorithms may result in significant performance degradation . Consequently, grid refinement and continuous CS approaches have been developed. In [10, 11], the Newtonalized orthogonal matching pursuit (NOMP) is proposed for -D frequency estimation, which first detects the grid and then iteratively refines the frequencies. Besides, the model order is estimated by the constant false alarm (CFAR) criterion with the knowledge of the noise variance. The computation complexity of NOMP is low, and the estimation performance is good. For the latter,  proposes a nuclear norm minimization over multi-fold Hankel matrices approach,  develops a reweighted trace minimization (RWTM) plus matrix pencil and (auto-)pairing (MaPP) method, where RWTM fully exploits the multilevel Toeplitz (MLT) structure and utilizes the sparsity to obtain the covariance estimates, and MaPP calculates the Vandermonde decomposition. As [12, 13] involve solving a semidefinite programming (SDP), their computation complexities are very high.
A more recent approach is the variational line spectral estimation (VALSE) algorithm, which uses a complete Bayesian treatment by imposing the frequencies as random variables. Such method naturally allows for representing and operating with the uncertainty of the frequency estimates, in addition to only that of the weights. Numerically, accounting for the frequency uncertainty benefits the line spectral estimation performance. VALSE automatically estimates the model order, noise variance, and provides the uncertain degrees of the frequency estimates. Besides, the computation is low, compared to the nuclear norm minimization based approaches. It was also extended to the multi-snapshot case in array signal processing  and quantized setting in high-speed sampling . In contrast, this work develops MDVALSE to deal with the multidimensional (MD) frequency estimation problem.
The main contribution of this work is as follows: Firstly, MDVALSE is rigorously developed, which is nontrival from the following two aspects: one is that the multifrequencies are treated and updated as a whole, another aspect is that even the probability density function (PDF) of the multifrequencies can be approximated as a multivariate von Mises distribution, computing expectations of is still hard. We novelly project the PDF of the multifrequencies as independent vonMises distribution, which allows computation of . Secondly, the -D fast Fourier transform (FFT) approach is adopted to reduce the computation complexity significantly. Finally, numerical experiments are conducted to demonstrate the excellent performance of MDVALSE.
Ii Problem Setup
Consider a measurement model for dimensional line spectra estimation problem expressed as
where is the noiseless line spectral composed of complex sinusoids
denotes the complex amplitude of the th frequency, and the th element is , is the th element of and , is the additive white Gaussian noise (AWGN) and , where
For simplicity, we define
Since the number of spectral is unknown, an overcomplete model is adopted, where the number of spectral is supposed to be and . Consequently, the binary hidden variables are introduced to characterize whether the th component is active or not. The probability mass function (PMF) is , where and
Besides, , where
follows a Bernoulli-Gaussian distribution
Let and be the model and estimated parameters. In general, one should first find the maximum likelihood (ML) estimates of by marginating the posterior PDF of , then find the maximum a posterior (MAP) or minimum mean squared error (MMSE) estimates of . However, the above approach is computationally intractable and thus an iterative algorithm is developed in the ensuing section.
Iii MDVALSE Algorithm
The variational approach tries to find a surrogate PDF via minimizing the Kullback-Leibler (KL) divergence between and , which is equivalent to maximize
Similar to , is imposed to have the following structures
Maximizing with respect to all the factors is also intractable. Similar to the Gauss-Seidel method, is optimized over each factor , and separately with the others being fixed. Let be the set of all latent variables. Maximizing (8) with respect to the posterior approximation of each latent variable yields
where the expectation is with respect to all the variables except and the constant ensures normalization of the PDF.
Due to the factorization property of (9), the frequencies can be estimated from the marginal distribution as
where returns the angle.
Given that , the posterior PDF of is
For the given posterior PDF , the mean and covariance of the weights is estimated as
Let be the set of indices of the non-zero components of , i.e., . In the following, we detail the procedures.
Iii-a Inferring the frequencies
where denotes the conjugate operation. While it is hard to obtain the analytical results of (11) for the PDF (14), is projected as the independent von Mises distribution. A Newton step is implemented to refine the previous estimates, i.e.,
where is the exponential part of and the detail of the inverse of is given in ,
returns a vector with elements being the diagonal of the matrix. As a result,and can be approximated as
Iii-B Inferring the weights and support
Next are fixed and is maximized w.r.t. . Define the matrices and as
where denotes the th element of . It can be shown that updating , and are the same as  and details are omitted here.
Iii-C Estimating the model parameters
Iii-D Summary of MDVALSE algorithm
The initialization of MDVALSE algorithm is presented. For the first step of initialization, we randomly extract one dimensional data from and use the method in  to initialize . For the later, we choose to initialize in a sequential manner. For the th step, the noncoherent PDF
is projected to independent von Mises distribution, where . Here N-D FFT is implemented to obtain the approximated independent von Mises distribution. Then the MDVALSE can be obtained and the procedures are similar to [14, Algorithm 3].
Now we analyze the computation complexity of MDVALSE. The complexity of initialization is dominated by dimensional FFT. For the th dimension, the number of grid is where can be viewed as the oversampling factor. Thus the initialization has complexity . For each iteration, the complexity is dominated by the estimate of and the projection of as independent VM. According to , the calculation of has complexity and in general . For the projection operation, the complexity is which are mainly dominated by the calculation and inverse of Hessian matrix in (16b). Therefore, the computational complexity of MDVALSE is where denotes the number of whole iterations. For , the computation complexity can be simplified as .
Iv Numerical Simulation
In this section, substantial numerical experiments are conducted to evaluate the performance of MDVALSE algorithm. The normalized mean squared error (NMSE) of , MSE of and the correct model order estimated probability are adopted as the performance metrics. The frequencies estimation error is averaged over the trials in which
for a given simulation point. We define nominal signal-to-noise ratio (SNR) as). The Algorithm stops when or , where is the number of iterations.
At first, the approximations of noncohent PDF (20) and coherent PDF (14) as independent von Mises distribution are validated numerically and results are shown in Fig. 1. Note that for , the approximation deviates a little away from the true PDF due to the inaccuracy of the FFT results. As increases to , the approximation is high and the two PDFs are indistinguishable around the peak.
Iv-a Estimation from 2D case
In this simulation, we consider a 2D case with signals and true frequencies are
The complex weight coefficients are generated i.i.d. from complex normal distribution. The number of measurements is . As for performance comparison, reweighted trace minimization (RWTM) algorithm proposed in  is chosen and the performances of RWTM and MDVALSE are investigated in Fig. 2. It can be seen that for noiseless case in Fig. 2, MDVALSE and RWTM algorithms both estimate well. Numerically, the MSE of frequencies of RWTM is lower than that of MDVALSE. For noisy case in Fig. 2, MDVALSE still performs well, while RWTM overestimates the model order. This demonstrates that MDVALSE is more robust against noise than RWTM.
Iv-B Estimation from 4D case
The performance of MDVALSE algorithm for 4D case versus SNR is investigated. The magnitudes and phases of the complex weight coefficients are generated i.i.d. from normal distribution, respectively. The number of sinusoids is and the wrap-around distance between any two frequencies is at least radians for each dimension. MC trails are performed and results are shown in Fig. 3. MDVALSE works well and performs better as SNR or number of measurements increases.
This paper develops the MDVALSE for multidimensional line spectral estimation, where the multidimensional frequencies are treated as a whole and their PDF are projected as independent von Mises distribution for tractable inference. MDVALSE inherits the advantages of automatically estimating the model order, number of nuisance parameters such as noise variance, and providing the uncertain degree of frequency estimates. Numerical results demonstrate the effectiveness of MDVALSE.
-  Y. Han, T-H. Hsu, C-K. Chao, and S. Jin, “Efficient downlink channel reconstruction for FDD multi-antenna systems,” IEEE Trans. Wireless Commun., vol. 18, no. 6, pp. 3161-3176, June 2019.
-  B. Jin, J. Zhu, Q. Wu, Y. Zhang and Z. Xu, “One-bit LFMCW radar: spectrum analysis and target detection,” to appear in IEEE Transactions on Aerospace and Electronic Systems.
-  Y. Hua, “A pencil-MUSIC algorithm for finding two-dimensional angles and polarizations using crossed dipoles,” IEEE Trans. Antennas Propag., vol. 41, no. 3, pp. 370-376, Mar. 1993.
-  J. Li and R. T. Compton, Jr., “Two-dimensional angle and polarization estimation using the ESPRIT algorithm,” IEEE Trans. Antennas Propag., vol. 40, no. 5, pp. 550-555, May 1992.
-  Y. Hua, “Estimating two-dimensional frequencies by matrix enhancement and matrix pencil,” IEEE Trans. Signal Process., vol. 40, no. 9, pp. 2267-2280, Sep. 1992.
-  X. Liu and N. D. Sidiropoulos, “Almost sure identifiability of constant modulus multidimensional harmonic retrieval,” IEEE Trans. Signal Process., vol. 50, no. 9, pp. 2366 C2368, Sep. 2002.
-  J. Liu, X. Liu, and X. Ma, “Multidimensional frequency estimation with finite snapshots in the presence of identical frequencies,” IEEE Trans. Signal Process., vol. 55, no. 11, pp. 5179-5194, Nov. 2007.
-  E. Candes and J. Romberg, “Sparsity and incoherence in compressive sampling,” Inverse problems, vol. 23, no. 3, p. 969, 2007.
-  Y. Chi, L. L. Scharf, A. Pezeshki, and R. Calderbank, “Sensitivity to basis mismatch in compressed sensing,” IEEE Trans. Signal Process., vol. 59, no. 5, pp. 2182-2195, 2011.
-  Z. Marzi, D. Ramasamy and U. Madhow, “Compressive channel estimation and tracking for large arrays in mm wave picocells,” IEEE J. Sel. Topics Signal Process., vol. 10, no. 3, pp. 514-527, 2016.
A. Gupta, U. Madhow and A. Arbabian, “Super-resolution in position and velocity estimation for short-range mm-wave radar,”50th Asilomar Conference on Signals, Systems and Computers, Nov. 2016, Pacifc Grove, USA.
-  Y. Chen and Y. Chi, “Robust spectral compressed sensing via structured matrix completion,” IEEE Trans. Inf. Theory, vol. 60, no. 10, pp. 6576-6601, Oct. 2014.
-  Z. Yang, L. Xie and P. Stoica, “Vandermonde decomposition of multilevel toeplitz matrices with application to multidimensional super-resolution,” IEEE Trans. Inf. Theory, vol. 62, no. 6, pp. 3685-3701, 2016.
-  M. A. Badiu, T. L. Hansen and B. H. Fleury, “Variational Bayesian inference of line spectra,” IEEE Trans. Signal Process., vol. 65, no. 9, pp. 2247-2261, 2017.
-  J. Zhu, Q. Zhang, P. Gerstoft, M. A. Badiu and Z. Xu, “Grid-less variational Bayesian line spectral estimation with multiple measurement vectors,” Signal Process., vol. 61, pp. 155-164, 2019.
-  J. Zhu, Q. Zhang and X. Meng, “Gridless variational Bayesian inference of line spectral estimation from quantized samples,” arXiv preprint, arXiv:1811.05680, 2019.
-  K. V. Mardia and P. E. Jupp, Directional Statistics. New York, NY, USA: Wiley, 2000.