A Fast Power Spectrum Sensing Solution for Generalized Coprime Sampling

Abstract

With the growing scarcity of spectrum resources, wideband spectrum sensing is necessary to process a large volume of data at a high sampling rate. For some applications, only second-order statistics are required for spectrum estimation. In this case, a fast power spectrum sensing solution is proposed based on the generalized coprime sampling. The solution involves the inherent structure of the sensing vector to reconstruct the autocorrelation sequence of inputs from sub-Nyquist samples, which requires only parallel Fourier transform and simple multiplication operations. Thus, it takes less time than the state-of-the-art methods while maintaining the same performance, and it achieves higher performance than the existing methods within the same execution time without the need to pre-estimate the number of inputs. Furthermore, the influence of the model mismatch has only a minor impact on the estimation performance, allowing for more efficient use of the spectrum resource in a distributed swarm scenario. Simulation results demonstrate the low complexity in sampling and computation, thus making it a more practical solution for real-time and distributed wideband spectrum sensing applications.

1. Introduction

The demand for spectrum resources is increasing due to the rapid development of low-orbit satellite constellation systems (e.g., SpaceX, OneWeb), 5G/6G networks, and the Internet of Things (IoT) [1,2]. These applications are driving an unprecedented increase in demand for wideband spectrum sensing. Correspondingly, direct sampling requires a high-speed analog-to-digital converter (ADC) [3] based on the Shannon–Nyquist sampling theorem, increasing data volume and energy consumption.

Currently, the most widely used methods are sweep-tune sampling and filter band sampling. Both of these methods fall under the category of low-speed sampling. However, the scanning scheme has a detection latency and may miss short-lived signals [4]. In addition, the filter band scheme has a complicated structure and is prone to serious channel crosstalk [5]. Consequently, there has been a trend towards using wideband spectrum sensing as a guide.

The recent compressive sensing (CS) theory provides a sub-sampling scheme that offers low-speed and large instantaneous bandwidth by utilizing the sparsity in the frequency domain [6,7]. The typical CS schemes include the analog-to-information converter (AIC), multi-coset sampling (MCS), and multi-rate sampling (MRS). The typical AIC architecture is the modulated wideband converter (MWC), which presents significant challenges for the implementation of the Nyquist-rate pre-randomizing [8]. The MCS architecture is currently a preferred scheme for the ADC with high-speed and high-precision, which cannot obtain a high significant bit because of the mismatching between multiple channels [9]. The MRS architecture is currently a preferred scheme for the sparse array signal acquisition, whose performance is limited by time synchronization accuracy [10]. Obviously, the engineering implementation of these schemes presents major difficulty, which restricts their applications. Furthermore, the CS-based methods [3,11] for sparse signal recovery have high computational complexity and are extremely sensitive to model matching and noise, which mainly include greedy methods, convex methods, and others. Therefore, reducing the dependence of signal processing algorithms on model mismatches between the digital and the analog world is an interesting research idea.

To overcome these difficulties, the approach to spectrum estimation has shifted from processing the original signal to analyzing its second-order statistics. In this context, the compressive covariance sensing (CCS) theory provides a wideband spectrum sensing scheme that operates at a low-speed and has a large instantaneous bandwidth. It is reliable even in environments with a low signal-to-noise ratio (SNR) and in non-sparse conditions. Generally, according to the methods of computation, the CCS-based wideband power spectrum sensing can be divided into the time-domain approach [12] and the frequency-domain approach [13]. The time-domain approach establishes a relationship between the original inputs and the output samples through the selection matrix with zero and one elements under the equivalent Nyquist-rate sampling. As well as this, the frequency-domain approach builds the relationship between both frequency representations of them. However, these studies mainly focus on the MWC and MCS schemes.

For the MRS scheme, the existing CCS-based methods include entropy function minimization [14], matrix norm minimization [15], and Toeplitz matrix completion [16], among others. However, all these methods are based on the reconstruction of the covariance matrix and the use of multiple signal classification (MUSIC) algorithm, which have high computational complexity and require pre-estimation of the number of signals. Additionally, they are also sensitive to model matching.

