Eigenvector-based initial ranging process for OFDMA uplink systems

At the contention-based synchronization, e.g., initial ranging process, it is crucial to identify multiple users through ranging codes and estimate the corresponding parameters such as timing offset and frequency offset. This paper presents an improved parameter estimation and pairing algorithm for initial ranging process in orthogonal frequency-division multiple access (OFDMA) uplink systems by exploiting multidimensional harmonic retrieval (HR) technique. Unlike most existing techniques that estimate each parameters independently and associate the ranging codes with the estimated parameters manually, the proposed method estimates multiple user’s ranging code, timing and frequency offsets simultaneously, and pair up all the multiple parameters automatically. The simulation results confirm that the proposed technique not only improves the ranging code detection capability and adjusts the acquisition range of the estimators but also increases the maximum number of resolvable users for given samples.


Introduction
In general orthogonal frequency-division multiple access (OFDMA) uplink systems, to maintain the orthogonality among multiple users, the signals from all active users should arrive at the base station (BS) synchronously [1]. In order to align the signals of multiple users, contentionbased random access procedure, called as initial ranging process, should be performed in the beginning in order to identify multiple users as well as to estimate the corresponding timing and/or frequency offsets for the adjustment/alignment [2,3]. So, code identification as well as multiuser timing estimation are the main tasks of the BS in contention-based synchronization processes.
On the contention-based synchronization process, the same time/frequency resources are shared by multiple users, so the multiple access interference (MAI) limits the performance of ranging detectors. Consequently, it is essential to deal with MAI of ranging subscriber stations (RSSs) to achieve the reliable initial ranging performance.
In the beginning, the existing ranging detection methods simply treat the MAI as noise, which results in *Correspondence: sungeun@gatech.edu School of Electrical and Computer Engineering, Georgia Institute of Technology, 777 Atlantic Drive NW, Atlanta, GA 30332-0250, USA the performance degradation as the number of RSSs increases due to the limitation by the amount of MAI [4][5][6]. To decrease the effect of the accumulated MAI on the ranging process, Ruan et al. [7] proposed a successive multiuser detection (SMUD) algorithm which detects the channel paths of active ranging signals and cancels their interference for further detection. In order to alleviate the high computational complexity of SMUD, the low-complexity method is proposed in [8], and in order to increase the resilience to MAI, a generalized likelihood ratio test (GLRT) approach is also exploited to derive a two-stage interference cancellation scheme [9]. However, it is shown that still SMUD and interference cancellation-based methods are vulnerable for MAI, and the code detection and parameter estimation performance are severely degraded as the number of active RSSs increases due to accumulated residual MAI [10].
In order to reduce the chance undergoing severe amount of MAI and exploiting the frequency selectivity of the fading channels, subchannel-based frame structure is proposed which allocate a smaller number of subcarriers to each ranging opportunity so that most of the RSSs are expected to transmit on disjoint sets of subcarriers with alleviated interference to each other [11][12][13][14][15]. http://asp.eurasipjournals.com/content/2015/1/1 With the improved parameter estimation capability, subspace decomposition-based approaches are exploited to accommodate the detection of multiple RSSs simultaneously [14,15]. However, in order to achieve the code detection and parameter estimation at the same time, the allowable parameter estimation range is strictly limited. In order to bind the detected ranging codes with the estimated timing or frequency offsets, the methods either sacrifice the estimation performance by inserting the code sequence in the estimated parameters [14] or perform exhaustive full search by correlating all sets of ranging codes and offsets at the BS [11,15]. Bacci et al. introduced a game-theoretic approach to derive an energy-efficient method for contention-based synchronization problem [16]. The ranging process is formulated as a noncooperative game for each terminal to determine the transmit power and detection strategy in a selfish way to maximize the detection probability with the minimum energy consumption.
The topic of multidimensional harmonic retrieval (HR) problems are encountered in a variety of signal processing applications including radar, sonar, seismology, communications, MIMO wireless channel, and the different types of schemes are utilized for the application [17][18][19][20][21]. Among these, principal-singular-vector utilization for modal analysis (PUMA) is used as the computationally attractive solution in HR [17]; however, it may not work properly when there are identical frequencies in one or more dimensions. On the other hand, the improved multidimensional folding (IMDF) scheme, which uses eigenvector instead of eigenvalue for estimating frequencies, can resolve identical frequency scenarios by introducing the randomness on the sample data [18,19]. These IMDF algorithms introduced more relaxed identifiability bound, i.e., increased number of resolvable frequencies for the given sample data [19]. In order to improve the performance, the approaches with the use of higher-order singular value decomposition (HOSVD) and/or structured least squares technique are introduced in [20,21]. Instead of using stacking matrix form, the HOSVD method captured the structure inherent in the received data at the expense of a high computational complexity [20,21]. It was shown in [20] that the general computational complexity of the tensor-based approach (HOSVD) is higher than that of the matrix based one (IMDF) even with the computationally efficient subspace estimation methods.
In order to improve the ranging process performance as well as to incorporate the superiority of HR technique with reasonable complexity, the preliminary design utilizing IMDF approach was first introduced in [22]. However, since the single snapshot IMDF approach [18] is simply exploited in [22], the automatic pairing property and multiple resource element utilization were not fully exploited in [22], which results in the limited number of resolvable RSSs as well as the limited detection and estimation performance improvement without proper two-dimensional code design.
Since the multiple snapshot approach is quite well aligned with the inherent OFDMA subchannel allocation concept, we propose an improved ranging technique using IMDF with finite snapshots [19] for the initial ranging process. It enables for all the detected ranging codes, the estimated timing offsets and frequency offsets of multiple users to be automatically paired, so extra ranging code and offset association process, which usually requires high computational complexity, can be avoided.
Even though our approach is based on the IMDF algorithm, but it is not a trivial application of the IMDF algorithm to the ranging process. On the top of the IMDF methodology, our novel contribution is twofold.
First, we properly formulate initial ranging signals into finite snapshots such that automatic pairing property can be fully exploited during the procedure, which brings the improved performance as well as the reduced pairing complexity. Moreover, the maximum number of resolvable RSSs that the BS can distinguish can be also increased by applying multidimensional folding (MDF) algorithm.
Next, with the virtue of the automatic paring property, two-dimensional ranging code is newly proposed. The new code index pair design, ranging code extraction and detection procedure, decision boundary for the circular code index pair, and redefined detection events using the code index pair with multiset theory are introduced. According to the new code index pair concept, the enlarged timing and/or frequency offset acquisition range is also derived.
The rest of this paper is organized as follows: Section 2 explains an OFDMA uplink system model in a contentionbased synchronization (initial ranging) process. Section 3 explains the proposed harmonic retrieval algorithm with multidimensional folding for extracting the code, timing offset, and frequency offset of active RSSs. In addition, the proposed two-dimensional code design and code detection/evaluation methodology is described. In Section 4, the statistical characteristics of the proposed algorithm, which utilizes IMDF method, is derived to show the identifiability for the number of resolvable RSSs and to confirm the acquisition range of the proposed algorithms on the timing offset and frequency offset. Section 5 evaluates the performance of the proposed method and investigates the comparison with the existing algorithms. Section 6 presents our conclusion.
Throughout this paper, upper (lower) boldface letters will be used for matrices (column vectors). A * , A T , A H , and A † denote the conjugate, transpose, Hermitian transpose, and pseudo inverse of A, respectively. We will use ⊗ for the Kronecker product, for the Khatri-Rao (column-wise Kronecker) product [23], I n for a n × n identify matrix, 0 m×n for an m×n zero matrix, [a] n for the nth element of a, and [A] m,n for the (m, n)th element of A, D(a) for a diagonal matrix with a as its diagonal. Figure 1 shows the brief signal flow on the initial ranging process. First of all, each RSS catches the downlink synchronization signal to detect its relative frame timing position and to estimate its absolute carrier frequency offset. Since every RSS can compare and compensate its frequency to the reference one from the BS, the coarse frequency synchronization can be completed before the initial ranging trial is started. Therefore, the residual frequency offset error due to the imperfect estimation is only remained for each RSS as shown in Figure 1.

