An Effective Reconstruction Algorithm Based on Modulated Wideband Converter for Wideband Spectrum Sensing

In recent years, wideband spectrum sensing combined with sub-Nyquist sampling and compressed sensing technology in the field of cognitive radio has received widespread attention. However, the existing broadband detection methods based on sub-Nyquist sampling do not fully consider the spectrum feature changes caused by the sampling structure, resulting in unnecessary computational complexity in the support reconstruction process. In this paper, a novel reconstruction algorithm, called nearest orthogonal matching pursuit (N-OMP), is proposed based on modulated wideband converter (MWC) sub-Nyquist sampling structure. This algorithm utilizes the special power spectrum slicing characteristics caused by pseudo-random sequence and low-pass filtering. After an occupied subband is detected, it calculates the correlation coefficient between the residual vector and the column vectors corresponding to two adjacent subbands, based on which we can directly judge the occupancy of two adjacent subbands by comparing the size of the two correlation coefficients, thereby reducing the number of iterations of the reconstruction algorithm. Theoretical derivation and simulation experiment results show that, compared with the orthogonal matching pursuit (OMP) algorithm, the proposed algorithm can reduce the computational complexity by up to 50%, while showing better support reconstruction accuracy.


I. INTRODUCTION
With the rapid development of Internet of Things (IoT) technology and its applications, the demand for Internet access has grown rapidly. A large number of emerging wireless access services cannot obtain suitable certified frequency bands, and crowd in unlicensed bands. Under this trend, the static radio spectrum resource management method is no longer effective enough [1]- [3].
Cognitive radio (CR) adopts a dynamic spectrum management method, which can provide unauthorized users (secondary users) with dynamic spectrum access opportunities without causing communication interference to authorized users (primary users) [4]. It has the potential to solve the current shortage of spectrum resources and is a key technology to promote the further development of IoT and 5G [5], [6].
The associate editor coordinating the review of this manuscript and approving it for publication was Hayder Al-Hraishawi .
Spectrum sensing is an important part of cognitive radio. Traditional narrow-band spectrum sensing technology allows secondary users to judge the presence or absence of the primary user signal in a target channel, while broadband spectrum sensing requires analysis of multiple channels once. However, in the process of wideband spectrum sensing, high-speed sampling may increase the hardware burden and energy consumption [7], [8]. Meanwhile, high computational complexity and unsatisfactory perception performance caused by low signal to noise ratio (SNR) have not been effectively resolved [9], [10], which requires further research and development by researchers.
In order to realize wideband spectrum detection, Pei et al. added an adjustable narrow-band filter at the RF front end to detect the occupancy of the primary user segment by segment [11]. It is simple and easy to implement, but the detection period is long, which may cause missed detection and communication interference. Utilizing the sparseness of the signal, compressed sensing (CS) technology provides a new way to break through the above Nyquist sampling bottleneck [12]- [14], which greatly promotes the development of broadband spectrum detection. Combined with CS technology, the current mainstream and highly recognized sub-Nyquist sampling circuits are: (1) Analog-to-information converter (AIC) [15], [16]; (2) Modulated wideband converter (MWC) [17], [18]; (3) Multi-coset sampling (MCS) [19]. To provide the possibility of reconstructing the signal based on sub-Nyquist sampling, both MWC and AIC use a mixed operation with a pseudo-random sequence to superimpose the multi-band signal spectrum information to the baseband domain with different weights. MCS realizes it through staggered sampling phase. Based on the MCS structure, Yen et al. directly estimated the power spectrum of the wideband signal by the statistical characteristics of the broadband signal [20]. Furthermore, in [21], Liu et al. designed a single-channel multi-coset sampling structure based on a single-period nonuniform sampling clock, a single sample holder and a single ADC, and proposed a spectrum sensing method. Zhang et al. proposed a distributed cooperative sensing scheme that uses nearby IoT devices to perform multi-coset sampling and reconstruction to avoid the problem of unsatisfactory reconstruction accuracy of a single node in [22]. On the other hand, on the basis of MWC, Xu et al. in [23] used orthogonal projection as a band-stop filter in the compression domain to eliminate the primary users (PUs) signal components adaptively that have been determined in the reconstruction process, reducing the complexity of the algorithm and improving the reconstruction accuracy. In [24], Li et al. proposed a dynamic compressive wide-band spectrum sensing method based on channel energy reconstruction. It only needs to perform partial reconstruction according to the change of the channel occupancy state. In addition, in order to effectively deal with the wireless fading problem, Gong et al. proposed a multi-antenna generalized modulated converter to achieve sub-Nyquist sampling, and further proposed two compressed subspace algorithms to achieve broadband spectrum sensing in [25]. However, the existing wideband spectrum detection methods based on sub-Nyquist sampling do not fully consider the changes of spectral characteristics caused by the sampling structure, and its potential for low computational complexity has yet to be explored, which is the focus of this paper.
In this paper, for the MWC sub-Nyquist sampling structure, we analyzed the special power spectrum characteristics caused by this sampling structure, and proposed a novel reconstruction algorithm, called nearest orthogonal matching pursuit (N-OMP). After an occupied subband is detected, the algorithm uses the correlation coefficient between the residual vector and the sensing matrix column vectors to detect the occupancy of two adjacent subbands directly, thereby reducing the number of iterations of the reconstruction algorithm. Furthermore, we deduced the spectrum occupancy discriminant based on the correlation coefficient and gave its proof process. Simulation experiments verify that our proposed algorithm can significantly reduce the computational complexity while showing better support reconstruction accuracy.
The remainder of this paper is organized as follows. In Section II, we present the system model. Reconstruction algorithm is studied in Section III. Next, we use simulation experiments to verify our proposed algorithm in Section IV. Finally, we come to the conclusions and future work in Section V.

