{'fraction_non_alphanumeric': 0.04153698826909265, 'fraction_numerical': 0.008331338281062963, 'mean_word_length': 4.183149274103487, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonlinear Distortion Reduction in OFDM from Reliable Perturbations in Data Carriers

A novel method for correcting the effect of nonlinear distortion in orthogonal frequency division multiplexing signals is proposed. The method depends on adaptively selecting the distortion over a subset of the data carriers, and then using tools from compressed sensing and sparse Bayesian recovery to estimate the distortion over the other carriers. Central to this method is the fact that carriers (or tones) are decoded with different levels of confidence, depending on a coupled function of the magnitude and phase of the distortion over each carrier, in addition to the respective channel strength. Moreover, as no pilots are required by this method, a significant improvement in terms of achievable rate can be achieved relative to previous work.


Introduction
Multicarrier signaling schemes such as orthogonal frequency division multiplexing (OFDM) are highly susceptible to nonlinear distortion at all stages of the transmission process. This is partly due to the impulsive nature of these signals in the time domain, where the superposition of modulated waveforms takes place. When nonlinear distortion is confined to the transmitter, many proposed methods in the literature use reserved carriers (or tones) to carry information about this distortion to the receiver, at the obvious cost of reducing data-rate [2][3][4][5]. The central idea in these techniques is to construct clipping (or peak-reducing) signals by performing a constrained search at the transmitter, one which confines the frequency support of the clipping signal to the reserved carriers, while reducing the peaks of the data signal in the time domain. These approaches are not generally robust, as they demand that the frequency support of the data and clipping signals remain strictly disjoint throughout the transmission process, and add significant complexity at the transmitter.
To combat this, techniques based on compressive sensing (CS) that were tuned to clipped OFDM models were proposed in [6]. These techniques removed the need for any constrained search at the transmitter, since the receiver could detect the entire clipping signal by observing a subset of its frequency components available on the reserved tones, provided that the signal is sparse in time. Consequently, the need to maintain orthogonality in frequency was completely relaxed, but the need for a significant amount of reserved carriers persisted. To avoid this loss in data-rate, the authors in [7] proposed using the channel estimation pilots for this purpose. Nonetheless, this approach severely limits the number of measurements available to the CS algorithm and hence its ability to deal with severe clipping scenarios. It also does not make use of available information such as clipping likelihood and phase resemblance in the time domain (i.e., phase resemblance between the clipped and clipping signal) as done in [6].
In this paper, a fundamentally different approach to these methods is pursued. Specifically, in contrast to the authors' previous work on the topic of using CS concepts in OFDM [6], the technique presented in this paper does not require any orthogonality between the frequency support of the data and the distortion, and no tones (null, edge, channel pilots, or otherwise) are needed either. The receiver is free to select which and how many data tones it will use to read off differential observations, and will use them to estimate and cancel the entire distortion over the tones. In addition, no data-rate is lost by employing the proposed strategy. Furthermore, a significant tradeoff also exists in regard to complexity, distortion tolerance, and robustness to channel estimation errors, so that the user has many algorithms to use within the proposed framework. Similarly, in contrast to another recent work [1], the current paper introduces an entirely new and rigorous way of analyzing the reliability of tones. In addition, it also introduces a method to finetune the performance of CS by minimizing the probability of incorrect measurements and maximizing the clipping-to-noise ratio (CNR). This framework is made possible by jointly taking three phenomena into account. The first is that not all tones carry correct decoding information to the receiver. The second is that the receiver can probabilistically assign levels of confidence to each tone, and the third is that the distortion is sparse in the time domain. These phenomena motivate us to employ CS and sparse recovery techniques much more effectively compared to previous techniques, as we can quantify the perturbations on data tones, select the most reliable ones, and then most-importantly, use the power of CS techniques to recover the significant time domain distortions. Our major contributions include: 1. Formulating CS models within a pilotless transmission framework. 2. Proposing systematic methods to adaptively select the tone subset to sense over, among the combinatorially large possibilities of pilotless CS models. 3. Developing a novel method for accurately assessing the reliability of estimated coefficients based on their symbol-wise magnitudes, phases, and relative locations to other constellation points, as well as channel strengths. 4. Deriving a closed-form expression that characterizes the modes of behavior of the reliability function, and devising geometrically-inspired approximations based on this expression for quick and efficient selection of the tone subset over each OFDM block. 5. Proposing dual-stage construction of the tone subset, where the first stage minimizes the probability of incorrect measurements, while the second maximizes an CNR metric to optimize CS performance. 6. Providing probabilistic upper bounds for choosing the number of tones in the CS model without risking incorrect measurements.
The remainder of the paper is organized as follows. Section 2 briefly describes the transmission and distortion models. Section 3 demonstrates how a pilotless CS model can be derived within the previous transmission model. Subsequently, Section 4, the heart of the paper, focuses on selecting the subset of tones used for CS. This includes developing reliability assessment criteria, deriving analytical approximations for quick and efficient assessment, selecting the number of tones, establishing dual-stage subset selection to maximize CS performance, and, finally, condensing the major results into an algorithm. Section 5 presents our simulations and Section 6 concludes the paper.

Notation
We use regular font for scalars and boldface letters for matrices and vectors. To distinguish between vectors in the time and frequency domains, we use boldface calligraphic notation for vectors in the frequency domain (e.g. X , X , Y) and boldface lowercase letters for their corresponding time domain representations (e.g. x, x, y). We use X (k) to denote the kth coefficient of X , or more simply X , when it is clear from the context. Moreover, we use X Ω to represent a vector formed by selecting the coefficients of X indexed by set Ω. Similarly, Y Ω is the vector formed by indexing the corresponding elements of vector Y Ω according to the index set Ω. We further define S Ω to be a diagonal binary selection matrix, with |Ω| number of 1's at locations along its diagonal specified by the tone set Ω.

Transmission and Clipping Model
In an OFDM system, serially incoming bits are mapped into an M -ary QAM alphabet A = {A 0 , A 1 , . . . , A M −1 } and concatenated to form an N -dimensional data vector, Since x has a high PAPR, the digital samples are subject to a magnitude limiter that saturates its operands to a value of γ. Hence, instead of x, we feed x to the power amplifier, where and where θ x(i) is the phase of x(i). This hard-limiting operation can be conveniently thought of as adding a peak-reducing signal c to x so that its low-PAPR counterpart x = x + c is transmitted instead. Furthermore, by setting a typical clipping threshold, γ, on x, c is controllably sparse in time by the impulsive nature of x, and dense in frequency by the uncertainty principle. We denote the temporal support of c by I c = {i : c(i) = 0} and always maintain the practical assumption that |I c | N . Subsequently, x is convolved with a channel of impulse response h ∼ CN (0, σ 2 h I L h ), and subjected to additive white Gaussian noise (AWGN) z ∼ CN (0, σ 2 z I N ), where L h is the length of channel impulse response. Equivalently, in the frequency domain, this translates to transmitting with complex coefficients that are now randomly pre-perturbed from the lattice A N , followed by additional multiplicative perturbations by the channel H and additive perturbations by the noise Z ∼ CN (0, σ 2 z I N ) at the receiver. By virtue of the added cyclic prefix (of length > L h ), the circulant channel matrix H can be decomposed and expressed as H = F H ΛF where Λ is an N × N diagonal matrix composed of the frequency-domain channel gains, {λ(k)} N k=1 . As a result the frequency domain received signal reads Y = ΛX + Z, where, for the moment, we make the practical assumption that the channel coefficients are known at the receiver. Consequently, X can be directly recovered scalar-wise from Y, i.e.,X where we use the notationX to represent the equalized estimate of X at the receiver. Writing (3) in vector notation yieldŝ Let D(k) C(k) + λ −1 (k)Z(k) denote the general distortion on X (k), We let f D to be the pdf of the general distortion D which we assume to be zero mean circularly symmetric Gaussian with variance σ 2 D . (Please refer to Appendix A for details regarding the derivation of σ 2 D .) Equation (4) could now be written asX Treating the clipping distortion as additive noise, an maximum likelihood (ML) decoder will recover X (k) by simply mappingX (k) to the nearest constellation In other words, the operation X (k) corresponds to roundingX (k) to the nearest constellation point. Note that X (k) does not need to be the true constellation point, i.e., X (k) might be different from X (k). Such a hard-decoding scheme is very efficient in the classical AWGN scenario for high signal-to-noise ratio (SNR). However, in our case, in addition to the additive noise λ −1 (k)Z(k), we have a γdependent source of perturbation C(k) which is independent of the SNR. CS and similar sparse recovery algorithms seem to be a very sensible solution towards recovery of C(k). Since c is sparse in the time domain, a partial observation of c in frequency domain is sufficient to estimate c and hence C in one shot. This would certainly get around the problem of unreliable perturbations as CS algorithms, for instance, can be totally blind to them and still offer near optimal signal reconstruction under mild conditions [8]. The main issue is to decide which partial observation to use. This will be the topic of the following section.

Development of Compressive Sensing Models with No Tone Reservation
With the addition of the general distortion vector D to the data vector X , we expect that part of the data samples will be severely perturbed such that they fall out of their true decision regions. Let X (k) denote the decoded data sample corresponding to X (k), then the true decision region for X (k) is defined as Q(k) { X (k) + U ∈ C : X (k) + U = X (k)} where U is a factor which when added to X (k) keeps it in its true decision region. Moreover, denote by Ω T = {k ∈ Ω : X (k) + D(k) = X (k)} the subset of data tones in Ω = 1 While Am refers to the mth constellation point (1 ≤ m ≤ M ), we reserve X (k) to denote the nearest constellation point corresponding to the kth received data sampleX (k). Furthermore, note that the true constellation point corresponding toX (k) is X (k).
{1, 2, . . . , N } in which the perturbations do not cause data samples to cross their true decision regions. LetΩ T = Ω\Ω T be its complement. Over the data tones of Ω T , the equality in X (k) = X (k) is true and hence from (2) and (5), More generally, we can write where S Ω T is an N ×N diagonal binary selection matrix, with |Ω T | number of 1's at locations along its diagonal specified by the tone set Ω T . It extracts the elements of the vectorX − X according to the tone set Ω T while nulling the others. The matrix SΩ T is similarly defined with 1's along the diagonal specified by the setΩ T . It is easy to see that S Ω T SΩ T = 0. Practically speaking, Ω T constitutes the larger part of the tone set Ω. An essential part of OFDM signal recovery obviously constitutes finding Ω T and correcting the distortions over Ω T to finally reach the state Ω T = Ω. From (5) we have at the receiver, which is the analog estimateX of the data vector X affected by the distortion D. We define E X − X to be a vector that is nonzero at locations where the decoded estimate at the receiver X differs from the data vector X . From the discussion above, we see that which allows us to write Note that we do not require all of Ω T to recover c. Rather, we only require an arbitrary subset Ω m ⊆ Ω T ⊆ Ω of cardinality m |Ω m | ≤ |Ω T | to correctly recover c by CS (m and |Ω m | will be used interchangeably as appropriate to denote the number of measurements). As a result, we can replace the equation above with S Ωm (X − X ) = S Ωm D = S Ωm Fc + S Ωm Λ −1 Z, where S Ωm is also a N ×N diagonal binary selection matrix defined in a similar way as S Ω T in the above. We write the above equation simply as Y Ψc+Z , where Ψ S Ωm F, Z S Ωm Λ −1 Z, and Y S Ωm (X − X ) which denotes the observation vector of the differences over the tones in Ω m , nulled at the discarded measurements.
This leads us to a pilotless CS model 2 where Y Ωm is the |Ω m |-dimensional vector composed of the nonzero coefficients in Y . By inspecting (10), we notice that c is an N -dimensional sparse vector in the time domain, corresponding to the difference between the time representations of the OFDM signal x and its clipped counterpart x. The matrix Ψ Ωm ∈ C m×N is obtained by m random row extractions from the N×N Fourier matrix according to Ω m (the cause of randomness is discussed later). The m-dimensional vector Y Ωm is the corresponding partial frequency-domain observation that we use to estimate c, contaminated by the Gaussian noise vector, Z Ωm .
This is a standard model in CS [8,9]. Note however that the parametrization by Ω m actually defines a huge set of 2 N possible models 3 . In the forthcoming sections we will discuss in detail how to determine a proper model from all these possibilities. For the time being, we assume that an appropriate Ω m is chosen, and c could therefore be recovered using any CS technique, be it convex programming, greedy pursuit, or iterative thresholding, and a very flexible region for tradeoff exists in regard to the performance and complexity of these techniques.
In this paper, we use two different schemes of CS to recover c from the developed CS model in (10), one from the convex relaxation group and the other from greedy pursuit methods. More specifically, the first is an adaptation of the least absolute shrinkage and selection operator (LASSO) [10] to this problem called the weighted and phase-augmented LASSO (WPAL) [6]. It incorporates phase information and clipping likelihood available from the data in the time-domain to improve distortion recovery performed in the frequencydomain. Specifically, we know that c is composed of just the clipped portions of the transmitted signal x, so at clipping locations ∠c(k) = −∠x(k). We also know that the closer |x(k)| to the value of the clipping threshold γ, the higher the likelihood that c had an active coefficient at k. This additional information is incorporated in the CS algorithm in the form of weighting to improve its performance. Therefore, we define w ||x|−γ| T to be such a weighting vector to the 1 -norm of c in the LASSO wherex refers to the estimated received clipped signal. We further define the diagonal phase matrixΘ c = − exp diag(θx) such that c WPAL =Θ c |c WPAL |. With these two variables defined, the optimization problem we solve becomes for some noise-dependent parameter . The other technique is the fast Bayesian matching pursuit (FBMP) by Schniter et al. [11] chosen for its superior performance and efficiency when a relatively large number of measurements is available, which is indeed the case compared to tone reservation cases proposed in [6] and [7]. Finally, once c CS -the CS estimate of c -has been obtained through any of the abovementioned schemes, we use (2) and (4) to obtain where E CS C − C CS and ∆ E CS + Λ −1 Z.X could be used to determine the data vector X exactly, provided that no ∆(k) causes crossing of X (k) out of its ML decoding region (this issue will come up in Section 4.3). Our subsequent objective is to scrutinize the general conditioning of the model itself by supplying our most reliable observations to the generic CS algorithm.

Cherry Picking Ω m
An essential question now is how to select among the 2 N possible Ω m (or N m if m is fixed) in order to compute c CS . In this connection we devise a reliability function which associates a reliability estimate with each tone and thus lets us determine the m most reliable tones to construct Ω m . A general strategy of CS techniques is to select these m tones randomly for near-optimal performance [8]. Although possible in our scenario, such a strategy neglects the fact that our observations vary in their credibility and attest to whether or not they represent true frequency-domain measurements of C. 4 Since we deal with each tone separately in what follows, we henceforth drop the k index while preserving the italic notation to emphasize the scalar-wise operations in this section. With the receiver risking faulty decisions, it must devise a procedure to select the most reliable set of observations over which to sense. To this end consider the estimateX and the nearest constellation point X . The latter is in general surrounded by eight points which either belong to the set of nearest neighbors (NN) or the set of next nearest neighbors (NNN) as illustrated in observations could be done based on the relative posterior probability that D equalsX − X to the probability that it equals some other difference vector Pr( X =ANN) . The higher the value of R, the higher the reliability that X = X as relative to the fact thatX should decode to A N N . From (7), we can see that R = Pr(D=X − X ) Pr(D=X −ANN) Now as mentioned earlier, we model D to be Gaussian circularly symmetric with variance σ 2 D , then f D (X ) = 1 and we can write, Intuitively, the minimum certainty occurs at the boundary of the decision region and attains R min = 1. At such tones, we would be highly skeptical of whether D =X − X or D =X − A NN , and we would hence supply a plausibly false measurement to the CS algorithm. To avoid such unreliable measurements, assume we only choose the tones with respective perturbationsX − X that are confined in the complex plane to a disk of radius r o (i.e., |X − X | ≤ r o ). In such a case, given the minimum distance between any two constellation points (d min ), the minimum reliability would increase to R min = for the next nearest neighbor A NNN when it points in the direction of a decision region's corner. So while bothX 1 andX 2 have the same distance r o from X 1 = X 2 ,X 1 has a higher reliability thanX 2 as it is farther from the nearest neighbor. Fig. 2 shows a part of constellation to illustrate the idea. This suggests a need to factor in the direction or phase of the perturbation, θX − X , in assessing its reliability in addition to the magnitude-dependent pdf f D . Defined axiomatically, the reliability of a measurement at each tone is then a function R that maps a 3-tuple (|X − X |, θX − X , λ −1 ) into R [1,∞] (i.e., a function of the magnitude of the observation, its phase and the channel gain at that tone).  (13) to sense over, where R i:N denotes the ith-order statistic in the vector R [12]. With this selection of Ω m , the locations of the measurement tones correspond to the indices of the highest m-order statistics of N random variables in R. As mentioned previously, each of these variables is a function of the 3-tuple above, and whereas the first two are uncorrelated across the tones, this does not generally hold for the third, i.e., λ(k).
In fact, assuming L h channel taps with a uniform power-delay profile, then the absolute autocorrelation of the channel across tones k and l can be expressed as |E[λ(k)λ H (l)]| = sinc (πL h (k − l)/N ) [13]. Hence, only for sufficiently large L h can we assume that the channel gains are uncorrelated. Otherwise, the set of reliable tones Ω m deviates from a uniformly random tone selection model typically assumed in the literature [8,9], and reliable tones would instead come in clusters corresponding to strong channel gains. The efficiency of CS in this case might be reduced 5 .

Criteria for Evaluating R
Using the reasoning based on the scalar-wise likelihood ratio defined in (13), an exact expression for the reliability could be a direct generalization of (13), namely, Unfortunately, this pursuit for exact reliability computation is inefficient since it requires NM evaluations of f D (·), which grows with the constellation size M , even though the probability of a perturbation exceeding the first tier of eight surrounding neighbors (i.e., the nearest and next nearest neighbors illustrated in Fig. 1) is insignificant. As such, we can truncate the computation above to the first tier, denoting its result by R trunc , with a minor effect on the performance. Two simpler reliability functions are also worth mentioning. The first is solely based on the probability of the perturbation (i.e., f D (|X − X |)) and is hence completely blind to the direction of the perturbation in the constellation's plane, while the other one intuitively takes this extra phase information into account by defining a square centered at X as the reliable region, hence having the ability to favor perturbations with larger magnitudes if θX − X were close to π/4, i.e., if they pointed to the next nearest neighbor. We denote these reliability functions by R • and R 2 , respectively, motivated by the geometric shape they define. In the next section, a more rigorous approach is taken to justify when such simpler methods can be used.

Analysis of Truncated and Approximate Reliability Criteria
Dropping the tone index, assume thatX is an observation that falls among four points in an M -QAM constellation such thatX ∈ Q. LetX − X re θ be the polar representation of this point with the origin at X , such that r |X − X | ∈ [0, 1 √ 2 d min ] and θ θX − X ∈ [0, π/2]. We are interested in a more abstract expression of the truncated Bayesian reliability function, R trunc , one that defines its output by only acting onX − X while taking the relative position in the constellation implicitly into account as well 6 .
By the definition of r, and by referring to Fig. 3 which again shows a part of the constellation diagram (similar to Fig. 2), the distances betweenX and the other three competing constellation points are and We neglect detailing their phases, θ 1 , θ 2 , and θ 3 , since they have no effect on  the results, although r 1 , r 2 , and r 3 are clearly functions of r and θ as portrayed in Fig. 3. In effect, R trunc reduces to Canceling out the common function f D (r) appearing in all terms above, and collecting common terms yields for the first quadrant, θ ∈ [0, π/2], in the complex plane, where α exp(2d min r/σ 2 D ) and β exp(−d 2 min /σ 2 D ). Although clearly a function of r, we will treat α as a constant (i.e., evaluated at a fixed magnitude, r) when we wish to focus on R trunc as an explicit function of θ, say g(θ) 7 . This function is symmetric about θ = π/4 and exhibits quite complex behavior with r and θ as indicated in Fig.  4.
Most importantly, our concern will be whether R trunc (r, θ) is convex or concave with respect to θ at different regions of r. This is because R trunc (θ) undertakes a fundamental shift in behavior as r varies from 0 to 1 √ 2 d min , and its approximation by basic trigonometric functions and geometric objects such as squares and circles depends on whether it is convex or concave with respect to θ. Notice first that when r σ 2 D /2d min , α ≈ 1 and hence R trunc ≈ 1 β[2+β] becomes relatively isotropic (i.e., independent of θ) and therefore akin to R • . Referring to Fig. 4 for example, the polar curve of the normalized reliability function R trunc (θ)| r=0.1dmin , evaluated at the smallest magnitude r = 0.1d min , confirms this observation. (Refer to the curve with blue circular markers in Fig. 4).
In fact, as will be shown shortly, R trunc (r, θ) will tend to even disfavor perturbations along π/4 (or π/4 + nπ/2, n = 1, . . . , 3 in general) until it shifts gears and takes on its concave behavior with respect to θ. For instance, the curve of the normalized function R trunc (θ)| r=0.2dmin appearing in Fig. 4 is actually convex (therefore assigning slightly higher reliability to perturbations in the direction of the next nearest neighbors having the same magnitude 0.2d min ), whereas the curves of the normalized function R trunc (θ)| r=0.3dmin and beyond are concave. To pinpoint the location of this behavioral shift, we need to find where ∂ 2 R trunc (θ)/∂θ 2 is expressed as ∂ 2 R trunc (θ) ∂θ 2 = 1 β 2 β ln α(cos θ − sin θ)α sin θ+cos θ − ln α sin θα cos θ + ln α cos θα sin θ 2 (βα sin θ+cos θ + α sin θ + α cos θ ) 3 − β ln 2 α(cos θ − sin θ) 2 α sin θ+cos θ + β ln α(− sin θ − cos θ)α sin θ+cos θ − ln α sin θα sin θ − ln α cos θα cos θ + ln 2 α cos 2 θα sin θ + ln 2 α sin 2 θα cos θ / βα sin θ+cos θ + α sin θ + α cos θ 2 (20) Clearly, this is a daunting task. However, we can reduce it to finding the root,r, which satisfies (19) when evaluated at θ = π/4, since our main concern is whether or not R trunc (r, θ) will be tapered along θ = π/4, and this will fortunately result in many cancelations in (20) due to symmetry about this point. Pursuing this reduces (19) to solving Expanding this into the original parameters implies that we need to findr such that √ 2d minr σ 2 which cannot be solved explicitly in terms ofr. Rather, by means of a proper substitution, it can be put in the implicit form g(r)e g(r) = q, where q is independent ofr and expressed using the primary branch W 0 of Lambert's W-function [14]. The explicit solution to the previous form can be expressed as g(r) = W 0 (q), and the desired explicit expressionr is obtained by back-substitution (Refer to Appendix B for details). Ultimately, we can show that Furthermore, as W 0 (0) = 0 and d min /2 > σ 2 D , it is clear that d 2 min /σ 2 D > 2d min , and that the argument of W 0 quickly approaches zero from the left as σ 2 D diminishes, resulting in the following accurate approximation of (23): for small σ 2 D relative to d min /2. Fig. 5 plots (23) and its approximation (24) as functions of σ 2 D . Using the approximation for simplicity,r then splits the behavior of R trunc (r)| θ=π/4 into two regions, supported by the intervals This result explains why R trunc (r, θ) in Fig. 4 first resembles a circular shape akin to R • and then inflates along the diagonals deforming into a square-shaped reliability region, R 2 , as can be seen when r = 0.3d min in Fig. 4 (see the black curve with square-shaped markers). Subsequently, as the magnitude, r, of the perturbation increases and approaches the decision boundaries, R trunc (r, θ) inflates outwards, resembling pointy leaves that can be modeled as µ + (1 − µ) cos(4θ + π), where µ ∈ [1/2, 1] > 1 − µ > 0. The analysis also provides restrictions for when square-like reliability regions suggested in the literature (such as [15]) can be justifiable.

Dual-stage construction of Ω m
The reader will notice that our primary objective so far in selecting Ω m was based on minimizing the probability of incorrect measurements, i.e., This is no doubt a necessary choice to preserve the success of the recovery algorithm as a whole, although we know that a more generic criterion, that is, one that is not at risk of using incorrect observations, would seek the tones with the maximum CNR 8 , i.e., (see (10)) Obviously, there is a conflicting interest between (25) and (26), as the former frequently seeks smaller perturbationsX − X , since they are the most likely to equal D, while the latter seeks the largest perturbations to maximize the CNR for enhanced estimation performance. This prompts us to consider a second recovery stage (i.e., another CS iteration) that uses (26) and produces a new subset of selected tones, denoted by Ω CS m . The second recovery stage takes in Ω m obtained from (25) and uses it with any of the mentioned CS recovery algorithms to get a CS estimate of C, denoted by C CS . This lets us achieve the corrected decoding decision X − C CS allowing us to have a higher confidence that D =X − X − C CS compared to the primary assumption that D =X − X . This is because the decoding error in X = X + D = X + C + λ −1 Z depends on the value of C, whereas the error in X − C CS = X + ∆ = X + E CS + λ −1 Z depends on the estimation error E CS = C − C CS of C which is expected to be smaller than C itself. These results follow from (3) and (12).
As illustrated in Figs. 6 and 7, it is possible now to use these carriers that have the largest values of the perturbations,X − X − C CS (or even the carriers with the largest values of C CS ), as the new CS measurements, without worrying much about how closeX is to the decision boundaries. Note, however, that we never have access to C, Z, or E CS , and therefore we can rely only on observable variables, such as D and ∆, to practically maximize the CNR. More importantly, these variables themselves are not always obtainable, since it is not necessarily the case that D =X − X − C CS or that ∆ =X − C CS − X − C CS (which is the main reason that we repeat CS over a subset of measurements). Instead, we have to rely on the observable variables,X − X − C CS and C CS , for this task. The availability of these two observable variables grants us the flexibility of computing CNR in two different ways, both of which are suitable for different scenarios. The first is suited for a high SNR and low CS estimation quality (see Fig. 6), while the second is suited for a low SNR and high CS estimation quality (see Fig. 7). Dropping the coefficient index, k, these are: 1. High SNR, large E cs : where the symbol means that the expression on the left hand side of this symbol is asymptotically equal to the expression on the right hand side. Figure 6: Computing CNR for the case of high SNR and low CS estimation quality 2. Low SNR, small E cs : Although the second ratio CNR E CS more vividly resembles the CNR defined in (26), the first ratio is more effective in this work as the inherent complexity of CS based methods justifies itself in severe clipping scenarios and hence expectedly higher CS error (i.e., large E CS ). Consequently, we select the differential measurements corresponding to the maximum |Ω CS m | ratios 9 , Figure 7: Computing CNR for the case of low SNR and high CS estimation quality and the new CS model is which produces an improved estimate, C rev . This new estimate of the clipping distortion from a different subset of reliable and stronger measurements can then be subtracted fromX , and a revised set of decoding decisions, X − C rev , can be made.

On selecting the cardinality |Ω m |
In the method we proposed, it is assumed throughout the above discussion that the observations supplied to the CS algorithm were true measurements of the actual perturbations caused by clipping and additive noise. Although this is a mild restriction in practice, a guarantee must nonetheless be established that any selected Ω m and |Ω m | according to Section 4 will not result in CS failure. This requires that the number of measurements |Ω m | be both small enough to minimize the probability of incorrect measurements, and also large enough relative to the sparsity level of c to ensure meeting recovery bounds of CS.
To this end, we derive a simple lower bound on the probability of Ω m ⊂ Ω T by deriving a lower bound on the probability that D is indeed equal toX − X . We study this for two cases: whenX − X is observed within a disk of radius r o from X , i.e., Ω m = {k ∈ Ω : |X (k)− X (k) | < r o }, and when it is observed within a square of side length 2r o centered at the QAM symbols A. We focus on the i=0S i . Our objective is to select |Ω m | such that Pr Ω m ⊂ Ω T |Ω m | is high given a minimum amount of required measurement for CS success.
Dropping the tone index, this will require finding Pr(D ∈ S), which requires evaluating an integral of f D over (non-centered) discrete disks in the complex plane (see Fig. 8). Since this is difficult, we will use an upper bound based on evaluating f D over centered disks that cover these regions and then slice out the irrelevant regions. More specifically, the integral D∈Si over the disk S i of a nearest neighbor will be bounded by the integral over the area highlighted by the shaded region in Fig. 8. This region covers the difference between the outer and inner sectors defined by radii d min + r o and d min − r o , respectively, and a common angle, θ s = 2 sin −1 ro dmin . In effect, the area ofS i is strictly less than θs for the nearest neighbor. Consequently, In the case where the average distortion is large and square regions S 2 = M −1 i=0S 2 i of side-length 2r o are used, pursuing the same logic above we just replace the ratio in (31) with and raise it to the power |Ω m | to obtain Pr(Ω m ⊂ Ω T |Ω m |), where Q is the familiar tail probability function. The user must then choose |Ω τ m | = arg max |Ωm| [Pr(Ω m ⊂ Ω T |Ω m |) > τ ] where τ is selected so as to supply as much information to the CS algorithm as possible while remaining in a safe region of correct measurements. Furthermore, given a clipping threshold, γ, we have an expected sparsity level, E[|I c |], and variance σ 2 |Ic| , which need to be taken into account when using sparse recovery techniques [8]. We will denote this minimum required number of frequency observations to recover an |I c |-sparse vector in time by |Ω γ m | to stress its strong dependence on γ, and take |Ω m | = max(|Ω τ m |, |Ω γ m |). The same can be done with the optional second stage |Ω CS m |. Suppose however, after taking all the protective measures thus far, that an incorrect measurement was nonetheless supplied to the CS algorithm. Does this result in CS failure? Luckily, in this application the answer is no. Recall first the decoding-error vector E X − X used to motivate (9). When a decoding error is made at the kth coefficient, an incorrect measurementX (k) − X (k) = X (k) + D(k) − X (k) D(k) + E(k) is supplied and it follows from (9) that the incorrect measurement has no impact on the performance of CS. Note that the nonzero entries of E are quantized and bounded since X and X ∈ A N . Furthermore, assuming most equalized coefficientsX (k) are decoded correctly, E is also sparse. The general differential model (i.e., the one that is not confined to the carriers in Ω m ) becomesX The first term added to C is the (dense) Gaussian noise vector, while the second is a structured -noise term. Equation (34) matches the model in [16], and it is demonstrated therein how it is possible to recover C from such noise via variations in the CS algorithm. The main results of the paper are summarized in the algorithm listed in Table 1.

Simulation Results
In this section, we perform several different experiments to show the effectiveness of the proposed technique. The methods proposed in this paper were tested on an 256 subcarrier OFDM signal, modulated using 64-QAM. The signal was subject to a block-fading, frequency-selective, Rayleigh channel model, subject to varying noise and clipping ratios (CR) defined as CR= γ/σ x [17]. Here, σ x is the standard deviation of the OFDM signal. Available packages for convex programming [18] and Fast Bayesian Matching Pursuit (FBMP) [11] were used to implement the sparse-recovery algorithms. The undistorted phase property (as utilized in [6]) is utilized while implementing FBMP and hence the modified version is termed phase augmented-FBMP (PAFBMP). Lastly, we refer to the second stage defined in 4.3 as a corrective action on the first estimation operation, and label its output by C-WPAL or C-PAFBMP for the cases when WPAL [6] and PAFBMP [11] are used, respectively.    Figure 10: Achievable Rate as a function of the clipping ratio CR= γ/σ x .

Comparison of Reliability Criterias
In the first experiment, the reliability criterias proposed in this work are compared, including R exact , R trunc , R • and R 2 . As a performance metric, we use normalized success rate (NSR) defined as |{k : X (k) = X (k), k ∈ Ω m }|/|Ω m |. The NSR depicts that among the |Ω m | tones favoured by a particular reliability criteria, how many were actually within their corresponding correct decision regions. Fig. 9 shows the result of this experiment. The results were plotted against a varying CR while 64 most reliable carriers were sought keeping E b /N 0 fixed at 20dB. It is expected that as CR is increased, all reliability criterias will tend to improve. The simulation results confirm this intuition and also confirm the conjecture that the truncated reliability computation comes at little cost compared to using the exact reliability function (14). Furthermore, it is shown for the parameters used for penalizing the circular reliability function, R • , by µ 2 = 0.95 (which results in a square-like function such as the curve plotted in Fig. 4 for r = 0.3d min ) was more effective than with a smaller value of µ 1 = 0.65, for milder CR levels. On the other hand, with severe CR, smaller values of µ were better. This observation is highlighted by showing enlarged version of the graph for severe and milder CRs.

Achievable Rate
In this experiment, the ultimate performance measure considered was the achievable rate with and without the proposed method. Assuming ergodicity over the subcarriers, this rate can be expressed as 1 for the un-mitigated case, and as 1 , for the case when an estimateĈ of C is obtained by an arbitrary method [2]. As a benchmark Oracle-LS is utilized, where the support of clipping signal is perfectly known and the active elements are estimated using LS estimate at the receiver. By using the WPAL (11) [6] and PAFBMP [11] sparse recovery algorithms over the exact reliability function R exact , we compared the clipping mitigation results. The results in Fig. 10 show the superior performance of C-PAFBMP [11] using the most reliable 39% of carriers and the ability of C-WPAL to beat all the techniques at CRs around 1.5 using the most reliable 39% carriers.

Complexity
A practical comparison of relative run times for the tested algorithms is reported in Fig. 11. All times are scaled and represented as percentage of the maximum time required for recovery. We observe that PAFBMP has the least complexity among the compared schemes and this complexity reduces with increased CR.

Conclusion
A novel method for nonlinear distortion mitigation using pilotless sparse recovery techniques was proposed. The method exploits the sparsity of the distortion in time domain to fully recover the signal without being influenced by incorrect ML-decoding decisions. The method adaptively senses over reliable subsets of observations of the distortion in frequency domain to perform the recovery. A new method of computing the reliability of each observation independently of the other M −1 candidates within a constellation was also proposed and tested. Through simulations, it is verified that the proposed scheme provides favourable results for clipping signal recovery and achieves a rate close to Oracle-LS based recovery.
Appendix A. Deriving f D Beginning with the fact that can be made, and we subsequently work out the energy of the sparse vector, c, which is a compound random variable. By total expectation we have where we have dropped the index in E |c | 2 to denote an arbitrary nonzero coefficient of c. Using this result, we can see that |x| . Substituting these last two expressions into (A.2) gives σ 2 C , which in turn produces (A.1), and this parameter characterizes f D .

Appendix B. Tailoring the Lambert W-Function
Given the equation y = xe x , which for reasons stated below we call the canonical form, it is desired to solve explicitly for x. Unfortunately, this could not be done using elementary operations. Instead, the solution can be expressed in terms of the Lambert W-Function, where x = W(y) is the solution to the canonical form, and where we can thus equivalently write y = W(y)e W(y) .
The function is generally multivalued. If we restrict its argument, y, to be real, then it produces two outputs for each point on the supporting interval y ∈ [−e −1 , 0], which is our interval of interest. However, one of the two outputs of W(y) is ≥ −1, and is referred to as the primary branch, W 0 (y), while the second is < −1, and is referred to as the secondary branch, W −1 (y). It will soon be clear that only the primary branch is needed (hence an injective mapping between y and x is retained). In any case, the function can be found iteratively by Newton's method; for instance, x j+1 = x j − xj e x j −y e x j +xj e x j . The problem at hand, as expressed in (22), is more complex than the canonical form. Nonetheless, it could be reduced to this form by a clever substitution [14].