Ranging process overview
On the other hand, the absolute frame timing of BS is unknown to each RSS since the propagation delay is unknown to each RSS (refer Figure 1). Therefore, each RSS only estimates its relative downlink frame timing position at first. Next, each RSS transmits its own initial ranging sequence, which is locked on the estimated frame timing, to the BS with the network access request. Once BS receives the sequences from multiple RSSs, the BS should identify multiple RSSs and estimate every RSS's timing offset and residual frequency offset independently. Finally, from the given estimate, the BS provides the timing and frequency feedback info to the identified RSSs with the network grant message. Consequently, it is crucial on the initial ranging process for BS to identify as many as RSSs and to estimate their timing and frequency offsets as accurate as possible.
The whole transceiver block diagram of the ranging process is shown in Figure 2, and the detailed procedure on how to construct the ranging sequence with the code and to perform the code detection and timing/frequency estimation is described in the following.

Resource allocation for ranging process
Let us consider an OFDMA uplink system employing N subcarriers. Among the whole subcarriers, virtual subcarriers are placed at both edges of the spectrum to prevent the spectrum aliasing. Except the virtual subcarriers, the useful subcarriers are grouped by multiple subchannels, and these subchannels are assigned to multiple users (subscribers) for transmission. Typically, each subchannel is divided into Q subbands composed of a set of V adjacent subcarriers [14,15,24]. When one transmission block consists of M consecutive OFDMA symbols, let us define a bunch of consecutive V subcarriers (one subband) and M OFDMA symbols as a tile. Now, among multiple subchannels, let us assume that each RSS only uses a set of R subchannels for ranging with M consecutive OFDMA symbols, and then the total number of subcarriers used for ranging is N R = RQV . After the fast Fourier transform (FFT) at the receiver, the received data at the qth tile of the ranging subchannel is written as: Since there is no big frequency offset difference between stations after the coarse frequency offset synchronization, the inter-carrier-interference (ICI) term caused by the residual frequency offset is ignored for the simplicity in (1) whereas the accumulated phase rotation term caused by frequency offset is still included.
In (1), the length of one OFDMA symbol is N T = N + N G with N G guard interval for ranging symbols, and i q is the starting subcarrier index for the qth tile, and θ k and k denote the timing and frequency offsets of the kth RSS, respectively. W q v,m is a complex additive white Gaussian noise (AWGN) with zero mean and variance σ 2 w . C k is the code sequence matrix for the kth RSS.
Here, let us assume that θ k ∈ [0, θ max ] and k ∈ [− max , max ] where θ max and max are defined as the maximum allowable timing offset and and absolute frequency offset of the system, respectively.