II. SYSTEM MODEL A. MWC SAMPLING
Assuming that the detection signal r(t)is a real-valued continuous time signal with a frequency range of f ∈ F = −f NYQ /2, f NYQ /2 , which conforms to the multi-band model; the number of subbands of r(t) is N ; and the maximum bandwidth of any subband does not exceed B max .
Signal r(t) is first modulated by m different pseudo-random sequences p i (t) with period T p at the same time, i ∈ {1, 2, · · · , m}. Formally, where α ∈ {−1, +1}, p(t) = p(t + nT p ), n ∈ Z, and M is the number of symbols ±1 in a period. And thus we can get x i (t) = r(t)p i (t). Obviously, p i (t) satisfies the Dirichlet condition, and its Fourier series is expressed as where Tp t dt. Let Further, the Fourier transform of x i (t) is where R(f ) = +∞ −∞ r(t)e −j2πft dt. It can be found that X i (f ) in (5) is the infinite weighted sum of R(f ) with B p = 1 T p as the step size and c i (n) as the weight coefficient. Therefore, we only need to use m low-pass filters to filter the signal x i (t) and ensure that the bandwidth B LPF of the low-pass filter is not less than B max , then all the spectrum information of the received signal r(t) can be retained. Subsequently, the output of the low-pass filter i can be expressed as where * means convolution operation.
When f / ∈ F, R(f ) = 0. We only need to consider the number of translation of the received signal spectrum R(f ) in the frequency range F. In other words, the infinite term summation of (7) can be simplified by the finite term summation. Thus we can further get where L 0 represents the number of translation required for shifting the rightmost subband in the frequency domain F to the baseband domain F LPF , which satisfies Finally, the output y i (t) of the low-pass filters i is sampled at Nyquist rate with a period of T s and thus we can obtain discrete samples y i [n]. Then the sampling rate of each sampler is B LPF , and the total sampling rate of the MWC system is mB LPF . If the received signal r(t) is directly sampled at the Nyquist rate, the required sampling rate is (2L 0 + 1)B p and thus sampling compression rate is η = As mentioned above, the MWC sampling block diagram is shown in Fig. 1.

