Wideband Spectrum Sensing Based on Single-Channel Sub-Nyquist Sampling for Cognitive Radio

Abstract

Spectrum sensing is an important task in cognitive radio. However, currently available Analog-to-Digital Converters (ADC) can hardly satisfy the sampling rate requirement for wideband signals. Even with such an ADC, the cost is extremely high in terms of price and power consumption. In this paper, we propose a spectrum-sensing method based on single-channel sub-Nyquist sampling. Firstly, a serial Multi-Coset Sampling (MCS) structure is designed to avoid mismatches among sub-ADCs in the traditional parallel MCS. Clocks of the sample/hold and ADC are provided by two non-uniform sampling clocks. The cooperation between these two non-uniform sampling clocks shifts the high sampling rate burden from the ADC to the sample/hold. Secondly, a power spectrum estimation method using sub-Nyquist samples is introduced, and an efficient spectrum-sensing algorithm is proposed. By exploiting the frequency-smoothing property, the proposed efficient spectrum-sensing algorithm only needs to estimate power spectrum at partial frequency bins to conduct spectrum sensing, which will save a large amount of computational cost. Finally, the sampling pattern design of the proposed serial MCS is given, and it is proved to be a minimal circular sparse ruler with an additional constraint. Simulations show that mismatches in traditional parallel MCS have a serious impact on spectrum-sensing performance, while the proposed serial MCS combined with the efficient spectrum-sensing algorithm exhibits outstanding spectrum-sensing performance at much lower computational cost.

1. Introduction

With the increase of wireless communication, spectrum resources have become increasingly scarce. However, most primary users, or authorized users, do not use the spectrum resources allocated to them most of the time [1]. For example, RFeye node [2] is used to monitor the wireless spectrum at Queen Mary University of London, and most of the spectrum resources are found to be idle [3]. TeleVision White Space (TVWS) in the Ultra High Frequency (UHF) band is the spectrum band not used by Digital Terrestrial Televisions (DTT), and Programme Making and Special Events (PMSE). The wireless communication network using TVWS can cover above 10 km in diameter, and it has strong penetration capabilities for buildings, mountains and forests. Cognitive radio devices working in TVWS can be used in rural areas where the communication infrastructure is imperfect to build a communication link among machines and improve residents’ access to the Internet. The communication regulator in the United Kingdom, the Office of Communications (Ofcom), announced that using TVWS does not need a communication license in December 2015 [4]. Therefore, TVWS is one of the most promising opportunistic spectrum access resources.

TVWS can be used by secondary users to transmit their information to enhance spectrum utilization efficiency. Geo-location database is a way to discover available TVWS [5]. However, this method requires a wired or wireless connection with the Geo-location database in advance [6]. The dynamic changes in the wireless environment pose a huge challenge to this database-based approach. Therefore, we need to quickly and accurately detect whether primary users exist before using spectrum resources in TV band. According to the well-known Nyquist sampling theorem, alias free sampling requires the system sampling rate be higher than twice the maximum frequency of the signal. The currently available Analog-to-Digital Converters (ADC) often cannot satisfy the sampling rate requirement for wideband signals. Even with such an ADC, the cost is usually extremely high in terms of price and power consumption. Thus, it is attractive to conduct spectrum sensing based on sub-Nyquist sampling techniques.

Some parallel sub-Nyquist sampling schemes have been developed in recent years, such as the Modulated Wideband Converter (MWC) [7] and the parallel Multi-Coset Sampling (MCS) [8]. A Distributed Modulated Wideband Converter (DMWC) scheme is proposed, and a sensor node is considered to be an MWC channel in DMWC [9]. DMWC combines cooperative spectrum sensing and MWC techniques for wideband spectrum sensing. A Random Triggering-based Modulated Wideband Compressive Sampling (RT-MWCS) is also designed for repetitive sparse wideband signals [10]. However, the input signal is multiplied with the periodic Pseudo-Random Binary Sequence (PRBS) by a mixer in MWC, and it is difficult to suppress the unwanted harmonics due to the usage of PRBS instead of a single-tone as the local oscillator signal. Parallel MCS is another promising sub-Nyquist sampling scheme which has several sub-ADCs working parallelly. The sampling clocks of these sub-ADCs have the same frequency but different phases. Time-Interleaved ADC (TIADC) platform can be used to implement parallel MCS to further improve the time resolution of TIADC. Just as mismatches among sub-ADCs, mainly offset mismatch, gain mismatch and sampling clock skew, degrade the Spurious Free Dynamic Range (SFDR) of a TIADC system [11,12], they can also severely deteriorate spectrum-sensing performance based on parallel MCS.

The mismatch problem in parallel MCS can be avoided in single-channel sampling structure. A single-channel Non-Uniform Sampling (NUS) structure is design, in which a custom-designed Sample/Hold (S/H) and a commercial ADC are used [13]. Then a nonlinear compressive sensing algorithm is used to reconstruct the signal. The computational complexity of reconstruction is reduced by emulating MCS in NUS [14]. The signal sparsity is exploited to reduce the sampling rate requirement in [7,8,13,14], and the Nyquist sampling rate is required when the signal is not sparse.

The above mentioned sub-Nyquist sampling methods concentrate on signal reconstruction, while one may only be interested in the power spectrum in some practical applications, such as cognitive radio and spectrum monitoring. A Single-Channel Modulated Wideband Converter (SCMWC) scheme is proposed, and experiment circuit is designed using commercial Integrated Circuits (ICs) to verify the feasibility of SCMWC [15]. Only power spectrum information is estimated to conduct spectrum sensing in [15]. The signal statistical property is exploited to estimate finite resolution power spectrum of signal without requiring the signal to be sparse [16]. The mean and variance of the finite resolution power spectrum estimator are investigated when the number of sub-Nyquist samples is limited [17]. The signal statistical property is also utilized to estimate power spectrum of signal in [18]. After the power spectrum is obtained, the primary user is detected using the energy detection method. However, the estimation of power spectrum requires a large amount of calculation. Compared with the multi-band model [18], authors in [19] consider both multi-band and multi-tone models, the case of sparse and non-sparse, and known and unknown carrier frequencies. They derive minimum sampling rate requirement for each setting without noise and provide power spectrum reconstruction techniques to achieve these minimum sampling rates. The least squares method is used to estimate the signal power spectrum, and the full rank condition of the system matrix is viewed as a circular sparse ruler design problem [20,21].

In this paper, we propose a spectrum-sensing method based on single-channel sub-Nyquist sampling. To begin with, to avoid the mismatch problem among sub-ADCs in parallel MCS [16,17,18], a serial MCS structure is proposed. The proposed serial MCS mainly consists of a S/H, an ADC and two different periodic non-uniform sampling clocks. The clocks of S/H and ADC are provided by these two different periodic non-uniform sampling clocks. The S/H functions similar to an analog signal register, which can shift the high sampling rate requirement for wideband signals from ADC to S/H. Although the sampling rate requirement for the ADC in serial MCS is higher than that in parallel MCS, there is no mismatch problem and the implementation size is smaller. Moreover, to reduce the computational complexity in [18], we propose an efficient spectrum-sensing algorithm. The frequency-smoothing property of the power spectrum estimator is exploited to estimate the frequency-smoothed power spectrum of the signal. Since the frequency bins in a frequency-smoothing window are highly correlated, spectrum sensing only needs estimate power spectrum at partial frequency bins, which will save a large amount of computational cost. Then finally, to use the least square algorithm to estimate the power spectrum, we must ensure that the system matrix has full column rank. Whether the system matrix has full column rank is determined by the sampling pattern design, which is further proved to be a minimal circular sparse ruler design problem with an additional constraint. The sampling pattern design method in [20,21] is only suitable for the traditional parallel MCS. We consider the sampling pattern design for the proposed serial MCS. Simulations show that mismatches among sub-ADCs deteriorate the spectrum-sensing performance seriously, especially the false alarm probability, while the proposed MCS combined with the proposed efficient spectrum-sensing algorithm not only has outstanding spectrum-sensing performance, but also requires less computational load.

This paper is organized as follows: Section 2 presents the signal model and forms the problem. The sampling structure and sampling process of the proposed serial MCS are described in Section 3. Section 4 firstly introduces a power spectrum estimation method and then proposes an efficient spectrum-sensing algorithm. The sampling pattern design of the proposed serial MCS is developed in Section 5. Simulation results are given in Section 6. Simulations take some practical parameters into account, such as Signal-to-Noise Ratio (SNR), compression ratio, and acquisition time. Section 7 discusses using the spectrum-sensing module to receive communication signals.