Code design
Let us design the code sequence for the kth RSS user as: Here, C T and C F are defined as the design parameters to determine the two-dimensional code index pair grid. Note that C T and C F can be properly chosen according to RSS's timing and frequency distributions. Let us http://asp.eurasipjournals.com/content/2015/1/1 define a k = [t k f k ] T as a two-tuple for the code index pair of the kth user where t k ∈ [0, C T ) is the code index at the timing offset domain and f k ∈ [0, C F ) is the one at the frequency offset domain, respectively. The tuple a k is chosen from the whole code index pair set A, i.e., a k ∈ A.
Since the subband size in one tile is designed to be smaller than the channel coherence bandwidth [14,15], H k [x] ≈ H k [x + 1] within the tile. Therefore, by averaging V subcarriers, the received output can be expressed as an average channel response as follows: where h q k is the average of the phase-rotated channel response for the kth user at the qth tile written as: In (3), the code sequence of the kth user, [C k ] v,m , is absorbed into the effective timing and frequency offsets as η k and ξ k , and this effective timing and frequency offsets of the kth user are given by:

MDF estimation using multiple tiles
In order to exploit automatic pairing property of HR, we formulate the received signal as a two-dimensional (2-D) mixture form and then apply MDF estimation method [18] for estimating timing and frequency offsets of RSSs as well as detecting RSSs codes simultaneously.

2-D mixture model
Let us recall the received symbols at the qth tile is X q ∈ C V ×M . Given (3), X q can be written as 2-D mixture matrix form as: where T ∈ C V ×K and F ∈ C M×K are the transpose of Vandermonde matrices with the common ratios e j2πη k and e j2πξ k for each column to represent the timing and frequency offsets of K RSSs, respectively.
For simplicity, let us define the (v, m)th element of the matrix X q as x q,v,m , i.e., x q,v,m = [X q ] v,m . Then, from the given Eq. 6, we can define the sample vector form x q ∈ C VM×1 , which is: where w q ∈ C VM×1 is the vector form of the noise matrix W q . Since the received signal is formed as undamped 2-D exponentials for each tile, we can apply MDF algorithm [18] for each tile at the initial ranging process.