B. COMPRESSED SENSING
Expanding (8) into matrix form, we can get (9), as shown at the bottom of the next page, which can be simplified as where Y(f ) represents the observation vector, G is the sensing matrix, and R(f ) represents the sparse vector and shares a sparse pattern. Since f ∈ F LPF = [−B LPF /2, B LPF /2] is continuous and uncountable, (10) is also called an infinite measurement vectors (IMV) system with sparsity K . That is, K = |S| = | supp(R(f ))|, where |S| = | supp(R(f ))| represents the support of the sparse vector R(f ). In order to discuss the restricted isometry property (RIP) of the sensing matrix, the sensing matrix R(f ) is decomposed. Substituting (1) into (3), we can get Therefore, G can be decomposed into (12), as shown at the bottom of the next page, where each row of matrix A is a pseudo-random sequence of alternating ±1, B is a Vandermonde matrix with maximum rank min (M , 2 L 0 + 1), and C is a diagonal matrix. Therefore, when the appropriate m is selected, the matrix G conforms to the RIP property in VOLUME 8, 2020 the sense of high probability. Under the condition of setting appropriate parameters, when the carrier frequency position information is unknown, in order to ensure signal reconstruction, m ≥ 2 N is required [17].
In order to avoid solving the IMV model in (10), [19] proposed to find the support S by solving the multiple measurement vectors (MMV) problem in the form of (13), and proved that (13) has a unique solution U 0 . The key is that the position and number of non-zero rows of U 0 are consistent with R(f ).
where V is the decomposition matrix of matrix Q, which satisfies Q = VV H and Finally, we can build a compressed sensing model of the MMV problem, as shown in Fig. 2.

III. RECONSTRUCTION ALGORITHM A. POWER SPECTRUM FEATURE ANALYSIS
Through the illustration in Section II, we can abstract the received signal power spectrum as a multi-band signal power spectrum segmentation diagram as shown in Fig. 3. In the figure, assuming that the number of effective translations in a single direction is L 0 = 9, the number of spectrum slices is L = 2L 0 + 1 = 19. The color blocks in the figure represent the power spectra of different primary users. There are 3 primary users in total, and the number of subbands N is 6. Fig. 3 shows the power spectrum of a multi-band signal for a specific situation, and we can find that the power spectrum of a multi-band signal based on MWC sampling has the following two special structural features. (i) According to the step size B p , the frequency domain F = −f NYQ /2, f NYQ /2 is fixedly divided into several slices. When the carrier frequency position of PUs is unknown, the power spectrum of PUs appears randomly for MWC, and this will lead to a high probability of green block (trapezoid), that is, a PU signal has a high probability of occupying two slices at the same time.

C(0)
. . . (ii) In addition, there is a certain probability that the blue block (triangle) in the figure occurs, that is, a PU occupies two slices at the same time, but only a small part of one (the 11th slice in the figure) of the slices. This will cause the 11th slice to be indistinguishable from noise, and will be wrongly judged as unoccupied.