2. Signal Model and Problem Formulation

The spectrum resource we are interested in is the UHF TV band, and its frequency range is from 470 MHz to 790 MHz. The bandwidth of a channel is 8 MHz, and there are 40 channels in total. Assume that $M$ users are active during spectrum sensing. $M$ users may comprise of both primary users and other secondary users. Without loss of generality, the k-th user can be represented as:

$$u_k\left(t\right)={{\displaystyle \sum }}_{n\in \mathbb{Z}}{d}_k[n]g_k\left(t-n{T}_k\right){\mathrm{e}}^{j2\pi f_{k}t} , k=1,2,\dots ,M,$$

(1)

where $\left\{d_k[n]\right\}$ is a sequence of transmitted symbol, $g_k\left(t\right)$ is the pulse-shaping function, $T_k$ is the symbol duration and $f_k$ is the carrier frequency.

The practical wireless transmission channel mainly has two effects on the signal, one is log-normal shadowing which is a typical large-scale fading, the other is Rayleigh fading which is a typical small-scale fading [22,23,24]. To begin with, Large-scale fading denotes the average signal power attenuation caused by a wide range of movements. This kind of fading is mainly caused by the prominent terrain between the transmitter and the receiver. Experimental results show that the large-scale fading can be represented by a log-normal random variable [23]. Thus, large-scale fading can be expressed as:

$$PL=P{L}_0+10\gamma \log \frac{d}{d_0}+X_{\rho }$$

(2)

where the reference distance $d_0$ is in the far field of the transmitting antenna, the reference power $P{L}_0$ is usually obtained by field test, the path loss exponent $\gamma$ depends on the carrier frequency, antenna height and transmission environment, $X_{\rho }$ denotes a Gaussian random variable with variance ${\rho }^2$ and unit dB. Thus, the received signal affected by large-scale fading can be represented as:

$$s_k\left(t\right)={10}^{-P{L}_k/20}·u_k\left(t\right)$$

(3)

where $P{L}_k$ denotes the attenuation between the $k$-th primary user and the cognitive radio device. Moreover, the primary user may arrive at the cognitive radio device via multiple paths, which can be modeled as Rayleigh fading. In other words, the envelope of the signal subjects to the Rayleigh distribution. The probability density function of Rayleigh distribution can be expressed as:

$$p_h\left(r\right)=\{\begin{array}{l}\frac{r}{{\lambda }^2}{\mathrm{e}}^{-r^2/2{\lambda }^2}, r\ge 0\\ 0, \mathrm{otherwise}\end{array}$$

(4)

where $r$ denotes the amplitude of the envelop of the received signal, $2{\lambda }^2$ denotes its average power. The $k$-th primary user $s_k\left(t\right)$ arrived at the cognitive radio device can be written as

$$s_k\left(t\right)={10}^{-P{L}_k/20}·u_k\left(t\right)\ast h_k\left(t\right),$$

(5)

where $\ast$ and $h_k\left(t\right)$ denote the convolution operator and channel impulse response.

Then the signal received by the cognitive radio device is:

$$x\left(t\right)=s\left(t\right)+n\left(t\right)={\sum }_{k=1}^M{s}_k\left(t\right)+n\left(t\right),$$

(6)

where $s\left(t\right)$ is the signal component and $n\left(t\right)$ is the additive white Gaussian noise. Figure 1 shows a wideband signal with 3 active users. They occupy the 23rd, 40th and 59th channel, while the other channels are vacant. These vacant channels can be used by secondary users to transmit their information.

After the spectrum range and channel partitioning format in the UHF TV band is given, the problem of wideband spectrum sensing is to detect which channels are occupied by the primary users and which channels are vacant. A band-pass filter is firstly used to filter out the noise outside the band [470 MHz, 790 MHz], and then a mixer is utilized to demodulate the filtered signal. The demodulated signal is a band-limited signal and its bandwidth is 320 MHz. A traditional method to solve this problem is to first sample the demodulated signal at the Nyquist sampling rate, and then estimate its power spectrum. Finally, the energy detection method is used to determine whether primary users exist. The Nyquist sampling theorem states that alias-free sampling requires the sampling rate be higher than twice the maximum frequency of the signal. Therefore, the sampling rate requirement for the demodulated signal is 640 MHz, which poses a great challenge to current available ADCs. Even with such a high sampling rate ADC, it is usually unacceptable in terms of price and power consumption for some sensors and battery-powered devices. In this paper, we consider using sub-Nyquist sampling technique to conduct spectrum sensing.

3. Single-Channel Sub-Nyquist Sampling Structure

An overview of the proposed single-channel sub-Nyquist sampling scheme is shown in Figure 2. The signal received by the antenna is first band-pass filtered to remove the out-of-band noise and interference, and then a mixer is used to demodulate the filtered signal. Finally, low-pass filter is used to remove the local oscillator signal. The demodulated signal still has a high bandwidth if we are interested in many channels. The dotted box shows the proposed serial MCS. The S/H is used to sample and temporarily store the input analog signal $x\left(t\right)$. Then the subsequent ADC can sample it at a lower rate. Two periodic Pseudo-Random Binary Sequences are used as the non-uniform sampling clocks of S/H and ADC. The periodic non-uniform sampling clock provided to S/H is determined by the sampling pattern of MCS, while the periodic non-uniform sampling clock provided to ADC is dependent on not only the sampling pattern but also the sampling rate of ADC. The spectrum-sensing procedure can be summarized as follows: (1) Classifying the periodic non-uniformly sampled data into several uniform sampling sequences; (2) Power spectrum at partial frequency bins is estimated using these uniform sampling sequences; (3) Compare the estimated power spectrum with a pre-defined threshold to determine whether the primary user exists.

In detail, let $x\left(t\right)$ denote the signal output from the low-pass filter. Let the Nyquist sampling rate and interval of $x\left(t\right)$ be $f_{\mathrm{N}}$ and $T$, respectively. Let the sampling rate and interval of ADC be $f_{\mathrm{ADC}}$ and $T_{\mathrm{ADC}}$, respectively. Let $d$ represent the sub-Nyquist sampling factor which is defined as $d=f_{\mathrm{N}}/f_{\mathrm{ADC}}$. $q$ samples are selected from $L$ consecutive Nyquist samples in MCS. The sampling pattern $C=\left\{c_1,c_2,\dots ,c_q\right\}$, where $1\le c_1<c_2<\dots <c_q\le L$ is used as the sample index. Let $N_0$ represent the number of samples of a coset, and the total number of periodic non-uniformly sampled data is $N=q\times N_0$. We will explain in Section 5, Sampling pattern design, that at least one pair of cosets with only $1\times T$ time difference should be included in the sampling pattern. The ADC cannot directly capture two samples $c_i$ and $c_{i+1}$ with only $1\times T$ time difference because the sampling interval of ADC $T_{\mathrm{ADC}}=d\times T$ is greater than $1\times T$. The function of S/H is to temporarily store the second coset $c_{i+1}$ until the ADC has finished capturing the previous coset $c_i$.

An example of the sampling process is given in Figure 3, where $q=6$, $L=24$, $f_{\mathrm{ADC}}=f_{\mathrm{N}}/4$ and $C=\left\{1,2,9,12,14,18\right\}$. The first coset $C\left(1\right)$ can be sampled directly, while the second coset $C\left(2\right)$ must be firstly sampled by the S/H and subsequently captured by the following ADC on the fifth Nyquist time grid. Thus, the third coset $C\left(3\right)$ must be greater than or equal to 5, otherwise, it will lose the last coset $C\left(2\right)$. The cooperation between the first and the second periodic non-uniform sampling clock ensures that two adjacent cosets with only $1\times T$ time difference can be captured by an ADC with sampling rate $f_{\mathrm{ADC}}$ which is much lower than $f_{\mathrm{N}}$, and this is the essence of the proposed serial MCS.

The periodic non-uniformly sampled data are classified into several uniform sampling sequences. Let $y[n],1\le n\le N$ represent the periodic non-uniformly sampled data output from the ADC and it can be written as:

$$y[n]=[y_{c1}[1],y_{c2}[1],\dots ,y_{cq}[1],y_{c1}[2],\dots ,y_{cq}[N_0]],$$

(7)

where $y_{ci}[k], 1\le i\le q, 1\le k\le N_0$ is the sample obtained from the k-th Nyquist sample block with length L using sample index $C\left(i\right)$. Please note that the number of $y[n]$ equals to the total number of periodic non-uniformly sampled data $N$. Then periodic non-uniformly sampled data $y[n],1\le n\le N$ is classified into $q$ uniform sampling sequences

$$y_{ci}[k]=y[\left(k-1\right)\times q+i], 1\le k\le N_0, 1\le i\le q.$$

(8)

Let ${\mathit{y}}_{ci}={[y_{ci}[1],y_{ci}[2],\dots ,y_{ci}[N_0]]}^T$ represent the i-th uniform sampling sequence which corresponds to the i-th coset $c_i$. The sampling rate and time offset of ${\mathit{y}}_{ci}$ are $f_{\mathrm{N}}/L$ and $ci\times T$, respectively. $q$ uniform sampling sequences can be written in matrix form concisely:

$${\mathit{y}}_c={[{\mathit{y}}_{c1},{\mathit{y}}_{c2},\dots ,{\mathit{y}}_{cq}]}^{\mathrm{T}}.$$

(9)

The parameter of current available ICs should be considered in practical implementation. The sampling rate requirement for S/H is high, but the sampling rate requirement for the following ADC is low. The high sampling rate burden for wideband signals is shifted from ADC to S/H. The sampling rate and input bandwidth of current state-of-the-art S/H are 4GSPS (Giga-Samples Per Second) and 18 GHz, respectively (HMC760LC4B, Analog Devices, Inc., Norwood, MA, USA). Moreover, the jitter of the periodic non-uniform sampling clock is also an important parameter in serial MCS. The alternating rate and jitter of current state-of-the-art PRBS generator are 80 GHz and 600 fs (femtosecond), respectively [25]. Periodic non-uniform sampling clock can also be generated using commercial parallel-to-serial converter. The alternating rate and typical jitter of current state-of-art parallel-to-serial converter are 3.2 GHz and 200 fs, respectively (MC100EP446, ON Semiconductor, Inc., Phoenix, AZ, USA).

4. Spectrum Sensing

After $q$ uniform sampling sequences are obtained using the proposed serial MCS, they are further processed to conduct spectrum sensing. We firstly introduce a power spectrum estimation method [18], and then an efficient spectrum-sensing algorithm is proposed.

4.1. Power Spectrum Estimation

This section mainly consists of two parts:

(a)

Establish the theoretical model for power spectrum estimation.

(b)

Develop an implementation method for the theoretical model.

(a). Establish the Theoretical Model for Power Spectrum Estimation

Let $X\left(f\right)$ denote the Fourier Transform of the demodulated signal $x\left(t\right)$ in Figure 2. We first analyze the relationship between the spectrum of the i-th uniform sampling sequence ${\mathit{y}}_{ci}$ and the original spectrum $X\left(f\right)$. Define $Y_{ci}\left({\mathrm{e}}^{j2\pi fLT}\right)$ as the Discrete Time Fourier Transform (DTFT) of the i-th uniform sampling sequence ${\mathit{y}}_{ci}$, and it can be written as:

$$\begin{array}{lll}{Y}_{ci}\left({\mathrm{e}}^{j2\pi fLT}\right)& =& {\sum }_{n=-\infty }^{+\infty }{y}_{ci}[n]{\mathrm{e}}^{-j2\pi fnLT}\\ & =& \frac{1}{LT}{\sum }_{l=-\frac{L}{2}}^{\frac{L}{2}-1}{X}_{ci}\left(f-\frac{l}{LT}\right)\\ & =& \frac{1}{LT}{\sum }_{l=-\frac{L}{2}}^{\frac{L}{2}-1}X\left(f-\frac{l}{LT}\right){\mathrm{e}}^{j2\pi c_{i}T\left(f-\frac{l}{LT}\right)}\\ & =& \frac{1}{LT}{\sum }_{l=-\frac{L}{2}}^{\frac{L}{2}-1}X(f-\frac{l}{LT}){\mathrm{e}}^{-j2\pi c_i\frac{l}{L}}{\mathrm{e}}^{j2\pi f{c}_{i}T}, f\in {\mathcal{F}}_0,\end{array}$$

(10)

where ${\mathcal{F}}_0=[0,1/\left(LT\right))$ and $L$ is assumed to be even. The band-limited property $X\left(f\right)=0$, $|f|>1/\left(2T\right)$ is used in the derivation from the first step to the second step. The relationship between the spectrum of the delayed signal $X_{ci}\left(f\right)$ and the original spectrum $X\left(f\right)$, $X_{ci}\left(f\right)=X\left(f\right){\mathrm{e}}^{j2\pi f{c}_{i}T}$, is used in the derivation from the second step to the third step. The term ${\mathrm{e}}^{j2\pi f{c}_{i}T}$ in Equation (10) can be moved to the left of the equation, and define

$$Y_i\left({\mathrm{e}}^{j2\pi fLT}\right)\triangleq Y_{ci}\left({\mathrm{e}}^{j2\pi fLT}\right){\mathrm{e}}^{-j2\pi f{c}_{i}T}.$$

(11)

Equation (10) can be rewritten as

$$Y_i\left({\mathrm{e}}^{j2\pi fLT}\right)=\frac{1}{LT}{\sum }_{l=-\frac{L}{2}}^{\frac{L}{2}-1}X(f-\frac{l}{LT}){\mathrm{e}}^{-j2\pi c_i\frac{l}{L}}, f\in {\mathcal{F}}_0.$$

(12)

At this point, we have established the relationship between the spectrum of the i-th uniform sampling sequence ${\mathit{y}}_{ci}$ and the original spectrum $X\left(f\right)$ in Equation (12) for single uniform sampling sequence ${\mathit{y}}_{ci}$. $q$ uniform sampling sequences, or the entire MCS system, can be expressed in matrix form concisely

$$\mathit{Y}\left({\mathrm{e}}^{j2\pi fLT}\right)=A\mathit{x}\left(f\right), f\in {\mathcal{F}}_0,$$

(13)

where $\mathit{Y}{\left({\mathrm{e}}^{j2\pi fLT}\right)}_i=Y_i\left({\mathrm{e}}^{j2\pi fLT}\right)$, $A_{il}={\mathrm{e}}^{-j2\pi c_i\frac{l}{L}}$, $\mathit{x}{\left(f\right)}_l=X\left(f-\frac{l}{LT}\right)$, $1\le i\le q$ and $-\frac{L}{2}\le l\le \frac{L}{2}-1$.

The demodulated signal $x\left(t\right)$ comprises of signal component $s\left(t\right)$ and noise component $n\left(t\right)$. Let $S\left(f\right)$ and $N\left(f\right)$ denote the Fourier Transform of $s\left(t\right)$ and $n\left(t\right)$, respectively. Define $L\times 1$ vector $\mathit{s}{\left(f\right)}_l=S\left(f-\frac{l}{LT}\right),-\frac{L}{2}\le l\le \frac{L}{2}-1$ and $L\times 1$ vector $\mathit{n}{\left(f\right)}_l=N\left(f-\frac{l}{LT}\right),-\frac{L}{2}\le l\le \frac{L}{2}-1$. Let ${\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)$ denote the covariance matrix of $\mathit{Y}\left({\mathrm{e}}^{j2\pi fLT}\right)$, and it can be written as:

$${\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)=\mathrm{E}\left\{\mathit{Y}\left({\mathrm{e}}^{j2\pi fLT}\right){\mathit{Y}}^{\mathit{H}}\left({\mathrm{e}}^{j2\pi fLT}\right)\right\}.$$

(14)

Similarly, the covariance matrix of $\mathit{x}\left(f\right)$, $\mathit{s}\left(f\right)$ and $\mathit{n}\left(f\right)$ can be expressed as ${\mathrm{R}}_{\mathit{x}}\left(f\right)=\mathrm{E}\left\{\mathit{x}\left(f\right){\mathit{x}}^{\mathit{H}}\left(f\right)\right\}$, ${\mathrm{R}}_{\mathit{s}}\left(f\right)=\mathrm{E}\left\{\mathit{s}\left(f\right){\mathit{s}}^{\mathit{H}}\left(f\right)\right\}$ and ${\mathrm{R}}_{\mathit{n}}\left(f\right)=\mathrm{E}\left\{\mathit{n}\left(f\right){\mathit{n}}^{\mathit{H}}\left(f\right)\right\}$. Equation (14) can be further written as:

$${\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)=A{\mathrm{R}}_{\mathit{x}}\left(f\right)A^H=A[{\mathrm{R}}_{\mathit{s}}\left(f\right)+{\mathrm{R}}_{\mathit{n}}\left(f\right)]A^H.$$

(15)

The following explanation shows that the covariance matrix ${\mathrm{R}}_{\mathit{x}}\left(f\right)$ is a diagonal matrix. The (i,j)-th element of ${\mathrm{R}}_{\mathit{x}}\left(f\right)$ can be expressed as:

$$\begin{array}{lll}{[{\mathrm{R}}_{\mathit{x}}(f)]}_{i,j}& =& \mathrm{E}[\mathit{x}{(f)}_i\mathit{x}{(f)}_j^{*}]\\ & =& \mathrm{E}\{[\mathit{s}{\left(f\right)}_i+\mathit{n}{\left(f\right)}_i]{[\mathit{s}{\left(f\right)}_j+\mathit{n}{\left(f\right)}_j]}^{*}\}\\ & =& \mathrm{E}[\mathit{s}{\left(f\right)}_i\mathit{s}{\left(f\right)}_j^{*}]+\mathrm{E}[\mathit{n}{\left(f\right)}_i\mathit{n}{\left(f\right)}_j^{*}]\\ & =& \mathrm{E}[{|\mathit{s}{\left(f\right)}_i|}^2]\delta (i-j)+{\sigma }^2\delta (i-j),\end{array}$$

(16)

where the signal component $s\left(t\right)$ and noise component $n\left(t\right)$ are assumed to be independent, different primary users are assumed to be independent, and the property of stationary process is used in these derivations. The first term $\mathrm{E}[{|\mathit{s}{\left(f\right)}_i|}^2]\delta \left(i-j\right)$ in Equation (16) corresponds to a diagonal matrix with only $M$ non-zero elements, while the second term ${\sigma }^2\delta \left(i-j\right)$ in Equation (16) also corresponds to a diagonal matrix but with $L$ non-zero elements. Therefore, the covariance matrix ${\mathrm{R}}_{\mathit{x}}\left(f\right)$ is a diagonal matrix with only $L$ non-zero diagonal elements. This characteristic is used to further develop the power spectrum estimation method.

Let $\mathrm{Vec}(·)$ and $\otimes$ represent the matrix vectorization operator and Kronecker product, respectively. Define ${\mathrm{r}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)\triangleq \mathrm{Vec}[{\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)]$, Equation (15) can be rewritten as:

$$\begin{array}{lll}{\mathrm{r}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)& =& \mathrm{Vec}[{\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)]\\ & =& [{\left(A^H\right)}^{\mathrm{T}}\otimes A]\mathrm{Vec}[{\mathrm{R}}_{\mathit{x}}\left(f\right)]\\ & =& \left(A^{*}\otimes A\right)\mathrm{Vec}[{\mathrm{R}}_{\mathit{x}}\left(f\right)],\end{array}$$

(17)

where the matrix identity $\mathrm{Vec}\left(AXB\right)=\left(B^{\mathrm{T}}\otimes A\right)\mathrm{Vec}\left(X\right)$ is used. The vector form of the diagonal matrix ${\mathrm{R}}_{\mathit{x}}\left(f\right)$ can be further written as $\mathrm{Vec}[{\mathrm{R}}_{\mathit{x}}\left(f\right)]={[P_1\left(f\right)e_1^T,P_2\left(f\right)e_2^T,\dots ,P_L\left(f\right)e_L^T]}^T$, where $P_i\left(f\right)$ is the i-th diagonal element and $e_i$ is selected from the i-th row of a $L\times L$ identity matrix. Further define ${\tilde{\mathrm{r}}}_{\mathit{x}}\left(f\right)={[P_1\left(f\right),P_2\left(f\right),\dots ,P_L\left(f\right)]}^T$ and a $L^2\times L$ selection matrix $B$ which has 1 at the j-th column and the $[\left(j-1\right)L+j]$-th row. The vector form of the diagonal matrix ${\mathrm{R}}_{\mathit{x}}\left(f\right)$ can be written as $\mathrm{Vec}[{\mathrm{R}}_{\mathit{x}}\left(f\right)]={\mathrm{B}\tilde{\mathrm{r}}}_{\mathit{x}}\left(f\right)$. Equation (17) can be further expressed as

$$\begin{array}{lll}{\mathrm{r}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)& =& \left(A^{*}\otimes A\right)\mathrm{Vec}[{\mathrm{R}}_{\mathit{x}}\left(f\right)]\\ & =& \left(A^{*}\otimes A\right){\mathrm{B}\tilde{\mathrm{r}}}_{\mathit{x}}\left(f\right)\\ & =& \Phi {\tilde{\mathrm{r}}}_{\mathit{x}}\left(f\right),\end{array}$$

(18)

where $\Phi =\left(A^{*}\otimes A\right)B$ is a $q^2\times L$ matrix. If $\Phi$ has full column rank and Equation (18) is an over-determined system, then least square algorithm can be used to estimate the power spectrum

$${\tilde{\mathrm{r}}}_{\mathit{x}}\left(f\right)={\Phi }^{†}{\mathrm{r}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right),$$

(19)

where ${\Phi }^{†}$ is the Moore-Penrose inverse of $\Phi$.

(b). Develop an Implementation Method for the Theoretical Model

The theoretical model of power spectrum estimation is established in the previous section. However, the Fourier transform involved in the theoretical model is implemented through Fast Fourier Transform (FFT) in practice. In this section, the implementation method of the theoretical model in the digital domain is discussed. The main contents include: (1) Compensate the fractional delay of the uniform sampling sequences in MCS. (2) Estimate the cross-power spectrum of the uniform sampling sequences using the frequency-smoothed cyclic cross-periodogram. (3) Detect the presence of the primary user using the estimated power spectrum.

To begin with, let $Y_{ci}[k]=Y_{ci}\left({\mathrm{e}}^{j2\pi fLT}\right){|}_{f=\frac{k}{LT{N}_0}},1\le k\le N_0$ represent the FFT of the i-th uniform sampling sequence. Define a constant phasor corresponding to the i-th uniform sampling sequence $P_i[k]={\mathrm{e}}^{-j2\pi f{c}_{i}T}{|}_{f=\frac{k}{LT{N}_0}},1\le k\le N_0$. Equation (11) in the theoretical model compensates the fractional delay of the i-th uniform sampling sequence of MCS. In practice, the implementation of fractional delay compensation in the digital domain can be represented as:

$$Y_i[k]=P_i[k]Y_{ci}[k], 1\le k\le N_0$$

(20)

where $Y_i[k]=Y_i\left({\mathrm{e}}^{j2\pi fLT}\right){|}_{f=\frac{k}{LT{N}_0}}$. Define $x[n]\triangleq x\left(nT\right)$ and its FFT $X[k]\triangleq \frac{1}{\sqrt{N_{0}L}}{\sum }_{n=1}^{N_{0}L}x[n]{\mathrm{e}}^{-j2\pi nk/N_{0}L}$. Partition $X[k]$ into $X_i[k]\triangleq X[k+\left(i-1\right)N_0]$, $1\le k\le N_0$ and $1\le i\le L$. Define $\mathit{x}[k]={[X_1[k],X_2[k],\dots ,X_L[k]]}^T$ and $\mathit{Y}[k]={[Y_1[k],Y_2[k],\dots ,Y_q[k]]}^T$. Equation (13) in the theoretical model can be written as $\mathit{Y}[k]=A\mathit{x}[k],1\le k\le N_0$.

Next, to estimate the signal power spectrum using Equation (19), we should first obtain ${\mathrm{r}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)=\mathrm{Vec}[{\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)],f\in {\mathcal{F}}_0$. The element at the i-th row and j-th column of matrix ${\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)$ is ${[{\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)]}_{i,j}=\mathrm{E}\left\{\mathit{Y}{\left({\mathrm{e}}^{j2\pi fLT}\right)}_i\mathit{Y}{\left({\mathrm{e}}^{j2\pi fLT}\right)}_j^{*}\right\},f\in {\mathcal{F}}_0,1 \le i,j\le q$. That is, in order to estimate the power spectrum of signal using Equation (19), we should first calculate the cross-power spectrum of $\mathit{Y}{\left({\mathrm{e}}^{j2\pi fLT}\right)}_i=Y_i\left({\mathrm{e}}^{j2\pi fLT}\right)$ and $\mathit{Y}{\left({\mathrm{e}}^{j2\pi fLT}\right)}_j=Y_j\left({\mathrm{e}}^{j2\pi fLT}\right)$ at frequency bin $f$, $f\in {\mathcal{F}}_0$. A cross-power spectrum estimator [26] is introduced to calculate the cross-power spectrum of $Y_i\left({\mathrm{e}}^{j2\pi fLT}\right)$ and $Y_j\left({\mathrm{e}}^{j2\pi fLT}\right)$. Define $S_{Z_i{Z}_j}\left(\alpha ;f\right)\triangleq \mathrm{E}\{Z_i\left(f\right)Z_j^{*}\left(f-\alpha \right)\}$ as the cyclic cross-spectrum of $Z_i\left(f\right)$ and $Z_j\left(f\right)$, where $\alpha$ is the cyclic frequency. The frequency-smoothed cyclic cross-periodogram is used as a consistent estimator of $S_{Z_i{Z}_j}\left(\alpha ;f\right)$, and it can be written as

$${\widehat{S}}_{Z_i,Z_j}\left(\alpha ;f\right)\triangleq \frac{1}{N_0}{\sum }_{n=1}^{N_0}w\left({\mathrm{e}}^{j2\pi \left(f-\frac{n}{N_{0}LT}\right)LT}\right)Z_i({\mathrm{e}}^{j2\pi \frac{n}{N_0}})Z_j^{*}({\mathrm{e}}^{j2\pi \left(\frac{n}{N_{0}LT}-\alpha \right)LT}),$$

(21)

where $Z_i({\mathrm{e}}^{j2\pi \frac{n}{N_0}})=Z_i({\mathrm{e}}^{j2\pi fLT}){|}_{f=\frac{n}{N_{0}LT}}$, $Z_i\left({\mathrm{e}}^{j2\pi fLT}\right)$ is the DTFT of $z_i[n]$ and $w\left({\mathrm{e}}^{j2\pi fLT}\right)$ is the spectral window function. In practical digital implementation, FFT is used to calculate the DTFT. The frequency-smoothed cyclic cross-periodogram can be correspondingly written as

$${\widehat{S}}_{Z_i,Z_j}\left(\alpha =\frac{m}{N_{0}LT};f=\frac{k}{N_{0}LT}\right)=\frac{1}{N_0}{\sum }_{n=1}^{N_0}w[k-n]Z_i[n]Z_j[n-m],$$

(22)

where $w[n]\triangleq w\left({\mathrm{e}}^{j2\pi fLT}\right){|}_{f=n/N_{0}LT}$ and $Z_i[n]\triangleq Z_i\left({\mathrm{e}}^{j2\pi fLT}\right){|}_{f=n/N_{0}LT}$. Define ${[{\mathrm{R}}_{\mathit{Y}}[k]]}_{i,j}$ as ${[{\mathrm{R}}_{\mathit{Y}}\left({\mathrm{e}}^{j2\pi fLT}\right)]}_{i,j}{|}_{f=\frac{k}{LT{N}_0}}=\mathrm{E}\{Y_i({\mathrm{e}}^{j2\pi fLT})Y_j^{*}({\mathrm{e}}^{j2\pi fLT})\}{|}_{f=\frac{k}{LT{N}_0}}$, and it can be estimated using

$${[{\mathrm{R}}_{\mathit{Y}}[k]]}_{i,j}={\widehat{S}}_{Y_i,Y_j}\left(\alpha =0;f=\frac{k}{N_{0}LT}\right)=\frac{1}{N_0}{\sum }_{n=1}^{N_0}w[k-n]Y_i[n]Y_j^{*}[n].$$

(23)

Thus, the frequency-smoothed cyclic cross-periodogram of ${\mathrm{R}}_{\mathit{Y}}[k]$ can be written as:

$$\begin{array}{lll}{\widehat{\mathrm{R}}}_{\mathit{Y}}[k]& =& \frac{1}{N_0}{\sum }_{n=1}^{N_0}w[k-n]\mathit{Y}[n]{\mathit{Y}}^H[n]\\ & =& A\frac{1}{N_0}{\sum }_{n=1}^{N_0}w[k-n]\mathit{x}[n]{\mathit{x}}^H[n]A^H\\ & =& A{\widehat{\mathrm{R}}}_{\mathit{x}}[k]A^H,\end{array}$$

(24)

where ${\widehat{\mathrm{R}}}_{\mathit{x}}[k]=\frac{1}{N_0}{\sum }_{n=1}^{N_0}w[k-n]\mathit{x}[n]{\mathit{x}}^H[n]$. Define the diagonal elements of ${\widehat{\mathrm{R}}}_{\mathit{x}}[k]$ as ${\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[k]=\mathrm{diag}({\widehat{\mathrm{R}}}_{\mathit{x}}[k])$. Similar to Equation (19), we can estimate the power spectrum of signal using least square algorithm

$${\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[k]={\Phi }^{†}{\widehat{\mathrm{r}}}_{\mathit{Y}}[k], k=1,2,\dots ,N_0,$$

(25)

where ${\widehat{\mathrm{r}}}_{\mathit{Y}}[k]=\mathrm{Vec}({\widehat{\mathrm{R}}}_{\mathit{Y}}[k])$.

Finally, the spectrum-sensing problem is subsequently regarded as a binary hypothesis testing problem,

$$\begin{array}{ll}{\mathrm{H}}_0:& {\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[k]~\mathcal{N}\left({\sigma }^2,{\rho }^2\right)\\ {\mathrm{H}}_1:& {\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[k]~\mathcal{N}\left({\sigma }^2+S[k],{\rho }^2\right).\end{array}$$

(26)

The noise level ${\sigma }^2$ and estimate variance ${\rho }^2$ are assumed to be known, which can be obtained from receiver calibration, but the power spectrum of signal $S[k]$ is unknown. Next, constant false alarm probability ${\mathrm{P}}_{\mathrm{FA}}$ is imposed across all frequency bins, and the detection threshold is determined as

$$\gamma =\rho {\mathrm{Q}}^{-1}\left({\mathrm{P}}_{\mathrm{FA}}\right)+{\sigma }^2,$$

(27)

where ${\mathrm{Q}}^{-1}(·)$ is the tail distribution function of the standard normal distribution.

4.2. Efficient Spectrum-Sensing Algorithm

The traditional spectrum-sensing method [18] first estimates the power spectrum of signal, and then uses the estimated power spectrum for primary user detection. However, the computational load of power spectrum estimation is heavy. The power spectrum estimation requires repeatedly solving Equation (25) $N_0$ times. Therefore, the power spectrum estimation requires a large amount of computation when the number of samples of a coset $N_0$ is large. Interestingly, note that the frequency-smoothed cyclic cross-periodogram is used to estimate ${\widehat{\mathrm{R}}}_{\mathit{Y}}[k]$ in Equation (24), and ${\widehat{\mathrm{R}}}_{\mathit{x}}[k]$ is also a frequency-smoothed result of the original ${\mathrm{R}}_{\mathit{x}}[k]$. If the width of the spectral window function $w[n]$ is wide enough, the estimated power spectrum ${\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[k]=\mathrm{diag}({\widehat{\mathrm{R}}}_{\mathit{x}}[k])$ at frequency bin $k$ will contain the power spectrum information from the surrounding frequency bins $k+1,k-1,k+2,k-2,\dots$. In other words, the estimated power spectrum ${\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[k]$ will reflect the outline information of the original power spectrum. Since the power spectrum information is highly correlated at frequency bins in the spectral window, it is not necessary to estimate the power spectrum information at all frequency bins. Therefore, to reduce computational complexity, only power spectrum at partial frequency bins is estimated to conduct spectrum sensing. The proposed efficient spectrum-sensing algorithm exploits the frequency-smoothing property of the power spectrum estimator to reduce computational complexity, and it is summarized as Algorithm 1.

Algorithm 1. Efficient Spectrum-Sensing Algorithm based on Power Spectrum at Partial Frequency Bins

Initialization: Classify the periodic non-uniformly sampled data $y[n]$ to $q$ uniform sampling sequences $y_{ci}[k],1\le i\le q, 1\le k\le N_0$, and calculate their FFT $Y_{ci}[k], 1\le i\le q, 1\le k\le N_0$. Step 1: Implement fractional delay compensation in the frequency domain, $Y_i[k]=Y_{ci}[k]·P_i[k], 1\le i\le q, 1\le k\le N_0$, $P_i[k]={\mathrm{e}}^{-j2\pi c_{i}k/\left(L{N}_0\right)}$. Step 2: The width of the spectral window $w[n]$ is chose to be $W_{\min }/△f$. Compute the frequency-smoothed cross-periodogram of $Y_i[m]$ and $Y_j[m]$ ${\{{\widehat{R}}_{\mathit{Y}}[m]\}}_{i,j}=\frac{1}{N_0}{\sum }_{n=1}^{N_0}w[m-n]Y_i[n]Y_j^{*}[n]$ at partial frequency bins $m=l{N}_s,1\le l\le N_0/N_s$. Step 3: Estimate the power spectrum at partial frequency bins ${\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[m]={\Phi }^{†}{\widehat{\mathrm{r}}}_{\mathit{Y}}[m], m=l{N}_s,1\le l\le N_0/N_s$. Decision: Compare ${\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[m], m=l{N}_s,1\le l\le N_0/N_s$ with the predefined threshold $\gamma$, and a sequence of binary results can be obtained. Then the result can be obtained by flipping consecutive 1s with length less than $W_{\min }/△f$ to 0s.

In the initializing stage, $q$ uniform sampling sequences $y_{ci}[k],1\le i\le q, 1\le k\le N_0$ can be obtained by classifying the periodic non-uniformly sampled data $y[n], 1\le n\le N$. Then the FFT of $q$ uniform sampling sequences $y_{ci}[n]$ are calculated and expressed as $Y_{ci}[k],1\le i\le q, 1\le k\le N_0$.

The first step achieves the function of compensating the fractional delay of $q$ uniform sampling sequences $y_{ci}[n]$ in the frequency domain. Compared with the traditional method which uses the fractional delay filter to compensate the fractional delay, this method has lower computational complexity. Detailed analysis of the reduction of computational complexity is beyond the scope of this paper.

In the second step, let $N_s$ denotes the interval of choosing partial frequency bins from all bins uniformly. The frequency-smoothed cross-periodogram of $Y_i[m]$ and $Y_j[m]$ at partial frequency bins $m=l{N}_s,1\le l\le N_0/N_s$ is calculated. We assume that $N_0$ is an integer multiple of $N_s$. Define $W_{\min }$ as the minimum modulation bandwidth of primary user, and $△f$ as the bandwidth of a frequency bin. The width of the spectral window function $w[n]$ and $N_s$ are both chosen to be $W_{\min }/△f$ to obtain a smoothed version of the original power spectrum. Please note that the width of the spectral window function is very important to obtain the outline information of the original power spectrum.

The third step uses the least square method to estimate the power spectrum at partial frequency bins ${\widehat{\tilde{\mathrm{r}}}}_{\mathit{x}}[m]$, which only needs repeatedly solving Equation (25) $N_0/N_s$ times. Compared with the traditional method which needs repeatedly solving Equation (25) $N_0$ times, the amount of calculation has been reduced by $N_s$ times.

Finally, the estimated power spectrum at partial frequency bins ${\widehat{\tilde{\mathrm{r}}}}_x[m]$ are compared with the predefined threshold $\gamma$. A sequence of binary results can be obtained after the comparison. The minimum modulation width of primary user $W_{\min }$ and the width of the frequency bin $△f$ can be exploited to further refine the binary sequence. Consecutive 1s with length less than $W_{\min }/△f$ are flipped to 0s, then channels with 1s are considered occupied by primary users.

Some applications are sensitive to the computational complexity of the spectrum-sensing module, such as sensor nodes and battery-powered devices. Compared with the traditional method which estimates the power spectrum at all frequency bins, the proposed efficient spectrum-sensing algorithm only estimates power spectrum at partial frequency bins, by exploiting the frequency-smoothing characteristics of the power spectrum estimator. It is reasonable to use the power spectrum information at partial frequency bins to determine whether the primary user exists because the estimate result in the third step of Algorithm 1 is the weighted average result of the original power spectrum in a spectral window function.

5. Sampling Pattern Design

The matrix $\Phi$ is assumed to have full column rank when the least square method is used to estimate the power spectrum in Equation (25). In fact, whether it has full column rank is determined by the sampling pattern of MCS. Let us explore the structure of matrix $\Phi =\left(A^{*}\otimes A\right)B$. Let $a_i^{*}$ and $a_i,1\le i\le q$, denote the i-th row of matrix $A^{*}$ and $A$, respectively. $A^{*}$ and $A$ can be further written as $A^{*}={\left(a_1^H,a_2^H,\dots ,a_q^H\right)}^T$ and $A={\left(a_1^T,a_2^T,\dots ,a_q^T\right)}^T$, respectively. The structure of $q^2\times L^2$ matrix $A^{*}\otimes A$ is firstly explored and it can be written as:

$$\begin{array}{lll}{A}^{*}\otimes A& =& ({\left(a_1^{*}\otimes a_1\right)}^{\mathrm{T}},{\left(a_1^{*}\otimes a_2\right)}^{\mathrm{T}},\dots ,{\left(a_1^{*}\otimes a_q\right)}^{\mathrm{T}},\\ & & {\left(a_2^{*}\otimes a_1\right)}^{\mathrm{T}},{\left(a_2^{*}\otimes a_2\right)}^{\mathrm{T}},\dots ,{\left(a_2^{*}\otimes a_q\right)}^{\mathrm{T}},\\ & & ⋮\\ & & {\left(a_q^{*}\otimes a_1\right)}^{\mathrm{T}},{\left(a_q^{*}\otimes a_2\right)}^{\mathrm{T}},\dots ,{\left(a_q^{*}\otimes a_q\right)}^{\mathrm{T}}{)}^{\mathrm{T}},\end{array}$$

(28)

where $a_i^{*}\otimes a_k, 1\le i,k\le q$ is the $\left(i-1\right)q+k$-th row vector of matrix $A^{*}\otimes A$. The $L^2\times L$ selection matrix $B$ can be written as:

$$B={\left(\begin{array}{lll}\begin{array}{ll}1& 0\\ 0& 0\end{array}& \cdots & \begin{array}{ll}0& 0\\ 0& 0\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{ll}0& 0\\ 0& 0\end{array}& \cdots & \begin{array}{ll}0& 0\\ 0& 0\end{array}\end{array}\begin{array}{lll}\begin{array}{ll}0& 0\\ 0& 1\end{array}& \cdots & \begin{array}{ll}0& 0\\ 0& 0\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{ll}0& 0\\ 0& 0\end{array}& \cdots & \begin{array}{ll}0& 0\\ 0& 0\end{array}\end{array}\dots \begin{array}{lll}\begin{array}{ll}0& 0\\ 0& 0\end{array}& \cdots & \begin{array}{ll}0& 0\\ 0& 0\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{ll}0& 0\\ 0& 0\end{array}& \cdots & \begin{array}{ll}0& 0\\ 0& 1\end{array}\end{array}\right)}^T.$$

(29)

Based on the structure of matrix $B$, it is not difficult to find that multiplying a $1\times L^2$ vector $a_i^{*}\otimes a_k$ with a $L^2\times L$ selection matrix $B$ is equivalent to $a_i^{*}.*a_k$, where $.\ast$ denotes the matrix dot product. Therefore, matrix $\Phi =\left(A^{*}\otimes A\right)B$ can be rewritten as:

$$\begin{array}{lll}\left(A^{*}\otimes A\right)B& =& ({\left(a_1^{*}.\ast a_1\right)}^{\mathrm{T}},{\left(a_1^{*}.\ast a_2\right)}^{\mathrm{T}},\dots ,{\left(a_1^{*}.\ast a_q\right)}^{\mathrm{T}},\\ & & {\left(a_2^{*}.\ast a_1\right)}^{\mathrm{T}},{\left(a_2^{*}.\ast a_2\right)}^{\mathrm{T}},\dots ,{\left(a_2^{*}.\ast a_q\right)}^{\mathrm{T}},\\ & & ⋮\\ & & {\left(a_q^{*}.\ast a_1\right)}^{\mathrm{T}},{\left(a_q^{*}.\ast a_2\right)}^{\mathrm{T}},\dots ,{\left(a_q^{*}.\ast a_q\right)}^{\mathrm{T}}{)}^{\mathrm{T}}.\end{array}$$

(30)

Since the row vector $a_i$ of matrix $A$ can be written as $a_i=[\exp \left(-j2\pi c_i\frac{1}{L}\right),\exp \left(-j2\pi c_i\frac{2}{L}\right),\dots ,\exp \left(-j2\pi c_i\frac{L}{L}\right)]$, the $\left(i-1\right)q+k$-th row vector of matrix $\Phi =\left(A^{*}\otimes A\right)B$ can be further written as

$$a_i^{*}.\ast a_k=[\exp \left(-j2\pi \left(c_k-c_i\right)\frac{1}{L}\right),\exp \left(-j2\pi \left(c_k-c_i\right)\frac{2}{L}\right),\dots ,\exp \left(-j2\pi \left(c_k-c_i\right)\frac{L}{L}\right)].$$

(31)

Please note that $a_i^{*}.\ast a_k$, the $\left(i-1\right)q+k$-th row of matrix $\Phi$, has some special structure. It can be viewed as the $c_k-c_i+1$-th row of the $L\times L$ Discrete Fourier Transform (DFT) matrix when $c_k-c_i$ is greater than or equal to 0.

Matrix $\Phi =\left(A^{*}\otimes A\right)B$ can be written as the product of two matrices using Equations (30) and (31),

$$\Phi =E\times F,$$

(32)

where $F$ is a $L\times L$ DFT matrix and $E$ is a $q^2\times L$ matrix. The row vector of matrix $E$ contains only one non-zero element 1 and the other elements are zero. The location of element 1 is $c_k-c_i+1$ when $c_k-c_i$ is greater than or equal to 0. The location of element 1 is $c_k-c_i+L$ when $c_k-c_i$ is less than 0. The row vector of matrix $E$ acts similar to a selection vector, and it selects which row vector of the DFT matrix $F$ as the current row vector of $\Phi$. Obviously, the matrix $\Phi$ has full column rank as long as the matrix $E$ has full column rank since the DFT matrix $F$ is a full rank square matrix.

Next, we briefly introduce the concept of the circular sparse ruler. A L-length circular sparse ruler can be viewed as a ruler with only $q$ marks which is less than L, but it is also able to measure all the length between 1 and L. Define the $q$ marks as $C=\left\{c_1,c_2,\dots ,c_q\right\}, C\subset \left\{1,2,\dots ,L\right\}$, and their difference set $\mathsf{\Omega }\left(C\right)=\left\{\left(c_k-c_i\right)\mathrm{mod} L|\forall c_k,c_i\in C\right\}$, where $\left(·\mathrm{mod}·\right)$ is the modulus operator. A set $C$ is a circular sparse ruler if its difference set $\mathsf{\Omega }\left(C\right)=\left\{1,2,\dots ,L\right\}$. An example of the circular sparse ruler is given in Figure 4, where the marks are $C=\left\{1,2,9,12,14,18\right\}$ and the measurement range is $L=24$. Although there is no mark 7 and 17 on the sparse ruler, length 7 and 17 can be measured using mark 2 and 9.

The definition of difference set of circular sparse ruler $\mathsf{\Omega }\left(C\right)=\left\{\left(c_k-c_i\right)\mathrm{mod} L|\forall c_k,c_i\in C\right\}$ can be rewritten as $\mathsf{\Omega }\left(C\right)={\mathsf{\Omega }}_1\left(C\right)\cup {\mathsf{\Omega }}_2\left(C\right)$, where ${\mathsf{\Omega }}_1\left(C\right)=\left\{c_k-c_i+1|\forall c_k,c_i\in C \mathrm{and} c_k-c_i\ge 0\right\}$ and ${\mathsf{\Omega }}_2\left(C\right)=\left\{c_k-c_i+L|\forall c_k,c_i\in C \mathrm{and} c_k-c_i<0\right\}$. That is, if the difference set of the sampling pattern of MCS includes the set $\left\{1,2,\dots ,L\right\}$, then each column vector of matrix $E$ contains at least one non-zero element 1. Combined with the fact that the row vector of $E$ has only one non-zero element 1, the matrix $E$ will have full column rank. Therefore, the sampling pattern design of the MCS can be viewed as a circular sparse ruler problem. Please note that circular sparse ruler uses the difference between two marks to measure any length between 1 and $L$. If we want to measure length 1 to length $L$, then at least one pair of marks, or cosets, should have length 1 difference, or $1\times T$ time difference. Further, to obtain the strongest compression rate, we would like to minimize the cardinality of $C$, which can be viewed as a minimal circular sparse ruler problem. Then finally, let $c_{i+2}$ and $c_i$ denote the i+2-th and i-th coset, respectively. $c_{i+2}$ should be greater than $c_i+d$ because the S/H in Figure 2 can only temporarily store one sample and the minimal sampling interval of ADC is $T_{\mathrm{ADC}}=d\times T$. Therefore, the design of sampling pattern turns into the following optimization problem:

$$\underset{C}{\min }|C| \mathrm{s}.\mathrm{t}. \{\begin{array}{l}\mathsf{\Omega }\left(C\right)=\left\{1,2,\dots ,L\right\}\hfill \\ c_{i+2}>c_i+d, 1\le i\le q-2\hfill \end{array},$$

(33)

where the sampling pattern $C\subset \left\{1,2,\dots ,L\right\}$.

6. Simulation Results

This section evaluates the spectrum-sensing performance of the proposed single-channel sub-Nyquist sampling by performing several numerical simulations. In practice, an antenna is used to receive the Radio Frequency (RF) signal, and a band-pass filter is utilized to filter the out-of-band noise and interference. Subsequently, the RF signal is demodulated to baseband signal using a mixer. Finally, the Local Oscillator (LO) signal mixed in the demodulated signal is filtered out using a low-pass filter. Since the receiver front-end circuit is not the focus of this paper, we assume that the proposed single-channel sub-Nyquist sampling directly deals with the demodulated signal output from the low-pass filter in Figure 2. Quadrature Phase-Shift Keying (QPSK) signal is used as the test signal, and it is generated by the following model:

$$\begin{array}{lll}x\left(t\right)& =& s\left(t\right)+n\left(t\right)\\ & =& E{\sum }_{k=1}^M\{{\sum }_{i=1}^{N_b}{I}_k[i]g_k\left(t-i{T}_k\right)\cos \left(2\pi f_{k}t\right)+{\sum }_{j=1}^{N_b}{Q}_k[j]g_k\left(t-j{T}_k\right)\sin \left(2\pi f_{k}t\right)\}+n\left(t\right),\end{array}$$

(34)

where $M$ is the number of primary users, $N_b$ is the number of random bits, $I_k[i]$ and $Q_k[j]$ are random bit streams, $g_k\left(t\right)$ is the pulse-shaping function, $T_k$ is the symbol duration, $f_k$ is the carrier frequency of the k-th primary user, and $n\left(t\right)$ is the additive white Gaussian noise. Simulation parameters are set as follows:

• The bandwidth of primary user is set as $\frac{1}{T_k}=8 \mathrm{MHz}$. The carrier frequency of primary user is randomly chosen from the center frequency of 40 channels. The pulse-shaping function is the root-raised cosine with roll-over factor 0.1. $I_k[i]$ and $Q_k[j]$ are generated from pseudo-random binary sequence generator.
• Additive white Gaussian noise $n\left(t\right)~\mathcal{N}\left(0,{\sigma }^2\right)$, where ${\sigma }^2$ is fixed to 1 and $E$ in Equation (34) is scaled to a certain SNR level.
• The minimum modulation bandwidth of primary user is set as $W_{\min }=500\mathrm{KHz}$.
• To simulate the mismatch behavior of parallel MCS, the clock timing skew is set as 2% of the Nyquist sampling interval $T$. The offset and gain mismatch are set as 2% and 2% of maximum amplitude of signal, respectively.
• The parameters of log-normal shadowing in Equation (2) are set as follows: reference distance $d_0=1\mathrm{km}$, the path loss at the reference point $P{L}_0=0 \mathrm{dB}$, the path loss exponent $\gamma =2.4$, the Gaussian random variable $X_{\rho }~\mathrm{N}\left(0,5\right) \mathrm{dB}$, the distance between the primary user and the cognitive radio device $d~\mathrm{N}\left(5,3\right) \mathrm{km}$. In addition, primary users are passed through Rayleigh channels with $2{\lambda }^2=2 \mathrm{MHz}$.
• In the second and third numerical simulation, 200 trials are performed, and averaged results are shown.

In the first simulation, we test the feasibility of spectrum sensing using power spectrum at partial frequency bins. To show the estimated frequency-smoothed power spectrum clearly, we only observed 12 of these 40 channels. In the following two simulations, 32 of these 40 channels are observed. The parameters of MCS are set as $L=24$, $q=6$ and the sampling pattern is $C=\left\{1,2,9,12,14,18\right\}$. The number of active primary users is set as $M=3$. The SNR of the test signal is set as 10 dB, and the acquisition time is set as 64 μS. Figure 5a shows the original power spectrum of the test signal. Figure 5b,c show the estimated frequency-smoothed power spectrum at partial frequency bins using the proposed serial MCS and the traditional parallel MCS, respectively.

Figure 5b reflects the profile information of the original power spectrum because frequency-smoothed cross-periodogram is used in Equation (24). Using the frequency-smoothed cross-periodogram is equivalent to taking a weighted average of the power spectrum in the frequency-smoothing window. The weight value in the middle of the window is the largest and decreases toward each side. The profile information of the original power spectrum is sufficient for the spectrum-sensing module to detect primary users. It is unnecessary to estimate power spectrum at all frequency bins because the power spectrum at two frequency bins within the frequency-smoothing window are highly correlated. The proposed efficient spectrum-sensing algorithm exploits the frequency-smoothing characteristic to reduce the computational complexity.

On the other hand, as shown in Figure 5c, mismatches among sub-ADCs in traditional parallel MCS cause large errors in the estimated power spectrum. Although these errors do not exist after calibration, the calibration will add a lot of additional burden to the system. Therefore, the proposed serial MCS combined with the efficient spectrum-sensing algorithm is a better choice.

In practice, the primary user signal arrives at the antenna of the cognitive radio device through multiple paths. Due to the different arrival time, the superposition of these signals will result in the distortion of the original signal. Rayleigh distribution is a model that can be used to describe this form of fading when multipath propagation exists. Rayleigh fading will result in Doppler shift, which equivalents to the expansion in the frequency domain. As shown in the green circle of Figure 5, although the baseband width of the primary user of the third channel is 8 MHz, the Doppler shift caused by the multipath effect causes the energy of the third channel leak into the second channel. As a result, cognitive radio devices may mistakenly regard the primary user of the second channel also exists. In addition, Log-normal Shadowing significantly attenuates the amplitude of signal. When the amplitude of primary user signal is close to that of noise, the proposed method is difficult to detect the presence of the primary user. Cooperative spectrum sensing based on decision fusion can efficiently solve the shadow fading problem.

In the second simulation, we evaluate spectrum-sensing performance in terms of different compression ratio. 32 channels are observed in this simulation. Its bandwidth and Nyquist sampling rate are 256 MHz and 512 MSPS (Mega-Samples Per Second), respectively. We consider several SNRs and they are 5 dB, 0 dB, −5 dB and −10 dB. The periodicity of the sampling instants is set as $L=64$ in MCS. Several compression ratio are considered, and they are 0.25, 0.375 and 0.5. The sampling pattern for compression ratio 0.25 is $C=\left\{1,5,9,10,14,18,25,26,33,37,38,43,49,51,54,61\right\}$, and the sampling rate requirement for ADC is 128 MSPS. The sampling pattern for compression ratio 0.375 is $C=\left\{1,3,5,7,9,10,12,14,16,18,21,25,26,29,33,37,38,41,43,45,49,51,54,61\right\}$, and the sampling rate requirement for ADC is 192 MSPS. The sampling pattern for compression ratio 0.5 is $C=\left\{1,3,5,7,9,10,12,14,16,18,19,21,25,26,27,29,33,35,37,38,39,41,42,43,45,47,49,51,53,54,61,62\right\}$, and the sampling rate requirement for ADC is 256 MSPS. The number of active primary users is set as $M=8$. The acquisition time is set as 64 μS.

Figure 6 shows the Receiver Operating Characteristic (ROC) curve of the detector for different compression ratios. It is noticed that spectrum-sensing performance under different compression ratios is different. The more compression, the worse the detection performance. The curve corresponding to ‘Compression Ratio = 0.375 (Mismatch)’ in Figure 6 is spectrum-sensing performance based on the parallel MCS with mismatches. Please note that the detection performance of the parallel MCS is better than that of the proposed serial MCS. The reason is that the mismatch error component is mistakenly regarded as primary users. The curve corresponding to ‘Compression Ratio = 0.375 (PSD)’ in Figure 6 is spectrum-sensing performance using the power spectrum at all frequency bins. The detection performance is slightly different from the detection performance using the power spectrum at partial frequency bins, while the later one reduces a large amount of computation.

Table 1. Computational Complexity Required.

Compression Ratio	Conventional	Proposed	Saved
0.25	512 1	32	93.75%
0.33	512	32	93.75%
0.5	512	32	93.75%

¹ The unit is the number of solving the Least Square Equation.

shows the reduced computational complexity. The proposed efficient spectrum-sensing algorithm uses the frequency-smoothed cyclic cross-periodogram to estimate the smoothed power spectrum. Since the power spectrum of the smoothed power spectrum at two adjacent frequency bins are highly correlated, it is not necessary to estimate the power spectrum at all frequency bins, which reduces the computational complexity of spectrum sensing.

In the third simulation, we evaluate spectrum-sensing performance in terms of different acquisition time. The proposed efficient spectrum-sensing algorithm requires the signal be stationary during the acquisition process. To consider the fast time-varying signal, the detection performance under short acquisition time condition is tested. The compression ratio is fixed to 0.375. Several acquisition time are considered, and they are 16 μS, 64 μS, and 128 μS. Other simulation parameters are the same with the second simulation.

As shown in Figure 7, spectrum-sensing performance gradually improves as the SNR increases. When the SNR is greater than 5 dB, the proposed efficient spectrum-sensing algorithm exhibits a perfect detection performance. The shorter the acquisition time, the worse the spectrum-sensing performance. Short acquisition time means that the number of samples used to calculate frequency-smoothed cross-periodogram in Equation (24) is small. Thus, the estimate accuracy of the frequency-smoothed cross-periodogram is affected. Ultimately, spectrum-sensing performance is influenced. In the case that the acquisition time is 16 μS, it is guaranteed that all primary users and spectrum holes can be detected with high probability when the SNR is higher than 5 dB.

The curve corresponding to ‘Acquisition Time = 64 μS (PSD)’ is spectrum-sensing performance using power spectrum at all frequency bins. It is shown that spectrum-sensing performance based on power spectrum at all frequency bins has little improvement. Therefore, spectrum sensing using power spectrum at partial frequency bins is preferable since it reduces a large amount of computation.

Table 2. Computational Complexity Required.

Acquisition Time	Conventional	Proposed	Saved
16 μS	128 1	8	93.75%
64 μS	512	32	93.75%
128 μS	1024	64	93.75%

¹ The unit is the number of solving the Least Square Equation.

compares the required computational complexity between the conventional method and the proposed method. The reason for the reduction of computational complexity is the same as that of the second simulation.

Compared with the detection performance, the false alarm probability is also very important. False alarm probability determines how many spectrum holes are available to secondary users. The curve corresponding to ‘Acquisition Time = 64 μS (Mismatch)’ is spectrum-sensing performance using the traditional parallel MCS with mismatches. Mismatches have little influence on the detection performance as shown in Figure 7a, while they have severe influence on the false alarm probability as shown in Figure 7b. The false alarm probability of the proposed serial MCS gradually decreases to 0 as the SNR increases. However, the false alarm probability of the traditional parallel MCS does not tend to zero even if the SNR is high. As shown in Figure 7b, all spectral holes can be detected with high probability when the SNR is higher than 8dB for serial MCS.

7. Discussion

The periodic non-uniformly sampled data obtained by the proposed serial MCS can be classified into several uniform sampling sequences. Waveform reconstruction can be achieved using MCS algorithm and compressive sensing [8]. The reconstructed waveform can be used for decoding, so the spectrum-sensing module can also be used as a receiver. The alternative function of the spectrum-sensing module saves hardware cost in terms of price and size. It is worth noting that the sampling pattern design for waveform reconstruction is different from the sampling pattern design for power spectrum estimation proposed in this paper, so the sampling pattern design should be given additional consideration.

8. Conclusions

This paper proposes a single-channel sub-Nyquist sampling structure, serial MCS, for wideband spectrum sensing. Although the requirement for the sampling rate of ADC is higher compared with the traditional parallel MCS, the proposed serial MCS cannot only avoid the mismatch among sub-ADCs in traditional parallel MCS, but also reduce the hardware implementation size. An efficient spectrum-sensing algorithm based on power spectrum at partial frequency bins is also developed, which could save a large amount of computational cost. A sampling pattern design method for the proposed serial MCS is given as well. Simulations show that the mismatch in traditional MCS has a serious impact on spectrum-sensing performance, and the proposed MCS combined with the efficient spectrum-sensing algorithm exhibits outstanding spectrum-sensing performance at much lower computational cost.