Smoothing operation
As shown in (6) and (7), the observed symbols at each tile can be described as the special form of the multiple harmonic combinations. Based on this structure, we even apply smoothing operator S(x) in order to fully utilize this special characteristics of the observed symbols and improve the estimation performance.
Let us define a set of two-dimensional selection matrices where V 1 , V 2 and M 1 , M 2 are positive integers satisfying Using (7) and (8), let us define smoothing operator S(x) and the smoothed matrixẊ q,S for the sample vector as: Here, we show an explicit expression ofẊ q,S using x q,v,m in order to figure out the smoothing operator: . http://asp.eurasipjournals.com/content/2015/1/1 In the absence of noise, it is already verified in [18], Lemma 2 that (9) can be decomposed into harmonic matrices as follows: where: with , are the transpose of Vandermonde matrices with the common ratios e j2πη k and e j2πξ k . For the noisy case, the perturbation analysis can be applied to derive the theoretical MSEs of one realization of timing and frequency offset estimation [19], Equations (33) to (47). However, since timing offset θ k and frequency offset k are all unknown random variables in the ranging process, the theoretical performance analysis on the deterministic realization of η k and ξ k would be meaningless on the whole. Therefore, for the simplicity, the noiseless case is used on the derivation while the noisy case is revisited on the simulation experiment.
In (10), in order to further explore the data structure, we operate the backward smoothing as well as forward smoothing on the sample vector x q in (7). Let us define a backward sample vector y q = VM · x * q where n is an n×n backward permutation matrix. By using the property of harmonics, it is easily verified that: whereh q = D e * β h * q and e β = e jβ 1 e jβ 2 · · · e jβ K T with β k = (V − 1)η k + (M − 1)ξ k [19]. Since (13) has the same harmonic matrix form as (7), we can use both forward and backward sample vectors to estimate the phase. In the similar way like (9),Ẏ q,S can be also obtained by backward smoothing operation:

Multiple tile combination
Let us define the forward and backward channel matrices as: where H ∈ C K×Q ,H ∈ C K×Q . Whereas [22] uses each tile independently for the code detection and timing/frequency offset estimation, now all channel responses and signals from multiple tiles are stacked and utilized together to create just one augmented matrixX.
Therefore, it would be more robust for the noisy scenario due to multiple tile usage. After calculating the smoothed matricesẊ q,S andẎ q,S for each tile, we collect all the smoothed matrices in order to perform MDF algorithm [19]: whereX ∈ C V 1 M 1 ×2QV 2 M 2 and B ∈ C 2Q×K as follows:

Eigenvector-based estimatioñ
X is of rank K almost surely. Therefore, the singular value decomposition (SVD) ofX yields: where U s has K columns that together span the column space ofX. Since the same space is spanned by the columns of G 1 from (19), there exists K × K nonsingular matrix L −1 such that U s = G 1 L −1 . U s can be partitioned into two submatrices U 1 and U 2 (refer (46) and (47) in Appendix). Then, L can be obtained from the eigenvalue decomposition of U † 1 U 2 up to column scaling and permutation ambiguity, i.e., L sp = L where L sp is a practically obtained eigenvector including column scaling and permutation ambiguity, is a nonsingular diagonal column scaling matrix, and is a permutation matrix [18,19]. This part is maybe most computationally complex part of the whole calculation, and the major computational load complexity to acquire L sp is calculated as [19][20][21].
Once we obtain L sp , we can retrieve P ∈ C (V 1 −1)M 1 ×K up to column scaling and permutation ambiguity as: where P sp is a practically calculated matrix which consists of Khatri-Rao product of two Vandermonde matrices with are K columns up to column scaling and permutation ambiguity. Here, the original stacking matrix P is defined as: where α η k = e j2πη k and α ξ k = e j2πξ k , i.e., [P] In P, the phase rotation components caused by timing offset and frequency offset appear in the same column of the P. In other words, for a fixed kth column, e j2πη k and e j2πξ k appear simultaneously in the same column of P. Consequently, thanks to this combined structure, we can estimate K RSS's timing offset and frequency offset at the same time by dividing suitably chosen elements of the P.
The detailed procedure on how to obtain the above P sp and L sp from U s is described in Appendix.