B. IMPROVED RECONSTRUCTION ALGORITHM N-OMP
As a classic reconstruction algorithm in the field of compressed sensing, the reconstruction steps of OMP are as follows: 1 calculate the inner product of each column vector of G and V according to (13), and return the index ξ corresponding to the maximum value of the inner product; 2 then update the residual vector to find the next maximum and update the support until the cut-off condition is satisfied. According to the analysis in Part A of Section III, one frequency band will occupy two slices at the same time with high probability. So the basic strategy is descripted as following: after detecting the occupation of the ξ th slice, the power spectrum slicing characteristics generated by the MWC sampling structure are utilized to directly detect the occupancy of slices ξ − 1 and ξ + 1, then the support and residual vector are updated, and the above process are repeated until the stop condition is met. In order to prove our idea, rewrite the (9) as where g i corresponds to the ith column vector of the sensing matrix G, i = 1, 2, · · · , L, and L = 2L 0 + 1. Unlike the single measurement vector model, we treat each row of the matrix U as an atom in the MMV model. Therefore, we expand the matrix U to u 1 · · · u ξ −1 u ξ u ξ +1 · · · u L T , where u ξ represents the ξ th row of U. Subsequently, substituting the expanded form of the matrices G and U into (13), we can get V = u 1 g 1 + · · · + u ξ −1 g ξ −1 + u ξ g ξ + u ξ +1 g ξ +1 + · · · + u L g L . (15) In order to determine the occupancy of slices ξ − 1 and ξ + 1, our idea is to compare the correlation coefficients between V and g ξ −1 and g ξ +1 . Introducing the Pearson correlation coefficient, we get where r ξ −1 represents the Pearson correlation coefficient between V and g ξ −1 , and r ξ +1 represents the Pearson correlation coefficient between V and g ξ +1 , Cov V, g ξ −1 represents the covariance of V and g ξ −1 , and σ v represents the variance of V. Further, the following decision expressions are established to determine the occupancy of the two positions g ξ −1 and g ξ +1 : where w is the threshold, and its value depends on the distribution of the difference between r ξ −1 and r ξ +1 .
Next, we will prove the first two decision inequalities in (17). Assume that slice ξ + 1 is occupied and ξ − 1 is free. Then u ξ −1 = 0, and u ξ +1 = 0. According to (15), g 1 , g 2 , · · · , g ξ −1 , g ξ , g ξ +1 , · · · , g L are the vector basis. Quoting the concept of cosine similarity, we can obtain The closer the cosine value is to 1, the closer the angle is and the closer the two vectors are. In terms of the unequal relation of (17), the cosine similarity is equal to the Pearson correlation coefficient. Therefore, the comparison between VOLUME 8, 2020 r ξ −1 and r ξ +1 can be equivalent to the comparison between the angles θ ξ −1 = V, g ξ −1 and θ ξ +1 = V, g ξ +1 . Because u ξ −1 = 0 and u ξ +1 = 0, θ ξ −1 = V, g ξ −1 > θ ξ +1 = V, g ξ +1 such that r ξ +1 > r ξ −1 . In order to increase the error tolerance, the threshold parameter w is introduced. When r ξ +1 − r ξ −1 ≤ w, no judgment is made; when r ξ +1 > r ξ −1 + w, ξ + 1 is occupied; and when r ξ −1 > r ξ +1 + w, ξ − 1 is occupied. Subsequently, we focus on the necessity of setting the third judgment inequality in (17). Leading to r ξ +1 − r ξ −1 ≤ w, two main reasons are as follows: 1 there is only noise in slices ξ −1 and ξ +1, resulting in r ξ +1 − r ξ −1 ≤ w; 2 there are primary user signals in both slices ξ −1 and ξ +1, and the power spectra in the two slices are relatively close, resulting in r ξ +1 − r ξ −1 ≤ w. It is impossible to distinguish between 1 and 2 by the difference r ξ +1 − r ξ −1 . Therefore, when r ξ +1 − r ξ −1 ≤ w, no judgment is made. For 2 , the occupied slices ξ − 1 or ξ + 1 can be determined by the maximum value of the inner product of the residual and each column of the sensing matrix in the following iteration, even though which cannot be determined by (17) immediately.
The support reconstruction process is shown in algorithm 1.

C. COMPUTATIONAL COMPLEXITY ANALYSIS
According to [23], we can get the computational complexity of the dth iteration of the OMP algorithm as With reference to this complexity calculation method, the complexity analysis of Algorithm 1 is as follows.
Corresponding to (19), the complexity of finding the maximum value can be ignored, so the calculation complexity of Step 3 is O(mL).
Step 4 is mainly to calculate the correlation coefficients of V and G(:, ξ − 1) and G(:, ξ + 1) respectively, therefore the complexity is O(2 L). When m 2, the complexity of Steps 3 and 4 is approximately O(mL). The complexity of calculating the residual in Step 6 is O mρ + mρ 2 + ρ 3 , where ρ is the length of support of the dth iteration, and its value is changed according to the judgment of (17).
Assuming that the cutoff number of iterations of the OMP algorithm is K , the set R = {1, 2, · · · , K }, and d ∈ R. In the iterative process of the OMP algorithm, the support length ρ of the dth iteration is equal to d, and takes each element in R in turn, so the complexity of the whole process is In our proposed algorithm, the support length ρ of the dth iteration does not need to take each value in R in turn, and it is a jump value. Then suppose that in the iterative process of the N-OMP algorithm proposed in this paper, the value set of the support length ρ is H, which must satisfy H ⊆ R.
Then the computational complexity of the entire process of