Furthermore, a computationally efficient method is developed, which is based on the relationship between the autocorrelation sequence and sub-sampling samples of the MCS scheme [17]. Building on this, a fast solution for generalized coprime sampling is introduced, which utilizes only parallel FFT and multiplication operations. As a result, it achieves a reduced time and low estimation error, presenting a trade-off between system performance and the number of degrees of freedom (DOFs). Moreover, model mismatch has minimal effect on performance, making a more practical solution for real-time and distributed wideband spectrum sensing applications.

The rest of this paper is organized as follows: Section 2 describes the signal model and the proposed fast power spectrum sensing solution. Section 3 conducts an analysis and validation of the proposed solution through simulation. The discussion is presented in Section 4.

Notations: The bold characters denote vectors. The notations $\mathbb{R}$, $\mathbb{N}$, and ${\mathbb{N}}^{+}$ represent the set of real numbers, nonnegative integers, and positive integers, respectively. The superscripts ${(·)}^T$ and ${(·)}^H$ indicate the transpose and conjugate transpose of a vector or a matrix, respectively. The operator ∘ signifies the Hadamard product, ${|·|}^2$ signifies the element-wise square modulus of a vector, and ceil(·) signifies round up to an integer. The symbols ${\mathit{F}}_a$ and ${\mathit{F}}_a^{-1}$ mean the a-point fast Fourier transofrm (FFT) and the inverse fast Fourier transofrm (IFFT), respectively.

2. Materials and Methods

2.1. Signal Model

The generalized coprime sampling architecture comprises two uniform sub-Nyquist sampling channels, whose sampling periods are coprime multiples of the Nyquist sampling period. The introduction of two additional operations, the multiple coprime unit factor $p\in {\mathbb{N}}^{+}$ and the non-overlapping factor $q\in {\mathbb{N}}^{+}$, enhances the number of DOFs and improves the estimation accuracy. Consequently, the coprime sampling scheme is presented with sampling intervals $r_0{T}_s$ and $r_1{T}_s$, as depicted in Figure 1. Without loss of generality, it is assumed that $r_0<r_1$ with $r_0,r_1\in {\mathbb{N}}^{+}$, where $r_0$ and $r_1$ are coprime. Additionally, the sampling interval $T_s$ corresponds to the Nyquist sampling rate $f_s$.

For a wide-sense stationary or cyclostationary process $X\left(t\right),t\in \mathbb{R}$, it consists of a number of frequencies that are confined within the bandwidth $B_s\le f_s/2$. The outputs of two uniform sub-Nyquist sampling channels can be expressed as

$$\begin{array}{cc}\hfill y_0\left[n_0\right]& =x\left[r_0{n}_0\right]=X\left(r_0{n}_0{T}_s\right),n_0\in \mathbb{N}\hfill \\ \hfill y_1\left[n_1\right]& =x\left[r_1{n}_1\right]=X\left(r_1{n}_1{T}_s\right),n_1\in \mathbb{N}\hfill \end{array}$$

(1)

where $x\left[n\right],n\in \mathbb{N}$ denotes the Nyquist sampling samples, and the highest sampling rate of the coprime sampling system is given by $1/\left(r_0{T}_s\right)=f_s/r_0$.

Accordingly, the elements of the sensing vector corresponding to the two coprime samplers can be denoted as