Effective timing offset and frequency offset estimation
Even though P sp contains the column scaling ambiguity , it is immaterial for estimating the timing and frequency offset of each user. It is because the timing and frequency offset of each user is obtained by phase differential term. Therefore, the scaling effect on the matrix is diminished.
In addition, each timing offset with frequency offset is even automatically paired without heavy computational search for the pairing. The effective timing offset on the kth column can be estimated as: and the effective timing offset from the kth column can be obtained by: Im ln e j2πη k .
In the similar way, the effective frequency offset on the kth column can be estimated as: and the effective frequency offset from the kth column can be obtained by: Im ln e j2πξ k .
Note that still we have the the permutation ambiguity issue on the estimation, i.e., unknown . Since the user (column) identification along with timing and frequency estimation should be performed accordingly in the ranging process, therefore it is crucial to identify each column correctly and to remove the permutation ambiguity clearly. This can be done by the proposed code extraction and code detection process.

Code extraction
Recall the effective timing and frequency offsets are composed of the code index as well as actual timing and frequency offsets as shown in (4) and (5). Sinceη k andξ k are automatically paired, the paired ranging code, timing offset, and frequency offset of the kth RSS can be detected and estimated simultaneously by using the effective timing and frequency offsets.
Followed by (4) and (5), let us define the extracted code index pair b k as: where T = θ max 2 C T N . By plugging (4) and (5) in (28), b k is derived as follows: where θ max Since the timing offset θ k only has the positive values, by introducing the bias T on the extracted code index pair, the disturbance due to timing offset can be spread out around the original code index point t k . Consequently, by observing b k , we can retrieve the code index pair (t k , f k ) accordingly.

Code detection
Now, the new code detection methodology and the decision boundary/criterion for the circularly repeated code index pair are introduced using multiset theory. In addition, the code detection events are newly redefined for the multidimensional code, e.g., phantom code, collided code, missing code, false-alarmed code, and correctly detected code in this section. Figure 3 illustrates the example of code index pair set A, decision boundary, and decision region R k for C T = C F = 5. As discussed before, from the automatic pairing property of HR, it is possible to build up the two-dimensional code index pair.
From b k , the code index pairâ k can be extracted as follows: where R i is the decision region for the code index pair a i , ∀i ∈ [1, |A|]. Since the code index is circularly repeated per design parameter C T and C F , the modified maximum likelihood approach which reflects the repeated lattice http://asp.eurasipjournals.com/content/2015/1/1 structure can be applied for extracting the code index pair as follows: where κ(k) is the detected code index for the kth received tuple (kth extracted column), Z is the set of all integers, and x 2 2 is the square of the two-norm. Once all of K code index pair is detected, and then, it should be investigated from the detected code index pairs whether or not it contains the phantom code index pair which means the index pair is duplicated more than one time.
In order to redefine the code detection events using code tuple, let us assume A T and A R as the transmitted and the detected code index pair multisets (or bag), respectively, i.e., a k ∈ A T andâ k ∈ A R . In contrast to set concept, the members of the multiset are allowed to appear more than once. Here is an example to figure out the different events on the detection procedure: A R = {a 1 , a 1 , a 2 , a 2 , a 3 , a 6 } = â k .
This illustrates the case that six active RSSs are participated on ranging process, and two of the RSSs pick up the same code index pair a 1 and transmit (A T ), and the BS detects the six code index pairs from the received tuple (A R ). The multiset indicator function of a subset of a multiset A is: defined by: where Z + is the nonnegative integer set and μ(x) denotes the number of occurrences of x in the set. Now, let us take a look at different events on the code detection procedure. http://asp.eurasipjournals.com/content/2015/1/1

Phantom code
Let us define phantom code that one code index pair is observed more than once at A R . The index set of the phantom code is defined as: with κ(k) defined in (31). Since a 1 and a 2 are detected more than once in (34), here, G ph = {1, 2}. The phantom code can be detected at the receiver, so this code index can be opted out during the detection procedure.