Algorithm 1 The Nearest Orthogonal Matching Pursuit Algorithm
Input: the sampling sample y = [y 1 y 2 · · · y m ] T , the sensing matrix G, and the threshold w. Initialization: number of iterations d = 1, restoring support = ∅, cut-off threshold of iteration number K = |supp(R(f ))|.

Steps:
1 Calculate matrix Q and its decomposition matrix V by (10)-(13); 2 Let the residual matrix r = V; 3 Find the maximum value of the modulus of each column of G T r, and assign the position of the corresponding column to ξ ; 4 Calculate the correlation coefficients of the residual matrix r and G(:, ξ − 1) and G(:, ξ + 1) respectively according to (16) to obtain r ξ −1 and r ξ +1 ; (Note: because the residual r is a matrix, r ξ −1 is obtained by averaging the Pearson correlation coefficients of G(:, ξ − 1) and the each column of matrix r. In the same way, r ξ +1 can be obtained.) 5 Judge the occupancy of ξ −1 and ξ +1 adjacent to ξ through discriminant (17), and determine the value of ξ ∈ {∅, ξ − 1, ξ + 1}. At the same time, update the support = ∪ ξ ∪ ξ in combination with Step 3; 6 Update the residual matrix r = V − G(:, )G(: , ) † V, and update the number of iterations d = d + 1; 7 Update the cut-off iteration number. If ξ ∈ {ξ − 1, ξ + 1}, the cut-off iteration number K = K − 1; if ξ ∈ {∅}, K remains unchanged; 8 If the number of iterations d > K , stop the loop; otherwise, return to Step 3; Output: output the finally restored support , then the free frequency band set is¯ = {l, l ∈ L, and l / ∈ }, L = {1, 2, 3, · · · , L}.

the N-OMP algorithm is
Because H ⊆ R, CC N−OMP − sum ≤ CC OMP − sum . Consequently, the computational complexity of each iteration of our proposed algorithm is close to that of OMP, but the N-OMP algorithm can greatly reduce the number of iterations, thereby reducing the complexity of the reconstruction process.

