Compressive cyclostationary spectrum sensing with a constant false alarm rate

Spectrum sensing is a crucial component of opportunistic spectrum access schemes, which aim at improving spectrum utilization by allowing for the reuse of idle licensed spectrum. Sensing a spectral band before using it makes sure the legitimate users are not disturbed. To that end, a number of different spectrum sensing method have been developed in the literature. Cyclostationary detection is a particular sensing approach that takes use of the built-in periodicities characteristic to most man-made signals. It offers a compromise between achievable performance and the amount of prior information needed. However, it often requires a significant amount of data in order to provide a reliable estimate of the cyclic autocorrelation (CA) function. In this work, we take advantage of the inherent sparsity of the cyclic spectrum in order to estimate CA from a low number of linear measurements and enable blind cyclostationary spectrum sensing. Particularly, we propose two compressive spectrum sensing algorithms that exploit further prior information on the CA structure. In the first one, we make use of the joint sparsity of the CA vectors with regard to the time delay, while in the second one, we introduce structure dictionary to enhance the reconstruction performance. Furthermore, we extend a statistical test for cyclostationarity to accommodate sparse cyclic spectra. Our numerical results demonstrate that the new methods achieve a near constant false alarm rate behavior in contrast to earlier approaches from the literature.


Introduction
The scarcity of radio spectrum constitutes a major roadblock to current and future innovation in wireless communications.To alleviate this problem, it has been proposed to make spectral resources, which are currently underutilized, available for reuse under a paradigm that goes by the name of opportunistic spectrum access (OSA) [1].Spectrum sensing is one of its core technologies.It allows an unlicensed transceiver, a so-called secondary user (SU), to access a licensed spectral band without interfering with the owner of the band's license, the so-called primary user (PU).The fundamental task in spectrum sensing is to decide between two hypotheses, the first of which states that the spectral band under investigation is free (H 0 ), while the second asserts that it is occupied (H 1 ).Considering the baseband signal x(t) observed at a secondary system receiver, the two hypotheses can be written as where η(t) denotes receiver noise and s (t) stands for a PU signal after propagation effects.
A number of spectrum sensing algorithms have been proposed in the literature [2][3][4].Broadly speaking, they can be divided into three major types, namely, energy detection, stochastic feature detection, and matched filter detection, where different types require different amounts of prior knowledge about the PU signal.While matched filter ( [5], Ch. 4.3) detectors require the knowledge of the exact waveform of at least a part of the PU signal, energy detection [6] does not require any prior knowledge.Feature detectors are an in-between as they only make assumptions about structural or statistical properties of the signal.One of the stochastic features which lets an SU receiver discriminate between pure stationary noise (H 0 ) and a communication signal contaminated with noise (H 1 ) is cyclostationarity.In contrast to pure stationary noise, most man-made signals vary periodically with time [7] and can thus be characterized as cyclostationary.Although the data contained in a modulated signal may be a purely stationary random process, the coupling with sine wave carriers, pulse trains, repeating, spreading, hopping sequences, and cyclic prefixes going along with its modulation causes a built-in periodicity [8].
The use of cyclostationarity for the purpose of spectrum sensing has been investigated from a variety of perspectives ranging from single node signal detection [9][10][11][12][13][14][15][16] to collaborative approaches that take use of the spatial diversity [17,18].One of the particularly well-known algorithms for cyclostationary spectrum sensing is the so-called time-domain test (TDT) as introduced in [9].The test can decide between the presence and absence of cyclostationarity for a pre-specified potential cycle frequency α.It operates on the cyclic autocorrelation (CA), which, given an observed signal x(t), is defined as [ for a potential cycle frequency α and a delay τ .For purely stationary signals, R α x (τ ) = 0 for all α = 0, while for cyclostationary signals, R α x (τ ) = 0 for some α = 0.The α with non-zero CA coefficients are called cycle frequencies.The set of cycle frequencies caused by one of potentially multiple incommensurate second-order periodicities in a cyclostationary signal comprises the periodicity's fundamental cycle frequency (the reciprocal of the fundamental period) as well as its harmonics.Given the above information, we can rewrite the hypothesis test (1) as It is important to note that practically, instead of the statistical CA (2), one normally operates on the sample CA obtained from a limited number of signal samples.The coefficients of the sample CA are not constant but rather follow different probability distributions, depending on whether H 0 or H 1 is true.To account for this, the hypothesis test (3) is modified by considering different test statistic under H 0 and H 1 .As a result, the TDT provides a typical constant false alarm rate (CFAR) performance.
Returning to (3), we note that the CA is zero on its whole support except the set of cycle frequencies and α = 0. Therefore, it can be called sparse.The exploitation of sparsity in signal processing has a long history [19].The recent years, however, have seen a vastly accelerated development of the field resulting in a new sampling paradigm called compressive sampling (CS) [20,21].It postulates that sparse or compressible signals, i.e., the signals that can be represented or well approximated by only a few non-zero coefficients in some domain, can be sampled and recovered from fewer number of measurements than traditionally required.A crucial observation here is that one can design an effective measurement strategy that is governed by the amount of the signals' information content, rather than its ambient dimension.To date, there is a large amount of powerful algorithms available that solve sparse recovery problems in the CS context ranging from optimization approaches to classical pursuits such as the orthogonal matching pursuit (OMP) [22] and more specialized algorithms such as the compressive sampling matching pursuit (CoSaMP) [23] for instance.
Multiple contributions have been made in the field of compressive cyclostationary spectrum sensing.The authors of [24] for instance formulate the estimation of the cyclic autocorrelation as a sparse recovery problem, which they solve using the OMP [22].Based on the sparse estimate of the CA, they propose two detection methods that exploit different CA properties.The first one, called slot comparison method (SCM), compares the biggest CA components OMP finds in two consecutive blocks of samples.If for both blocks, the same discrete cycle frequencies are chosen, H 1 is selected; otherwise, H 0 is selected.The second detection algorithm is called symmetry method (SM).It exploits the fact, that for certain types of signals, the CA is symmetric around the direct current (DC) component.Although both of them present blind detectors meaning that they operate without prior knowledge of the cycle frequencies present in the signal, they do not allow for a CFAR performance, which is considered a desired detector feature [4].Instead of the CA, the authors of [25] use the spectral correlation (SC), which is the Fourier transform of the CA over τ , for detecting multiple transmitters in a wideband signal using compressive sampling.In order to estimate the SC from compressed samples via CS, they established a direct linear relation between the compressed samples and the SC.Based on [25], the authors of [26] derive a method for recovering the SC from sub-Nyquist samples using a reduced complexity approach, for which they provide a closed-form solution.In [27], the modulated wideband converter (MWC) [28] is used to obtain the SC from sub-Nyquist samples to then apply cyclostationarity detection.Furthermore, it has been recently shown that under certain conditions, the CA can be efficiently recovered from a low number of samples even without enforcing the sparsity property [29][30][31].This can be done by exploiting cross-correlations between different outputs of the compressive sampler.The main drawback of the aforementioned works is that providing the estimate of the entire cyclic spectrum, they still require the knowledge of the cycle frequencies for the detection step.
In this work, we employ a composite approach that combines the sparse recovery of CA from its compressive measurements for blind cycle frequency estimation with a CFAR TDT detection.This said the contribution of this paper is manifold.We propose two novel sparsity-aided CA estimation algorithms, both of which exploit further prior information about the CA in addition to its sparsity: the simultaneous OMP-based (SOber) and the dictionary assisted (Dice) compressive CA estimator.The first one exploits the joint sparsity of the CA vectors with regard to the time delay in order to recover the CA matrix for all delays simultaneously, while the second one takes advantage of the signal-induced structure of the CA by introducing structure dictionaries into the recovery process.In order to evaluate the performance of the proposed CA estimators, we derive a closed-form expression of the CA of sampled linearly modulated signals with rectangular pulse shape.Furthermore, we show how this expression can be used as prior information in the dictionary assisted approach.Note that the use of sparse recovery in the novel CA estimation approaches results in the automatic detection of signal's cycle frequencies.This in turn allows blind spectrum sensing by eliminating the integral need of the classical TDT for the perfect knowledge of the said cycle frequencies.However, the resulting sparse structure of the compressive CA estimates does not allow for the application of the traditional TDT since the noise statistics are missing.To compensate for this phenomenon, we develop a modified TDT and thus enable blind compressive cyclostationary spectrum sensing.Numerical tests show that the proposed method achieves a near CFAR behavior.
The remainder of this paper is structured as follows.Section 2.1 introduces the signal model and presents the classical method for CA estimation, while Section 2.2 presents the time-domain test based on the classical CA estimation.A CA estimator based on joint sparsity of multiple vectors is introduced in Section 3 and the CA estimator exploiting additional prior knowledge is described in Section 4.An extension of the TDT to accommodate sparse CA estimates is developed in Section 5.The numerical evaluation of the proposed estimation and detection approaches as well as the interpretation of the results is given in Section 6. Section 7 concludes the paper.

System model and CA estimation
Consider a secondary system receiver that needs to decide whether a certain spectral band is occupied or free.It samples the baseband signal x(t) uniformly with a sampling period T e .This results in the vector of discrete samples x t 0 ∈ C N , where We assume that the vector x t 0 is discrete and zero mean, and due to the nature of man-made signals, it represents an (almost [32], Ch. 1.3) cyclostationary process [9].The presence of stochastic periodicity in the samples and thus the presence of a man-made signal can be revealed by applying a detection algorithm such as the TDT to the CA of the samples.There are different ways of obtaining the CA from the baseband samples, one of which is the following (classical) estimator Evaluating this function results in the CA coefficient for the cycle frequency α = a NT e and the time delay τ = νT e , where a stands for the discrete cycle frequency and ν denotes the discrete time delay.Note that the factor e −jπ a N ν remains constant throughout the sum.It is a phase shift necessary to maintain compatibility with the symmetric CA (2).The estimator ( 5) is biased but exhibits a smaller estimation variance than an unbiased one [9].
We define rν x as an N length CA vector whose nth element is Rn x,t 0 (ν), i.e., Subsequently, we re-write the estimation of the CA vector as a matrix-vector product.To do so, we need the (N element) delay product with time delay τ = νT e , which is given by where • denotes component-wise multiplication.Note that since the receiver only takes N samples, x * t 0 +νT e is zero-padded at the end while y ν N is a vector of length N. The CA vector is now given by where F denotes the (N × N) discrete Fourier transform (DFT) matrix.The N × n ν CA matrix for time delays ν 1 T e , . . ., ν n ν T e is given by with

The time-domain test (TDT) for cyclostationarity
Given the statistical CA, one could decide between H 0 and H 1 by testing it for being non-zero at the signal's inherent cycle frequencies according to (3).However, as mentioned in Section 1, instead of the statistical CA, we only have access to its estimation, the sample CA (which asymptotically converges to the statistical CA).This hinders the direct applicability of (3) for signal detection as coefficients of the sample CA are not constant anymore.
In the seminal work [9], the probability distributions that the sample CA coefficients follow under H 0 or H 1 have been identified and a test for cyclostationarity based on this knowledge has been designed.The test is briefly described in the following.
Consider the 2n ν × 1 vector which represents the concatenation of the real and the imaginary part of the row of Rx corresponding to the discrete cycle frequency a 0 .The frequency a 0 is the cycle frequency of interest, i. e., the one for the presence of which we want to test the signal.Given this vector, we can formulate the following non-asymptotic hypotheses where r xx * (a 0 ) is the deterministic but unknown asymptotic counterpart of rxx * (a 0 ) and xx * (a 0 ) is the estimation error.Note that in contrast to the hypotheses from Eq. ( 3), this formulation considers the presence of cyclostationarity in the received signal for one fixed cycle frequency a 0 .Since r xx * (a 0 ) is non-random, the distribution of rxx * (a 0 ) under H 0 and H 1 only differs in mean.As shown in [9], the estimation error xx * (a 0 ) asymptotically follows a Gaussian distribution, i. e., lim where xx * (a 0 ) is the statistical covariance matrix of rxx * (a 0 ) and D = denotes convergence in distribution.The 2n ν × 2n ν covariance matrix xx * (a 0 ) can be computed as [9] xx where the (m, n)th entries of the n ν × n ν matrices Q and Q * are given by respectively.The term S y νm N y νn N (•, •) denotes the unconjugated, while the term S * y νm N y νn N (•, •) denotes the conjugated cyclic spectrum of a signal.One way to estimate these is to determine the following frequency-smoothed periodograms: where W is a normalized spectral window of odd length L. Looking at the Eqs.( 16) and ( 17), it becomes clear why the cyclic spectrum is often referred to as the spectral correlation.Given the estimated quantities described above, the following generalized likelihood ratio (GLR) test statistic can be derived [17] The test statistic can be interpreted as a normalized energy.The inverse of the covariance matrix scales rxx * (a 0 ) such that under H 0 , its entries follow a standard normal distribution.Thus, under H 0 , the test statistic asymptotically follows a central chi-squared distribution with 2n τ degrees of freedom i. e., lim = χ 2 2n ν , while under H 1 , the test statistic asymptotically follows a non-central chi-squared distribution with unknown noncentrality parameter λ, i. e., lim = χ 2 2n ν (λ).Based on the above test statistic, we can design a CFAR detector with some false alarm rate P fa by finding the corresponding decision threshold in the χ 2 2n ν tables.We cannot design a test based on a desired detection rate P d since although r xx * (a 0 ) is deterministic, it depends on the type of signal emitted by the transmitter as well as the signal to noise ratio (SNR) at the receiver, both of which are assumed to be unknown.
The classical approach for cyclostationary spectrum sensing is to apply the TDT to the CA estimate from (5).However, to do so, one needs to know which cycle frequency to test beforehand, which eliminates the possibility of true blind spectrum sensing.One could sequentially test the received signal for all possible cycle frequencies.However, with high probability, the estimation noise at some cycle frequency would have a value above the decision threshold, leading to a false alarm.

Sparsity-aided CA estimation: simultaneous OMP-based estimator
As discussed in Section 1, for most man-made signals the CA is (asymptotically) sparsely occupied, containing spikes only at the DC component as well as the cycle frequencies of inherent signal periodicities and their harmonics.In this section, we take advantage of this inherent sparsity and cast the CA estimation as a joint sparse recovery problem.Since this method is able to detect the CA's support, it removes the traditional approach's requirement of knowing the cycle frequencies beforehand, enabling thus blind cyclostationarity-based spectrum sensing possible.We begin by rewriting Eq. ( 9) as where where Y m contains a selection of m coefficients of the delay products for different delays ν i T e .Note that ν i ∈ [0, N − 1] are chosen such that xt 0 + ν i T e is non-empty.We now want to recover Rx from Y m by solving the underdetermined inverse problem (20).To do so, we exploit our knowledge about the CA's sparsity.
The straightforward solution would be to solve the following optimization problem where • 0 denotes the 0 "norm" [20], which is the number of non-zero entries in a vector, and vec{•} stands for the vectorization of a matrix, i. e., the concatenation of its columns to a single vector.Eq.21 is known to be a nonconvex combinatorial problem [20].One way to solve it within a practically feasible amount of time is to substitute the 0 -"norm" by its tightest convex relaxation, the 1 norm.With high probability, this produces the same result since for most large underdetermined systems of linear equations, the minimal 1 -norm solution is also the sparsest solution [33].Another way of solving (21) efficiently is by applying one of the many greedy sparse recovery algorithms that have been developed in the field of CS, such as, e.g., orthogonal matching pursuit (OMP) [22].
OMP is a greedy algorithm that iteratively determines a vector's support from an underdetermined system of linear equations and subsequently recovers the vector by solving a least-square problem.Using it, we could solve (21) for each column of Rx individually (as in [24]), i. e., we could solve min rν for each ν.Instead, we notice that the vectors rν are jointly sparse with regard to the time delay.Therefore, stacking rν x in Rx results in a row-sparse matrix whose rows are (asymptotically) non-zero only at the indices corresponding to the cycle frequencies.In order to exploit this additional structure, we propose to use an extension of OMP called simultaneous orthogonal matching pursuit (SOMP) [34] to recover the CA matrix Rx at once.The CA estimation based on SOMP is summarized in Algorithm 1, which we further refer to as SOber.
The goal of SOber is to find the indices of the atoms contained in Y m , i. e., the support of the columns of Rx , and subsequently recover the identified non-zero rows of Rx by solving least-square problems.We start with an empty support S 0 .Each iteration, one atom (a column of the matrix A = NMF −1 ) index is added to the support.The index is selected according to the sum of the absolute correlation values between the corresponding atoms and the delay products of different time delays (lines 3-4).Using the new support set S i , a least-square problem is solved for each column in Rx (lines 5-6).In each iteration, the atom index to be added to the index set is chosen according to the correlation between the residuum of Y m and the atom set.Since every iteration adds one index to the support set, one usually chooses n iter greater than or equal to the sparsity of the signal to be recovered.The difference between OMP (used in, e. g., [24]) and SOMP can be found in line 4, where SOMP jointly considers the amount of correlation between atoms and the delay products of multiple delays, while OMP would select the support of rν

Sparsity-aided CA estimation: dictionary-assisted estimator
In Section 3, we have described a SOMP-based algorithm that estimates the cycle frequencies and the CA from fewer samples than required using the classic approach by taking into consideration the inherent sparsity of the CA.In this section, we develop an algorithm that makes use of additional prior knowledge about the signal's structure in the form of structure dictionaries to further enhance the cycle frequency and CA estimation.Like SOber, the new algorithm does not require the prior knowledge about the cycle frequencies contained in the signal.
One fact about the CA that could be exploited is that using a rectangular pulse shape, a linearly modulated signal's CA exhibits spikes not only at the signal's fundamental cycle frequency but also at the harmonics thereof.Another one is the symmetry of the CA around the DC component.First steps in this direction showing promising results have been taken in [35].The drawback of the solution proposed in [35] is that the convex optimization problem used to recover the CA becomes huge for practical parameter choices, which results in a prohibitively large computational complexity.To circumvent this, we propose an OMP-based greedy algorithm that takes advantage of the additional prior knowledge while featuring a much smaller complexity than the optimization problem.
The proposed Dice algorithm (Algorithm 2) follows the same idea as the SOber algorithm (Algorithm 1) in that, it iteratively determines the support of the sparse CA and subsequently recovers it by solving an overdetermined least-square problem.However, in contrast to SOber, Dice facilitates the use of further prior knowledge in addition to the CA's sparsity in the recovery process.Thus, in addition to the inputs received by SOber, Dice needs a set of structure dictionaries | n ν l=1 , one dictionary for each delay value ν l that is to be considered in the recovery process.Since the structure dictionaries do not necessarily model the DC component of the CA, it is added to the support set in the initialization phase in Dice (line 1).Instead of working with the amount of correlation between the residuum and the atoms directly as in SOber, the Dice algorithm computes combinations of these as dictated by the structure dictionaries in use (lines 3, 4).This way, the decision about the non-zero cycle frequencies (line 5) takes into account the structure of the CA.Additionally, instead of adding a single element to the support set per iteration, Algorithm 2 adds all indices to the support set that have a non-zero value in the selected dictionary word.The recovery step (cf.lines 6, 7) remains unchanged.Note that in Algorithm 2, the abs(•) operator stands for the element-wise absolute value of a matrix, while the selection operator [•] l: denotes the lth row of a matrix.
In the following, we introduce two particular structure dictionaries that can be used with the proposed algorithm: (i) the dictionary that accounts for the symmetry of the CA and (ii) the dictionary that describes the harmonic structure of the CA as well as its shape.

Symmetry dictionary
Let D 2 denote the symmetry dictionary.Its columns represent possible cycle frequencies contained in the set a ∈ {1, . . ., N 2 }.For simplicity, this set is chosen such that the frequencies contained in it lie at the center frequencies of the CA's DFT bins.An entry of the dictionary covers elements 1 to N 2 of rν x which is indexed from 0 to N − 1.The symmetry dictionary is simply given by the identity matrix, i. e., D To model the whole vector rν x , the dictionary is extended to include the DC component, which is set to zero, as well as the negative cycle frequencies.Note that the DC component is set to zero because its value is independent of the presence of cyclostationarity.The resulting full dictionary is exemplarily given by The circle above the symbol indicates that it is the full version of the dictionary, i. e., the one spanning the whole Fourier range.The ones in the matrix specify the locations of the non-zero coefficients in the CA fitting the format of (9).Note that in the case of the symmetry dictionary, all of sym .

Asymptotic CA and asymptotic dictionary
The symmetry structure dictionary exploits one of the facts we know about the CA.In order to explore an extreme in terms of prior knowledge, we create a dictionary that contains the maximum possible amount of prior information about the CA, i. e., the one containing the asymptotic CA itself.This requires knowledge of the analytic expression for the discrete asymptotic CA vector, which we derive in the following.
To assess the performance of different CA estimation algorithms, we employ common linearly modulated signals with symbol length T s as described by the following equation ( [7], Eq. 73) Here, p(t) is a deterministic finite-energy pulse, φ represents a fixed pulse-timing phase parameter and c n stands for the nth symbol to be transmitted.We are now interested in an expression for the discrete asymptotic CA vector of the above signal type.
We consider the case where c n is a purely stationary random sequence.Thus, its autocorrelation where σ 2 c is the average power of c n .In the following, we assume a rectangular pulse shape of length T s , i. e., p(t) = rect( t T s ), which leads to p(t . Thus, applying the Fourier transform to (26) yields for |τ | ≤ T s where sinc(x) = sin(πx) πx .Note that the use of the absolute value of the delay stems from the fact that for a real symmetric pulse shape p(t), the expression p(t + τ 2 )p * (t − τ 2 ) is symmetric with respect to τ .Equation 28 represents the CA of the continuous-time signal described by (24).The CA of the sampled version of ( 24) at its fundamental cycle frequency and the harmonics thereof is given by The derivation of this expression can be found in the appendix.
The coefficients of the closed-form expression (29) together with the alternative case R a s,n s (ν) a =k N ns = 0 at different discrete cycle frequencies a are arranged in a vector r ν s,n s [a] matching the format of the DFT matrix, such that Note that adding purely stationary noise to the signal s(t) does not change its asymptotic CA (with the exception of (a, ν) = (0, 0), at which point the CA's value is the average power of signal and noise, cf. ( 2)) since the noise exhibits no inherent periodic behavior.Due to this fact, (30) can also be used as a reference for the CA of signals contaminated with additive white Gaussian noise (AWGN) with the exception mentioned.Given (30), we can now construct the asymptotic dictionary: Note that in contrast to the single symmetry dictionary, there is a whole set of asymptotic dictionaries, one for each delay value of interest.The columns of the dictionaries correspond to actual symbol lengths i. e., actual cycle frequencies.Thus, each column contains the absolute value of the normalized asymptotic CA of a cycle frequency candidate where the discrete symbol lengths n s ∈ N 1 , . . ., N N/2 correspond to the discrete cycle frequencies a ∈ {1, . . ., N/2}.It is worth noting that in addition to its role as the basis of the second structure dictionary for Algorithm 2, the expression (30) serves as a reference for the direct comparison of different CA estimation methods in Section 6.

Cyclostationarity detection from sparse cyclic spectra
Both the SOMP-based (Algorithm 1) and the dictionaryassisted CA estimation (Algorithm 2) are able to recover the CA without knowing which cycle frequencies are contained in the signal beforehand.Furthermore, since the the (row) support of Rx corresponds to the candidate cycle frequencies, its identification in Algorithm 1 (line 4) or in Algorithm 2 (line 5) can be interpreted as blind cyclostationary spectrum sensing by itself.However, under practical limitations on the number of samples available for CA estimation and attainable SNR levels, the support estimate is likely to contain errors, e.g., missed and/or falsely identified support entries.This calls for a further testing of the candidate cyclofrequencies that can be performed by applying the TDT method described in Section 2.2.However, the obtained sparse CA is not directly compatible with the traditional TDT because, since only few of the coefficients of Rx are recovered and all other coefficients are set to zero, it is not possible to reliably estimate the covariance matrix ˆ xx * under H 0 .To tackle this problem, we present a modification of the traditional TDT.The traditional TDT is a CFAR detector, i. e., the probability density function (PDF) of its test statistic under H 0 is asymptotically independent of any signal parameters like, e. g., the noise power.To achieve this, the TDT first estimates the CA noise covariance and then rescales the original CA by this estimate so that the scaled CA follows a standard Gaussian distribution.This is where the problem occurs.Although, we are ultimately only interested in the CA coefficients that are located at the signal's cycle frequencies, for the estimation of the noise covariance, we need the coefficients lying between the cycle frequencies, which only carry estimation noise.SOber and Dice do not recover these.Thus, we propose an extension to the TDT, the sparse TDT, to bridge this gap in the following.
To obtain optimal CA recovery performance, one would choose the sensing matrix A with minimum structure, i. e., the selection of the m entries of the delay product would be completely random.However, to tackle the aforementioned problem, we choose a combination of consecutive and random delay product elements.The consecutive part comprises the first βm rows of Y m , where β ∈ [0.01, 0.5] and • denotes the ceiling operation.The remainder of the rows of Y m is a random selection of the remaining rows of Y N .The first step of the sparse TDT is to determine the classical CA estimation of the consecutive block of delay product elements.In the next step, the cycle frequency of interest a 0 is determined using either Algorithm 1 or Algorithm 2. Next, the covariance matrix for the cycle frequency a 0 corresponding to the N size CA ˆ (N) xx * (a 0 ) needs to be determined, where the superscript (N) indicates the corresponding CA size.It is obtained as where is the covariance matrix corresponding to the βm size CA estimated from the consecutive samples in the first step.The test statistic is subsequently evaluated as (cf.( 18)) The consecutive sample ratio β is a trade-off parameter.The optimal sparse recovery performance is to be expected for the case that A = NMF −1 has the smallest possible amount of structure, which here corresponds to the case where the set of known delay product elements is chosen completely at random i. e., for β = 0. Contrarily, the best estimation quality for the CA covariance matrix ˆ xx * is achieved when all known delay product elements are consecutive i. e., for β = 1.

Numerical evaluation
In this section, we compare the performance of the methods presented in the preceding sections.The parameters used throughout this section are given in Table 1.
We begin by investigating the influence of the consecutive sample ratio β on the spectrum sensing performance.Figure 1 shows how the detection rate changes with β for an SNR of 0 dB and different false alarm rates.For all methods but the OMP, β = 0.15 seems to be a good choice.For the OMP, the detection rate increases monotonically with β.However, as can be seen below, even for the OMP, a high β is no good choice regarding other performance categories.In Fig. 2, the best achievable detection rate, i. e., the detection rate for the individual best choice of β, of the different detectors is plotted over the receiver SNR for different false alarm rates.The term oracle expresses that a method has prior knowledge about the exact cycle frequencies contained in the signal.The classic method depends on this knowledge while for the sparse recovery, it reduces the CA recovery to solving the overdetermined Fig. 2 Maximum detection rate (optimal individual consecutive sample ratio selection) over SNR for different false alarm rates Fig. 3 False alarm rate that has to be selected according to the chi-squared distribution to obtain different actual false alarm rates over the consecutive sample ratio least-square problem for the given support (cf. lines 5 and 6 in Algorithm 1 or lines 6 and 7 in Algorithm 2).As expected, the oracle methods outperform the methods which have to determine the CA support themselves by a large margin.Regarding the case of missing support knowledge, the Dice algorithm clearly outperforms the SOber algorithm as well as OMP.It is to be noted that both, Fig. 1 as well as Fig. 2 do not show a significant performance advantage of exploiting the full knowledge of the asymptotic CA (Dice (asy)) over just exploiting its symmetry property (Dice (sym)) for a sensible choice of β.
The lines in Fig. 3 show which false alarm rate according to the ideal chi-squared distribution has to be set in order to achieve 1, 3, 5, and 10% false alarm rate in the actual system.The dashed lines cross at the desired false alarm rate with β = 0.15.While the two Dice methods roughly keep within a 1% offset, OMP, and SOber show a decreasing degree of equivalence for an increasing false alarm rate.This indicates that using the chi-squared distribution for setting the decision threshold of the Dice algorithm is viable, which is an important observation.It means that in contrast to many other spectrum sensing  algorithms, Dice approximately possesses a desirable feature called constant false alarm rate i. e., its test statistic is independent of system parameters like the receiver noise power.
Figure 4 shows how well the support of the CA is recovered by the different methods.Since different types of communication signals feature different cycle frequencies, this information can be used for system identification.The hit rate is the chance of exactly recovering the correct support while the absolute index error is the mean recovery error in terms of CA bins.Obviously, the Dice methods have superior support recovery capabilities.
The final performance category we evaluate is the CA estimation quality achievable by sparse recovery methods measured by the mean squared error (MSE).In the left graph of Fig. 5, the MSE over the whole CA is plotted while the right graph shows the MSE at the spikes of the CA, i. e., the MSE at the actual cycle frequencies.To determine the error, we use the analytic expression for the asymptotic CA vector as derived in Section 4.2, i. e., the MSE is defined as where • F is the Frobenius norm.The sparse recovery method has a much lower overall MSE.This is caused by the fact that it sets all CA coefficients but the detected support to zero while the classical method results in a CA that features estimation noise between the spikes.Regarding the spike MSE, both methods seem to perform roughly equivalently.

Conclusions
Blind operation and constant false alarm rate (CFAR) are desirable characteristics of spectrum sensing algorithms.Unfortunately, cyclostationarity-based approaches typically only feature either one or the other.We showed that this can be changed by using sparse recovery methods in the CA estimation.Subsequently, we developed a way to use further prior knowledge in addition to sparsity for superior CA estimation.We derived a closed-form expression of the CA of sampled linearly modulated signals with rectangular pulse shape to be used both as prior information for the CA estimation and as a reference for comparison.Finally, we extended a well-known statistical test for cyclostationarity to accommodate sparse input.
The results allow us to conclude that the proposed Dice algorithm in combination with the symmetric structure dictionary constitutes a viable alternative to the classical TDT for the case of missing prior information about the cycle frequencies contained in the signal.

Endnote
1 Note that although we use the vector of Nyquist rate samples x t 0 to calculate xt 0 , we do so for notational convenience only.In practice, one can directly obtain the sub-Nyquist samples xt 0 by means of non-uniform sampling for instance.

Discrete asymptotic CA
The relation between the CA of the continuous timedomain signal s(t), and its sampled counterpart {s(nT e )} is given by ( [36], Ch. 11, Sec.C, Eq. ( 111)) The sum over l reflects the infinite aliasing caused by the sampling.In the next step, we insert (28) into (35).Also, we express quantities in terms of the sampling period T e , i. e., T s → n s T e , α → a NT e , φ → d φ T e , with n s , a, N ∈ Z.This leads to •e j2πld φ sinc(( a N + l)(n s − |ν|)) otherwise, (36) for |ν| ≤ n s , where n s is the oversampling factor.In this step, we used the fact that for our assumptions all aliases of the fundamental cycle frequency and its harmonics lie on top of the actual fundamental cycle frequency and its harmonics, i. e., (α + l T e )T s ∈ Z iff αT s ∈ Z. Inserting the discrete quantities given above, we get ( a N + l)n s ∈ Z iff a N n s ∈ Z. Since n s ∈ Z and l ∈ Z, this always holds.To rule out any spectral leakage, we choose N as an integer multiple of n s , since then, a = k N n s is also an integer and thus the fundamental discrete cycle frequency and its harmonics hit center frequencies of frequency bins.
For a = k N n s expression (36) can be shown to be To obtain (37) we used the definition of the sinc and exploited the facts that e jπk = (−1) k for k ∈ Z and that sin(x + kπ) = (−1) k sin(x) for k ∈ Z.The pulse timing phase parameter d φ was set to n s +1 2 .This has the following reason.In order to simplify the numerical evaluation, we want to choose φ such that the beginning of the observed receiver signal is aligned with the rectangular pulse shapes, i. e., we would set φ = T s 2 .However, doing so would lead to the need to sample at the discontinuities caused by the instant change in amplitudes at the transition between symbols.To avoid this, we choose φ = T s 2 + , where ∈ (0, T e ).Note that (37) is the same for any ∈ (0, T e ).In order to ease the derivation, we can thus choose = T e 2 , i. e., d φ = n s +1 2 .The infinite series in (37) can be expressed as and ( [37], Eq. (6.3.5)) respectively, we obtain and Inserting ( 43) and ( 44) into (40) results in Finally, substituting (45) into (37) gives us the expression (29).

Fig. 4
Fig.4 Left hit rate over SNR, right absolute index error over SNR.Both at consecutive sample ratio 0.15

Table 1
System parameters Detection rate over consecutive sample ratio for different false alarm rates at 0 dB SNR