Collided code
Collided code is defined such that one code index pair is utilized more than once at A T , which is the case when more than two RSSs choose the same code for the ranging. The index set of the collided code is defined as: From the example, G co is defined as G co = {1}. Typically, it is hard to detect the collided code at the receiver side, but with the virtue of automatic pairing property in HR, if each RSS has different timing offsets and/or frequency offsets, it has a chance to remove this collided code as phantom code at the receiver side (like the example). Therefore, HR can avoid the false detection of the collided code effectively.

Missing code
Let us redefine missing code such that the transmitted code index pair is not detected at the receiver side except collided code, which is denoted as: From the example, the missing code is shown as G m = {4, 5}.

False-alarmed code
False-alarmed code means that the wrong code index pair, which is never transmitted, is detected at the receiver side. There is no way to classify the false-alarmed code at the receiver, so it usually degrades the detection and estimation performance of the system. The false-alarmed code can be calculated as: Here, G f = {6}.

Correctly detected code
After excluding the phantom code, all the remained codes are treated as detected codes including the false-alarmed code. However, in order for the accurate evaluation of the code detection performance, let us newly define the correctly detected code which only counts the code truly transmitted and accurately received as: From the example, the correctly detected code is calcu- On the whole, the index set of detected codes is described as {3} while the index multiset of misdetected codes including both G co and G m is {1, 1, 2, 4, 5}. Finally, the index set of false-detected codes is defined as {6}. Note {1} appears twice on the misdetected multiset since A T already carries the two collided code.

Actual timing and frequency estimation
Among K extracted code index pairs, the phantom code G ph is firstly excluded, and only remaining code index the actual timing and frequency estimation. So, the timing and frequency offsets are calculated only for the κ(k * )th (4) and (25), the estimated actual timing offset for the κ(k * )th RSS from the k * th received tuple is given by: and it would be compared with θ k * if it is successfully identified without false-alarm (k * ∈ G d ) and the offset is accurate estimated. In the same way, the estimated actual frequency offset for the kth RSS is written by: and it would be also compared with k * if it is successfully identified without false-alarm (k * ∈ G d ).

Identifiability
When the estimator uses the timing offset dimension as a primary dimension, the maximum number of RSSs, which is almost surely uniquely identifiable, is derived as [19], Equation (31): Actually, (44) can be easily induced from (19) since the matrix size of G 1 andĜ 2 would determine the overall rank.
By adjusting the smoothing factors V 1 , V 2 , M 1 , and M 2 , the code detection and timing and frequency offset http://asp.eurasipjournals.com/content/2015/1/1 estimation performance can be improved as well as K uid can be increased.

Acquisition range
Basically, the transmitted code index pair a k * is twisted by the timing and frequency offsets. So, large timing and/or frequency offsets can cause the misdetection and/or false detection. Therefore, for the reliable code detection, it is essential to limit the allowable timing and/or frequency offset ranges. In principle, once b k * resides on the correct decision region, i.e., b k * ∈ R κ(k * ) , the timing and frequency offsets of the κ(k * )th RSS can be successfully recovered from the k * th received tuple, (refer Figure 3), and depending on the individual θ k * , ξ k * , and R κ(k * ) , the maximum allowable timing and frequency offset ranges, θ max and max , can be changed. Meanwhile, from (4) and (5), the guaranteed acquisition range not to make ambiguity at the code detection without any noise term can be derived as: where θ max and max stand for the maximum allowable timing and frequency offsets, respectively. So, from the inequality, it is available to adjust the maximum ranges of timing offset and frequency offset accordingly by enlarging one of the allowable ranges while shrinking another ranges at the same time. Remark that even though the acquisition range is derived as a lower bound by taking only a portion of the decision region R in (45), the maximum allowable ranges θ max and max are larger than the one introduced in MUSIC [15] and ESPRIT [14]. For example, with the same C T = C F = 5, the maximum ranges of MUSIC and ESPRIT were defined as θ max < N C T and max < N 2N T C F while the ranges of the proposed HR can be picked up as θ max < (whose range is enlarged by more than 10%) when we just treat the importance of both timing offset and frequency offset ranges equivalently with d min = √ 5 from the Figure 3.