IV. SIMULATION EXPERIMENTS
In the Gaussian white noise scene, the received signal is r(t) = s(t)+n(t), where n(t) represents the additive Gaussian white noise, and s(t) is a multi-band signal composed of PUs In (22), N denotes band number, and sinc(t) = sin(πt)/(πt). In addition, E i , α i , and f i represent the energy coefficient, time offset and carrier frequency of the PU signal respectively. And all the three parameters are generated randomly and follow the uniform distribution of interval (0, 3), (0, 0.9 T r ) and (B max , f NYQ /2 − B max ) respectively, where B max = 50MHz represents the maximum bandwidth of the PUs signal; T r = 1.97µs, denotes the duration of the signal; and f NYQ = 10GHz. The target wideband spectrum is divided into L = 195 slices such that the step size B p = f NYQ L ≈ 51.3MHz. Then the bandwidth of the low-pass filter is B LPF = B p ≈ 51.3MHz. In addition, the low-speed sampling period T s and the period T p of the mixed signal p i (t) are set as T s = T p = L/f NYQ = 0.02µs. Fig. 4 is a unilateral spectrogram of the detection signal in SNR = 10dB and 13 PUs (N = 26), which is used to simulate the actual received signal, where SNR is defined as 10 log s(t) 2 2 / n(t) 2 2 . Monte Carlo method is used to verify the detection performance of the algorithm, and each of the following data points is the statistical mean obtained by running the program for 500 times independently. To ensure signal reconstruction, m ≥2 N is required when the carrier frequency position is unknown [17]. In this paper, the simulation range of parameter N is 2-30. In order to meet the basic conditions for successful recovery of support, we fix m to 100. Fig. 5 shows the comparison of the accuracy of support reconstruction between the proposed algorithm N-OMP and the OMP algorithm when w = 0 and SNR = 10dB. As shown in the figure, compared with the OMP algorithm, the N-OMP algorithm has a certain improvement in the accuracy of support recovery. The reason is that the N-OMP algorithm can avoid the missed detection caused by the feature (ii) of Part A in Section III to a certain extent by comparing the correlation coefficients of the two sub-channels adjacent to ξ . Fig. 6 shows the comparison of the number of iterations between N-OMP and OMP when w = 0 and SNR = 10dB. Because the complexity of the N-OMP and the OMP in a  single iteration process are similar, the number of iterations reflects the total computational complexity of the two algorithms. It can be found from Fig. 6 that when w = 0, compared with the OMP algorithm, the N-OMP algorithm can reduce the computational complexity by about 50%.
Combining Fig. 5 and Fig. 6, we can further find that the N-OMP algorithm can reduce the computational complexity by up to 50% while showing better reconstruction accuracy comparing with the OMP algorithm.
In order to find the optimal threshold w, Fig. 7 depicts the success recovery probability of the support and number of iterations required under different threshold w, when SNR = 10dB and N = 10 are fixed. According to feature (i) of Part A in Section III, a PU signal has a high probability of occupying two slices at the same time. So when the threshold w approaches to 0, it shows a higher probability of success. Even when w = 0, the best recovery accuracy is shown. At this time, whether the correlation coefficient r ξ −1 or r ξ +1 is larger, the corresponding position is judged to be occupied. Therefore, when w = 0, the number of iterations only needs half of the number of frequency bands. This is also the reason why the threshold w of the N-OMP algorithm in Fig. 5 and Fig. 6 is 0. With the increase of w, when w ≥ 0.03, the N-OMP algorithm evolves into the OMP algorithm. The detection performance of the two algorithms is very close,  which is not only reflected in the reconstruction accuracy but also in the number of iterations.
Looking back at Fig. 5, our first impression may be that the two algorithms have a low probability of completely restoring each element of the original support S. The reason is that the power of PUs signals randomly fading and appearing randomly in the frequency domain B max , f NYQ /2 − B max . Therefore, in each process of searching for the support, there will always be a few subbands that the signal energy cannot be distinguished from the noise, resulting in a low probability of the overall successful recovery of the support. However, in each reconstruction process, the ratio of the number of correctly restored atoms to the length of the original support S is relatively high, as shown in Fig. 8. Compared with the OMP algorithm, the N-OMP algorithm also has a certain improvement in the ratio of the number of correctly recovered atoms to the original support length.

V. CONCLUSION AND FUTURE WORK
In this paper, combining MWC sub-Nyquist sampling and compressed sensing technology, we propose a novel compressed sensing reconstruction algorithm. The algorithm uses the power spectrum slice characteristics brought by MWC sampling to determine the occupancy of adjacent sub-channels directly. Theoretical derivation and simulation experiment results show that, compared with the OMP algorithm, the algorithm proposed in this paper can greatly reduce the computational complexity of the compressed sensing reconstruction process and provide better support reconstruction accuracy. However, the N-OMP algorithm is still based on energy detection in essence, and cannot provide sufficient reconstruction accuracy in the case of low signal-to-noise ratio and high sparsity. In future work, we will consider combining the power spectrum features brought by MWC sampling and transform domain analysis to obtain higher resolution features for the determination of the occupancy of certain low SNR channels in the wideband spectrum sensing process, and thus improve the overall support reconstruction accuracy.