$$a_0\left[i\right]=\left\{\begin{array}{cc}1,\hfill & i=r_0{l}_0+k{r}_0{r}_1\hfill \\ 0,\hfill & \mathrm{elsewhere}\hfill \end{array}\right.$$

(2)

and

$$a_1\left[i\right]=\left\{\begin{array}{cc}1,\hfill & i=r_1{l}_1+(k+q)r_0{r}_1\hfill \\ 0,\hfill & \mathrm{elsewhere}\hfill \end{array}\right.$$

(3)

where $l_0=0,1,\dots ,r_1-1$, $l_1=0,1,\dots ,r_0-1$, and $k=0,1,\dots ,p-1$.

From a data acquisition perspective, the output samples obtained from the generalized coprime sampling scheme are a subset of the Nyquist samples, positioned at

$$\mathbb{P}=\{r_0{l}_0+k{r}_0{r}_1\}\cup \{r_1{l}_1+(k+q)r_0{r}_1\}$$

(4)

Based on the sensing vectors, the relationship between the elements of the generalized coprime sampling vector and the Nyquist sampling vector can be expressed as

$$y\left[n\right]=\left\{\begin{array}{cc}x\left[n\right],\hfill & n\in \mathbb{P}\hfill \\ 0,\hfill & \mathrm{elsewhere}\hfill \end{array}\right.$$

(5)

and the elements of the associated sensing vector are defined as

$$a\left[n\right]=\left\{\begin{array}{cc}1,\hfill & n\in \mathbb{P}\hfill \\ 0,\hfill & \mathrm{elsewhere}\hfill \end{array}\right.$$

(6)

As a result, there is

$$\mathbf{y}\left[n\right]=\mathbf{a}\left[n\right]\circ \mathbf{x}\left[n\right]$$

(7)

where $\mathbf{x}\left[n\right]={\left[x\left[0\right],x\left[1\right],\dots ,x[N-1]\right]}^T$, $\mathbf{a}\left[n\right]=\left[a\right[0],a[1],$ $\dots ,a[N-1]{]}^T$, and $\mathbf{y}\left[n\right]={\left[y\left[0\right],y\left[1\right],\dots ,y[N-1]\right]}^T$ are all vectors of size $N\times 1$.

2.2. Proposed Fast Solution

Considering the widely-used unbiased estimation of the autocorrelation sequence for the output of generalized coprime sampling, the elements $r_y\left[m\right]$ can be expressed as

$$r_y\left[m\right]=\frac{1}{N}·\mathbf{y}\left[n\right]{\mathbf{y}}^H[n-m]$$

(8)

Substituting Equation (7) into Equation (8) results in

$$\begin{array}{cc}\hfill r_y\left[m\right]& =\frac{1}{N}·(\mathbf{a}\left[n\right]\circ \mathbf{x}\left[n\right])({\mathbf{a}}^H[n-m]\circ {\mathbf{x}}^H[n-m])\hfill \\ \hfill & =\frac{1}{N}·\left(\left(\mathbf{a}\left[n\right]{\mathbf{a}}^H[n-m]\right)\circ \left(\mathbf{x}\left[n\right]{\mathbf{x}}^H[n-m]\right)\right)\hfill \\ \hfill & =r_a\left[m\right]\circ r_x\left[m\right]\hfill \end{array}$$

(9)

where $\left|m\right|\le N-1$. Thus, the power spectrum can be obtained by performing FFT on the autocorrelation sequence $\left\{r_x\left[m\right]\right\}$, which is derived from the autocorrelation sequence $\left\{r_y\left[m\right]\right\}$ and $\left\{r_a\left[m\right]\right\}$. Therefore, a computationally efficient practical solution is employed to obtain the estimation of the autocorrelation sequences. The steps are as follows:

Step 1: Pad vectors $\mathbf{a}\left[n\right]$ and $\mathbf{y}\left[n\right]$ with an additional N zeros.

$$a_{2N}\left[n\right]=\left\{\begin{array}{cc}{a}_N\left[n\right],\hfill & 0\le n\le N-1\hfill \\ 0,\hfill & N\le n\le 2N-1\hfill \end{array}\right.$$

(10)

and

$$y_{2N}\left[n\right]=\left\{\begin{array}{cc}{y}_N\left[n\right],\hfill & 0\le n\le N-1\hfill \\ 0,\hfill & N\le n\le 2N-1\hfill \end{array}\right.$$

(11)

Step 2: Calculate the autocorrelation sequence based on the power spectrum estimation of vector ${\mathbf{a}}_{2N}\left[n\right]=[a\left[0\right],a\left[1\right],$ $\dots ,a[2N-1]{]}^T$ and ${\mathbf{y}}_{2N}\left[n\right]={\left[y\left[0\right],y\left[1\right],\dots ,y[2N-1]\right]}^T$ by involving FFT and IFFT.

$${\widehat{\mathbf{r}}}_a^{\prime }\left[k\right]={\mathit{F}}_{2N}^{-1}{\left|{\mathit{F}}_{2N}{\mathbf{a}}_{2N}\right|}^2/N$$

(12)

and

$${\widehat{\mathbf{r}}}_y^{\prime }\left[k\right]={\mathit{F}}_{2N}^{-1}{\left|{\mathit{F}}_{2N}{\mathbf{y}}_{2N}\right|}^2/N$$

(13)

where $k=0,1,\dots ,2N-1$.

Step 3: Truncate the autocorrelation sequence of interest according to the frequency resolution of the system $\mathsf{\Delta }f$.

$${\widehat{r}}_y\left[m\right]=\left\{\begin{array}{cc}{\widehat{r}}_y^{\prime }\left[m\right],\hfill & 0\le m\le M-1\hfill \\ {\widehat{r}}_y^{\prime }[m+2N],\hfill & -M+1\le m\le -1\hfill \end{array}\right.$$

(14)

and

$${\widehat{r}}_a\left[m\right]=\left\{\begin{array}{cc}{\widehat{r}}_a^{\prime }\left[m\right],\hfill & 0\le m\le M-1\hfill \\ {\widehat{r}}_a^{\prime }[m+2N],\hfill & -M+1\le m\le -1\hfill \end{array}\right.$$

(15)

where $M=\mathit{ceil}(fs/2/\mathsf{\Delta }f)+1$.

Step 4: Compute the autocorrelation sequence of the inputs using the obtained sequences $\left\{{\widehat{r}}_y\left[m\right]\right\}$ and $\left\{{\widehat{r}}_a\left[m\right]\right\}$.

$${\widehat{r}}_x\left[m\right]={\widehat{r}}_y\left[m\right]/{\widehat{r}}_a\left[m\right]$$

(16)

where $m=-M+1,\dots ,-1,0,1,\dots ,M-1$.

Step 5: Obtain the power spectrum estimation by taking the FFT of the vector ${\widehat{\mathbf{r}}}_x\left[m\right]={\left[{\widehat{r}}_x[-M+1],\dots ,{\widehat{r}}_x\left[1\right],{\widehat{r}}_x\left[0\right],{\widehat{r}}_x\left[1\right],\dots ,{\widehat{r}}_x[M-1]\right]}^T$.

$${\widehat{\mathbf{S}}}_x\left(\omega \right)=\left|{\mathit{F}}_{2M-1}{\widehat{\mathbf{r}}}_x\left[m\right]\right|$$

(17)

The block diagram that illustrates the fast power spectrum sensing solution for generalized coprime sampling is depicted in Figure 2. The proposed solution is efficient because it only involves FFT/IFFT operations and some basic multiplication operations. Moreover, as shown in the red section of Figure 2, the autocorrelation sequence $\left\{{\widehat{r}}_a\left[m\right]\right\}$ of the sensing vector can be pre-calculated offline. This calculation solely depends on the generalized coprime sampling scheme.

Consequently, the computational complexity of the proposed solution involves performing the FFT on a $\left(2N\right)$-point sequence twice, resulting in $\left(2N\right)\log \left(2N\right)$ floating-point operations according to (13). In addition, the FFT is performed on a $(2M-1)$-point sequence once, leading to $(2M-1)\log (2M-1)$ floating-point operations according to (17). After incorporating $2M-1$ multiplication calculations, the total computational complexity requires $\left(4N\right)\log \left(2N\right)+(2M-1)\log (2M-1)+(2M-1)$ floating-point operations. This leads to a lower computational complexity compared to the state-of-the-art methods. Meanwhile, it is feasible to effectively compute the FFT operations in parallel, making it a more suitable practical solution for real-time wideband power spectrum sensing applications.

As researched in the state-of-art, there is an example of wideband spectrum sensing based on the MWC architecture with a 1 GHz bandwidth, which requires a spectrum resolution of 10 kHz. The time-domain approach involves at least ${10}^{14}$ floating-point operations in total, assuming that the number of sampling branches is set to 8, the downsampling factor sets to 25, then the number of output samples need to be 4000, and the number of samples used to calculate the correlation matrix is set to 100. Moreover, the more efficient time-domain approach has the same computational complexity as the frequency-domain approach, which involves more than ${10}^7$ floating-point operations in total under the same assumption. Additionally, regarding the MCS schemes, $1.08\times {10}^7$ floating-point operations are needed. In comparison, the proposed method involves $8.805\times {10}^6$ floating-point operations in total.

3. Results

In the experiments, it is assumed that there are I inputs with identical powers, which are distributed in the frequency band $[2,18]$ GHz. Subsequently, the coprime integers $r_0=3$, $r_1=4$ and the Nyquist sampling rate $f_s=32$ GHz are set. Furthermore, the relative root mean square error (RMSE) is adopted to evaluate the performance of the proposed fast power spectrum sensing method, which is defined as follows:

$$\mathrm{Relative}\mathrm{RMSE}\left(f_i\right)=\frac{1}{f_s}\sqrt{\frac{1}{500I}\sum _{j=1}^{500}\sum _{i=1}^I{({\widehat{f}}_i\left(j\right)-f_i)}^2}$$

(18)

where ${\widehat{f}}_i\left(j\right)$ is the estimation of $f_i$ from the jth Monte Carlo trial, and five hundred Monte Carlo trials are conducted.

3.1. Estimated Power Spectrum Performance

Herein, the estimated power spectrum results are initially displayed, with $p=3000$ and an input SNR of 15 dB. As shown in Figure 3, there are $I=50$ mono-frequency pulse (MP) signals, which are randomly distributed. Figure 4 depicts $I=20$ binary phase shift keying (BPSK) signals for 1 M symbols per second with random frequency and code. Figure 5 presents $I=2$ linear frequency modulation (LFM) signals with 10 GHz bandwidth under $\pm 6$ GHz initial frequencies. Figure 6 shows $I=21$ mixture signals of three types. As can be observed, all frequencies are estimated accurately with the proposed method.

Furthermore, there are multiple LFM signals that share the same carrier frequency but have different quadratic modulation coefficients. The frequency and bandwidth of these signals are randomly selected. Compared with Figure 7, the power spectrum estimation shows an increase in the number of pseudo-spectra as more signals are aliased, and as the number of aliased signals increases, as shown in Figure 8 and Figure 9, respectively.

3.2. Relative RMSE Performance

As depicted in Figure 10, the RMSE results are compared as a function of the input SNR, where $I=18$ MP signals are utilized and the frequency is randomly selected. It is observed that the RMSE tends to stabilize when SNR is greater than −2 dB for $p=300$. Under the same conditions, the original method exhibits better performance at a low SNR by utilizing the Toeplitz matrix completion. Meanwhile, the estimation performance improves as p increases, due to the fact that the DOF increases with p, leading to improved resolution. Clearly, the performance of the MCS scheme is the worst with the same number of samples as $p=3000$, which is due to the probability of signal loss during the process of power spectrum estimation. As expected, Figure 11 presents the same result. However, the selection of coprime sampling rate makes less of a difference to performance when p is greater than 1000. This is because the system redundancy under simulation is sufficient.

3.3. Execution Time Performance

Furthermore, the multiple coprime unit factor p not only affects the resolution, but also determines the execution time of algorithms. Consequently, the execution time results, as a function of p, are compared as illustrated in Figure 12, where $I=10$ MP signals are utilized and the frequency is randomly selected at 0dB SNR. The computing environment is based on Windows 11, equipped with an AMD Ryzen 5 3500U processor, Radeon Vega Mobile Gfx at 2.10 GHz, and 20.0 GB of RAM from Lenovo in Beijing, China. It is evident that the execution time of the proposed method has significant advantages over the matrix completion method. Analogous to Figure 13, the proposed method surpasses the matrix completion method under identical execution time. Additionally, the proposed method requires less time than the matrix completion method under the same performance. However, both methods depend on a larger number of samples. Therefore, a tradeoff exists between execution time and system performance.

3.4. Model Mismatching Performance

Finally, the influence of the model matching degree between the sensing vector and measurements on the performance is discussed in Figure 14 with the 5 dB SNR. Here, $I=18$ frequencies are randomly selected for MP signals, and different time delays are used, which are unknown to the sensing vector. As a result, the influence of the model mismatch causes minimal fluctuation in RMSE within 200 $\mathsf{\mu }$s. This is interesting, as it suggests that the fact potentially enables more efficient utilization of the spectrum resource in a distributed swarm scenario.

4. Discussion

A fast power spectrum sensing solution for generalized coprime sampling is proposed that only uses the parallel FFT and simple multiplication operations. It has obvious advantages over existing methods in terms of spectrum estimation performance and execution time. Moreover, there is no need to pre-estimate the number of inputs. Furthermore, the influence of the model mismatch has minimal impact on the spectrum estimation performance. That makes it more suitable for further discussion on its application in the distributed swarm scenario.