Simulation results
The overall system model and simulation parameters are based on the IEEE 802.16m and 3GPP LTE environments in [2,3]. The total bandwidth is 10 MHz, and N = 1024. Let us assume that each ranging subchannel is composed of Q = 3 tiles with V = 6 consecutive subcarriers, while M = 5 OFDMA blocks are presented in ranging subchannels. In order to corporate with wide bandwidth, extended vehicular A (EVA) channel model is applied to evaluate the performance [25].
The timing and frequency offset ranges are designed for θ k ∈ [0, 114) and k ∈[ −0.02, 0.02] which correspond to the system more than 1 km cell radius considering the round-trip delay, e.g., dense small cell, and 135 km/h Doppler frequency at f c = 2.4 GHz. Note that all the simulation results reflect the ICI caused by residual frequency offset even though it is neglected in (1) for the simplicity.
The code design parameters C T and C F are set as C T = C F = 7 for the proposed HR algorithm.
To verify the superiority of the proposed algorithm in terms of the code detection and offset estimation performance, the ranging algorithms using MUSIC [15] and ESPRIT [14] are also compared under a common simulation setup. In addition, the ranging algorithm exploiting Zadoff-Chu sequence structure [2,26] is also compared as a reference under almost the same environment. Here, to evaluate the code detection performance solely, no collided code is assumed at the transmitter for all the algorithms, i.e., G co = {} like [14,15], which means all transmitted code index pair is exclusive to each other. The overall simulation parameters and frame structure are shown in Table 1, Table 2, and

Smoothing factor decision and identifiability comparison
Let us recall that the smoothing operator S(x) enlarges the observed data matrix X q ∈ C V ×M to the smoothed data matrixẊ q,S ∈ C V 1 M 1 ×V 2 M 2 by choosing appropriate smoothing factors V 1 , V 2 and M 1 , M 2 with the condition V 1 + V 2 = V + 1 and M 1 + M 2 = M + 1. Therefore, the size of smoothed data matrixẊ q,S depends on the smoothing factors. In fact, this smoothing factor could affect the code detection and offset estimation performance as well as the maximum number of distinguishable RSSs. As shown in (44), K uid is determined by (V 1 , V 2 ) and (M 1 , M 2 ) pair values, and Figure 5a shows the variation of this K uid according to V 1 and M 1 for given Q = 3, V = 6, and M = 5. First, we investigated most popularly obtained number in K uid distribution, which is K uid = 12 as shown in Figure 5a. Then, from the given four cases with K uid = 12, i.e., {(V 1 , M 1 ) |(4, 4), (5,3), (6,4), (5, 5) }, the smoothing factor combination which shows better performance on the whole is picked up, which is (V 1 , M 1 ) = (4, 4).
In addition, Figure 5b shows the trend on the maximum number of identifiable RSSs according to the number of subcarriers for fixed number of OFDM symbols M = 5. As seen in the figure, the maximum number of RSSs which MUSIC and/or ESPRIT can support is limited by M even though we allocate more subcarriers for the ranging whereas the number of RSSs for HR can keep increasing without limitation. This addresses that the proposed HR technique can be popularly exploited on the scenarios such that high volume of users and/or devices should be separated at the same time at the dense small cell. Figure 6 describes the misdetection probability for four code detection algorithms with respect to signal-to-noise ratio (SNR) with the number of active RSSs R · K = 12 and 24. The notation for each code detector and offset estimator is expressed in the Table 3.

Misdetection probability
As seen in the figures, the proposed HRCD shows good P miss performance. It is because automatic pairing property in HR improves the performance of code indices detection. Even though ESPRIT code detector (ECD) also exploits two independent estimates of the code indices to make a decision, the misdetection frequently occurs for ECD due to the discrepant code indices for two independent estimations. In addition, it is shown in Figure 6a,b that ZCCD has very poor performance regardless of SNR. It is because the orthogonality property of Zadoff-Chu sequence can sustain only with small number of RSSs. Different multipath fading channels as well as asynchronous timing/frequency offsets among multiple RSSs deteriorate the ZCCD detection performance dramatically.
Note that in Figure 6b, MUSIC and ESPRIT misdetection performance cannot be evaluated for R · K = 24 since it exceeds the maximum number of RSSs which MUSIC and/or ESPRIT can support, which is the limiting factor of MUSIC and ESPRIT algorithms. is very poor, we expect that the overall estimation performance can be also degraded by the false detection accordingly. Figures 8 and 9 show the error variance performances of the timing and frequency offset estimation given the correctly detected codes, respectively, i.e., Var

Timing and frequency offset estimation accuracy
The performance is evaluated with different numbers of active RSS users in the system. Note that the performance evaluation on the timing and frequency estimation is only performed for the corrected detected code index pairs, i.e., κ(k * ) ∈ G d . Unlike other methods, the proposed algorithm with HR can link the k * th received tuple to the code index κ(k * ) directly without any brute-force code index search effort due to automatic pairing property.
As seen in the Figure 8, HRTE outperforms other timing estimators over the whole SNR ranges regardless of the number of active users. Since it is crucial to obtain the timing misalignment information of each RSS on the initial ranging process, these curve confirm that the proposed HRTE can sustain the great performance even with high contention volume.
In Figure 9, it is shown that HRFE outperforms MFE and EFE over the whole range, and HRFE becomes better than ZCFE with SNR increase. Even though HRFE shows worse performance on low SNR ranges compared to ZCFE, but please remark that this estimation error variance is only evaluated for the successfully corrected code sets, i.e., Var ˆ κ(k * ) − k * κ(k * ) ∈ G d . Therefore, if we consider the poor false-alarm and missing probability performance of ZCCD shown in Figure 7 (more than 10% of falsealarm and missing), the overall ZCFE performance will be severely degraded on the whole.
On the whole, the performance of other algorithms undergo severe performance degradation because wrongly detected code index disturbs timing and frequency offset estimation. However, the proposed HRTE and HRFE algorithm more robust to the wrongly detected code effect.  Table 3 Code detection and offset estimation notation

Conclusions
First, the proposed method utilizing IMDF scheme supports the flexible subcarrier allocation usage. Basically, the proposed algorithm can be applied for the subcarrier allocation based on small tile and/or resource block structure whereas Zadoff-Chu sequence algorithm requires only the subband allocation, which should consist of consecutive subcarriers. Since the proposed algorithm also supports frequency noncontiguous combination of multiple tiles, actually it would even bring the performance improvement on the proposed scheme due to the enriched frequency diversity experience on the sample data. In addition, tile concept is compliant to the state-of-the-art wireless cellular systems, e.g., IEEE 802.16m and 3GPP LTE-Advanced [2,3], so this scheme can be adopted directly to the practical system without any painful remedy on the frame and resource allocation structure.
Next, the proposed method using IMDF scheme fully utilizes the resources to distinguish multiple RSSs as many as possible. In the practical scenario which only a few OFDMA symbols are utilized in the ranging process, the proposed HR algorithm can still increase K uid by adjusting the number of tiles Q and the number of subcarriers V. Note that the number of RSSs in MUSIC and ESPRIT is easily restricted by the limited number of OFDMA symbols. Consequently, it confirms that the HR algorithm can smartly utilize the given data matrix to distinguish many RSSs simultaneously.
From the simulation results, it is shown that the proposed method outperforms the other algorithms in terms of the code detection and offset estimation while supporting more number of RSSs. Especially, the newly designed and introduced two-dimensional code index pair with the aid of automatic pairing property enables for the algorithm to increase the minimum Euclidean distance between two code index pairs and improve the acquisition range and detection/estimation performance accordingly. Consequently, this HR method can be proper for the dense network ranging process such as small cell, Machine-type communications, device-to-device communications. http://asp.eurasipjournals.com/content/2015/1/1 where P = T (V 1 −1) F (M 1 ) and e η = e j2πη 1 , e j2πη 2 , . . . , e j2πη K T . Assuming that the elements of e η are distinct, L can be obtained from the eigenvalue decomposition of U † 1 U 2 up to column scaling and permutation ambiguity [18], i.e., the eigenvalue decomposition of U † 1 U 2 gives: where L sp = L and is a nonsingular diagonal column scaling matrix and is a permutation matrix. Once we obtain L sp , we can retrieve P up to column scaling and permutation ambiguity: