# Patent

## Title
System and method for initial ranging in wireless communication systems

## Abstract

A system and method for initial ranging in wireless communication systems is provided. A plurality of orthogonal frequency division multiplexing (OFDM) blocks are received by an OFDMA base station transceiver from a plurality of remote user devices in wireless communication with the base station. A ranging subchannel is extracted from the OFDM blocks. The number of active codes in the ranging subchannel is determined, active codes are identified, and carrier frequency offsets (CFOs) are estimated for each active code. Timing delays and power levels for each active code are then estimated. The estimated CFOs, timing delays, and power levels are broadcasted by the base station to the remote user devices, so that the user devices can utilize same to adjust transmission parameters to optimize power levels and synchronize communication with the base station.

## Background
STATEMENT OF GOVERNMENT INTERESTS
The present invention was made with government support under National Science Foundation Grant Nos. ANI-03-38807 and CNS-06-25637. Accordingly, the Government has certain rights to the present invention.


BACKGROUND OF THE INVENTION
1. Field of the Invention
The present invention relates to wireless communication systems, and more particularly, to a system and method for initial ranging in wireless communication systems.
2. Related Art
In the field of wireless communication technology, recent demands for high data rates have stimulated intense research activity in multicarrier modulation techniques. Particular interest has been directed to orthogonal frequency-division multiple-access (OFDMA) wireless communication systems, which have become part of the Institute of Electrical and Electronics (IEEE) 802.16 family of standards for wireless metropolitan area networks (WMANs). OFDMA systems allow for the transmission of digital information between a base station (BS) and a plurality of remote user devices, such as cellular telephones, wireless handheld computers, etc. OFDMA systems split the available bandwidth into smaller subchannels composed by a set of orthogonal subcarriers, which are assigned to different users to simultaneously communicate with the BS. The user signals are received at the BS as a series of orthogonal frequency division multiplexing (OFDM) blocks, from which information data are extracted by the BS. Unfortunately, OFDMA systems are extremely sensitive to timing errors and carrier frequency offsets (CFOs) that may occur between uplink signals and the local references at the BS. Timing errors give rise to interblock interference (IBI), while inaccurate compensation of the CFOs destroys orthogonality among subcarriers and produces interchannel interference (ICI) as well as multiple access interference (MAI).
In an attempt to alleviate these drawbacks of OFDMA systems, the IEEE 802.16 standards specify a synchronization procedure called Initial Ranging (IR), wherein users that intend to establish a link with the BS adjust their transmission parameters so that uplink signals arrive at the BS synchronously and with approximately the same power level. In its basic form, the IR process develops through the following steps. First, a remote user computes frequency and timing estimates on the basis of a downlink control channel. The estimated parameters are then used in the uplink phase, during which each remote user transmits a randomly chosen code over a ranging time-slot that comprises a specified number of adjacent OFDM blocks. However, as a consequence of the different users' positions within the cell, uplink signals arrive at the BS at different time instants. Furthermore, since the ranging time-slot is randomly selected, several users may collide over a same ranging subchannel. After separating colliding codes and extracting timing and power information, the BS broadcasts a response message indicating the detected codes and giving the necessary instructions for timing and power adjustment.
The main functions of the BS during the ranging process may be classified as multi-user code detection, and multi-user timing and power estimation. In existing IR techniques, a long pseudo-noise (PN) code is transmitted by each RSS over all available ranging subcarriers. Code detection and timing recovery is then accomplished on the basis of suitable correlations computed in either the frequency or time domains. These approaches require huge computational complexity since each correlation must be evaluated for each possible ranging code and hypothesized timing offset. Moreover, separation of users is achieved by exploiting only the correlation properties of the employed code set. In the presence of multipath distortions, however, ranging subcarriers are subject to different attenuations and phase shifts, thereby leading to a loss of the code orthogonality. This gives rise to MAI, which severely degrades the system performance.
Alternative solutions involve replacing the PN ranging codes with a set of modified generalized chirp-like (GCL) sequences and mitigating the effects of channel distortion through differential detection of the ranging signals. Unfortunately, this approach is still plagued by significant MAT and, furthermore, it cannot be applied to interleaved OFDMA systems since ranging subcarriers must be assigned adjacently. Some attempts at reducing the system complexity rely on the decomposition of the multi-parameter estimation problem into a number of successive steps. However, they are only suitable for flat channels and fail in the presence of multipath distortions. Still another approach involves a reduction of the computational burden, but this is achieved at the price of an increased overhead using a dedicated common ranging code.
A common drawback, among others, of existing IR methods is their poor performance in the presence of frequency selective channels. For example, in one existing method a signal design is proposed which is robust to multipath distortions, wherein ranging signals are divided into several groups and each group is transmitted over an exclusive set of subcarriers with a specific timing delay. This approach leads to a significant reduction of MAI as signals of different groups are perfectly separable in either the frequency or time domain. Additionally, in existing approaches, phase rotations induced by residual CFOs between the uplink signals and the BS frequency scale are ignored even though such rotations may become significant over ranging periods spanning several adjacent blocks. In such cases, the received ranging signals are no longer orthogonal and CFO compensation is necessary to avoid a serious degradation of the system performance. Unfortunately, frequency recovery schemes for OFDMA systems are normally derived under the assumption that uplink signals are transmitted over disjoint subcarriers and, as a result, cannot be applied to a scenario where multiple codes share the same subchannel.
Accordingly, what would be desirable, but has not yet been provided, is a system and method for initial ranging in wireless communication systems, which addresses the foregoing shortcomings of existing initial ranging approaches.
SUMMARY OF THE INVENTION
The present invention relates to a system and method for initial ranging in wireless communication systems. OFDM blocks are received by an OFDMA base station transceiver from a plurality of remote user devices. After passing the received signal into a discrete Fourier transform (DFT) unit, ranging subchannels are extracted from the received blocks, and a ranging subchannel is selected for processing. The number of active codes that are present in the selected subchannel is then determined. Once this information is obtained, CFOs are subsequently estimated together with the active codes using a multiple signal classification (MUSIC) algorithm. Timing delays and power levels are then estimated for each active code using maximum likelihood (ML) and least-squares (LS) estimation algorithms. The BS then broadcasts a response message indicating the detected codes and their corresponding instructions for adjusting synchronization parameters and power levels. A reduced complexity algorithm is also provided by the present invention for estimating timing delays and power levels with reduced computational complexity. In another embodiment of the present invention, CFOs, timing errors, and power levels are estimated using an Estimation of Signal Parameters via Rotational Invariance Techniques (EPSRIT) algorithm.



BRIEF DESCRIPTION OF THE DRAWINGS
The foregoing features of the invention will be apparent from the following Detailed Description of the Invention, taken in connection with the accompanying drawings, in which:
FIG. 1 is a flowchart showing processing steps according to the present invention for initial ranging in a wireless communication system;
FIGS. 2-5 are flowcharts showing, respectively, processing steps 18-24 of FIG. 1 in greater detail;
FIG. 6 is a flowchart showing processing steps according to the present invention for estimating power levels and timing delays during initial ranging with reduced computational complexity;
FIG. 7 is a block diagram showing a sample hardware and software/firmware arrangement for implementing the present invention;
FIGS. 8-14 are graphs showing performance results of the present invention during software simulation testing;
FIG. 15 is a flowchart showing another embodiment of the present invention, wherein CFOs and timing delays are estimated using an ESPRIT algorithm; and
FIGS. 16A-16C are graphs showing performance results of the embodiment of the present invention shown in FIG. 15, during software simulation testing.



DETAILED DESCRIPTION OF THE INVENTION
The present invention relates to a system and method for initial ranging in wireless communication systems, wherein a specified number of OFDM blocks are received by the BS from a plurality of remote user devices. The received signal is passed to a DFT unit which extracts ranging subchannels from the received blocks, and a ranging subchannel is selected for processing. After determining the number of active codes in the considered subchannel, CFOs are estimated for each active code using a MUSIC algorithm. Timing delays and power levels are then estimated for each active code using ML and LS estimation algorithms. The estimated CFOs, timing delays, and power levels are broadcast by the base station to the remote user devices, which can optimize their power levels and synchronize with the BS reference scales. In another embodiment of the present invention, CFOs, timing errors, and power levels are estimated using an ESPRIT algorithm.
As used herein, matrices and vectors are denoted by boldface letters, with IN being the identity matrix of order N. The expression A=diag{a(n); n=1, 2, . . . , N} denotes an N×N diagonal matrix with entries a(n) along its main diagonal, while B−1 is the inverse of a square matrix B. The expressions E {•}, (•)*, (•)T and (•)H are used for expectation, complex conjugation, transposition, and Hermitian transposition, respectively. The notation ∥•∥ represents the Euclidean norm of the enclosed vector, while |•| represents the modulus. The expression └c┘ indicates the largest integer smaller than or equal to c, and [•]k,l denotes the (k,l)th entry of the enclosed matrix.
The present invention is operable with OFDMA wireless communication systems, in which remote wireless communication devices (referred to herein as “ranging subscriber stations” (RSSs)) engage in an IR procedure with an OFDMA BS transceiver (e.g., a transceiver at a cellular tower location, a building, etc.). The equations/algorithms according to the present invention can be utilized in OFDMA systems where M denotes the number of OFDM blocks in a ranging time-slot and each ranging subchannel can be accessed by no more than M−1 RSSs. The latter are separated by means of orthogonal codes which are randomly selected from a set {c1,c2, . . . ,cM−1}, with ck[ck(0),ck(1), . . . ,ck(M−1)]T and |ck(m)|=1 for 0≦m≦M−1. Each OFDM block employs N subcarriers with frequency spacing Δf and indices in the set {0, 1, . . . , N−1}. R denotes the number of ranging subchannels, each of which is divided into Q subbands uniformly spaced over the signal bandwidth at a distance (N/Q)Δf from each other. A given subband is called a tile, and is composed by a set of V adjacent subcarriers. The subcarrier indices in the qth tile (q=0, 1, . . . , Q−1) of the rth subchannel (r=0, 1, . . . , R−1) are collected into a set Jq(r)={iq,v(r)}v=1V−1 with entries




  
    
      
        
          i
          
            q
            ,
            v
          
          
            (
            r
            )
          
        
        =
        
          
            qN
            Q
          
          +
          
            rN
            QR
          
          +
          
            v
            .
          
        
      
    
    
      
        (
        1
        )
      
    
  




The rth subchannel is thus composed of subcarriers with indices taken from J(r)=∪q=0Q−1Jq(r),, while a total of NR=QVR ranging subcarriers are available with indices in the set JR=∪r=0R−1J(r),. The remaining N−NR subcarriers are used for data transmission and are assigned to data subscriber stations (DSSs), which have successfully completed their IR process and are assumed to be perfectly aligned to the BS time and frequency scales.

FIG. 1 is a flowchart showing processing steps according to the present invention for initial ranging in a wireless communication system (e.g., an OFDMA system), indicated generally at 10. The steps shown in FIG. 1 can be programmed into any suitable OFDMA base station transceiver, and allow such a transceiver to conduct initial ranging processes with RSSs (e.g., cellular telephones, etc.) in accordance with a suitable communication protocol, including, but not limited to, the standard set forth in the Institute of Electrical and Electronics (IEEE) 802.16e standard.
Beginning in step 12, an OFDMA base station transceiver receives M OFDM blocks from a plurality of RSSs that intend to establish a link with the BS. In step 14, each OFDM block is processed by a discrete Fourier transform (DFT) unit to extract the QV subcarriers corresponding to the R ranging subchannels. As a general rule, different RSSs in a given subchannel employ different codes. In step 16, a ranging subchannel is selected and modeled in accordance with the present invention, as follows.
For the purposes of this description, we concentrate on the rth subchannel, where K(r)≦M−1 denotes the number of simultaneously active RSSs. To simplify the notation, the subchannel index (r) is not referred to herein. The waveform transmitted by the kth RSS propagates through a multipath channel of order L with impulse response hk=[hk(0),hk(1), . . . ,hk(L−1)]T. At the BS, the received samples are not synchronized with the local references. θk denotes the timing error of the signal received from the kth RSS expressed in sampling intervals, while εk denotes the corresponding frequency offset normalized to the subcarrier spacing.
We denote by Ym a QV-dimensional vector which stores the DFT outputs of the considered subchannel during the mth OFDM block. Since DSSs are perfectly aligned to the BS references, their signals do not contribute to Ym. In contrast, the presence of uncompensated CFOs destroys orthogonality among ranging signals and gives rise to ICI. Recalling that the ranging subchannels are interleaved in the frequency domain at a distance of NΔf/(QR) from each other, interference on Ym arising from subchannels other than the considered one can be neglected. Under the above assumptions, it is possible to write:




  
    
      
        
          Y
          m
        
        =
        
          
            
              ∑
              
                k
                =
                1
              
              K
            
            ⁢
            
              
                
                  c
                  k
                
                ⁡
                
                  (
                  m
                  )
                
              
              ⁢
              
                e
                
                  j
                  ⁢
                  
                    
                  
                  ⁢
                  m
                  ⁢
                  
                    
                  
                  ⁢
                  
                    ω
                    k
                  
                  ⁢
                  
                    N
                    T
                  
                
              
              ⁢
              
                A
                ⁡
                
                  (
                  
                    ω
                    k
                  
                  )
                
              
              ⁢
              
                
                  S
                  k
                
                ⁡
                
                  (
                  
                    θ
                    k
                  
                  )
                
              
            
          
          +
          
            n
            m
          
        
      
    
    
      
        (
        2
        )
      
    
  




where ωk=2πεk/N is the angular CFO of the kth RSS, NT=N+NG denotes the duration of the cyclically extended blocks and nm is a Gaussian vector with zero mean and covariance matrix σ2IQV. In addition, the following equation has been defined.

A(ω)=FW(ω)FH   (3)

where W(ω) is given by:

W(ω)=diag{ejnω;n=0,1, . . . ,N−1}  (4)

while F=[F0HF1H . . . FQ−1H]H, with each Fq denoting a V×N matrix with entries:





  
    
      
        
          
            
              [
              
                F
                q
              
              ]
            
            
              v
              ,
              n
            
          
          =
          
            
              
                1
                
                  N
                
              
              ⁢
              
                e
                
                  
                    -
                    j
                  
                  ⁢
                  
                    
                  
                  ⁢
                  2
                  ⁢
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  
                    i
                    
                      q
                      ,
                      v
                    
                  
                  ⁢
                  
                    n
                    /
                    N
                  
                
              
              ⁢
              0
            
            ≤
            v
            ≤
            
              V
              -
              1
            
          
        
        ,
        
          0
          ≤
          n
          ≤
          
            N
            -
            1.
          
        
      
    
    
      
        (
        5
        )
      
    
  




Vector Sk(θk) in Equation 2 can be partitioned into Q segments [SkT(θk,i0)SkT(θk,i1) . . . SkT(θk,iQ−1)]T each corresponding to a distinct tile. The qth segment Sk(θk,iq) has length V and elements given by:





  
    
      
        
          
            
              S
              k
            
            ⁡
            
              (
              
                
                  θ
                  k
                
                ,
                
                  i
                  
                    q
                    ,
                    v
                  
                
              
              )
            
          
          =
          
            
              e
              
                
                  -
                  j
                
                ⁢
                
                  
                    2
                    ⁢
                    
                      
                    
                    ⁢
                    π
                    ⁢
                    
                      
                    
                    ⁢
                    
                      θ
                      k
                    
                  
                  N
                
                ⁢
                
                  i
                  
                    q
                    ,
                    v
                  
                
              
            
            ⁢
            
              
                H
                k
              
              ⁡
              
                (
                
                  i
                  
                    q
                    ,
                    v
                  
                
                )
              
            
          
        
        ,
        
          0
          ≤
          v
          ≤
          
            V
            -
            1
          
        
      
    
    
      
        (
        6
        )
      
    
  




where Hk(θk,iq,v) is the kth channel frequency response over the iq,vth subcarrier, which can be written as:





  
    
      
        
          
            H
            k
          
          ⁡
          
            (
            
              i
              
                q
                ,
                v
              
            
            )
          
        
        =
        
          
            ∑
            
              l
              =
              0
            
            
              L
              -
              1
            
          
          ⁢
          
            
              
                h
                k
              
              ⁡
              
                (
                l
                )
              
            
            ⁢
            
              
                e
                
                  
                    -
                    j
                  
                  ⁢
                  
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      l
                    
                    N
                  
                  ⁢
                  
                    i
                    
                      q
                      ,
                      v
                    
                  
                
              
              .
            
          
        
      
    
    
      
        (
        7
        )
      
    
  




The above equation indicates that θk appears only as a phase shift across the DFT outputs. The reason for this is that the CP duration is longer than the maximum expected propagation delay, thereby making the ranging signals quasi-synchronous at the BS. The power of the kth RSS over the ranging subcarriers is defined as:





  
    
      
        
          P
          k
        
        =
        
          
            1
            QV
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  0
                
                
                  V
                  -
                  1
                
              
              ⁢
              
                
                  
                    
                    
                      
                        S
                        k
                      
                      ⁡
                      
                        (
                        
                          
                            θ
                            k
                          
                          ,
                          
                            i
                            
                              q
                              ,
                              v
                            
                          
                        
                        )
                      
                    
                    
                  
                  2
                
                .
              
            
          
        
      
    
    
      
        (
        8
        )
      
    
  




In step 18, the number of currently active codes in the selected subchannel is calculated, as discussed below in greater detail in connection with FIG. 2. Then, in step 20, CFOs are estimated for each active code, as discussed below in greater detail in connection with FIG. 3. In step 22, after CFOs have been estimated, timing delays/errors for each active code are estimated. This process is discussed below in greater detail in connection with FIG. 4. In step 24, power levels for each active code are estimated, as discussed in connection with FIG. 5.
In step 26, a determination is made as to whether all subchannels have been processed. If a negative determination is made, control returns to step 16, so that an additional subchannel can be processed as described above to provide estimates of CFOs, timing delays, and power levels for the corresponding RSSs. If a positive determination is made, step 28 occurs.
After CFOs, timing delays, and power levels have been estimated for each RSS, this information is broadcast in step 28 to the remote users (i.e., the RSSs), so that each user device can adjust its transmission parameters and align its own uplink signal to the BS reference scales. Specifically, users that intend to access the network for the first time compute initial frequency and timing estimates on the basis of a downlink control signal broadcast by the BS. The estimated parameters are then employed by each RSS as synchronization references for the uplink ranging transmission. This means that, during the IR process, the CFOs are induced only by Doppler shifts and/or estimation errors and, as a consequence, they will be quite small. For example, consider the IEEE 802.16 standard for WMANs with subcarrier spacing Δf=11.16 kHz and carrier frequency f0=5 GHz. In the case of perfect downlink synchronization, the maximum CFO in the uplink is 2f0νk/c, where νk denotes the terminal speed while c=3×108 m/s is the speed of light. Letting νk=120 km/h yields εk≦0.10, which means that the normalized CFO lies within 10% of the subcarrier spacing. Timing errors depend on the distances of the RSSs from the BS and their maximum value θmax corresponds to the round trip propagation delay for a user located at the cell boundary. This parameter is known and is given by θmax=2R/(cTs), where R is the cell radius and Ts=1/(NΔf) is the sampling period. A simple way to counteract the effects of timing errors relies on the use of sufficiently long CPs comprising NG≧θmax+L sampling intervals. This leads to a quasi-synchronous network in which timing errors do not produce any IBI and result only in phase shifts at the output of the receive DFT unit. The CP of data blocks should be made just greater than the channel length to minimize unnecessary overhead. It follows that accurate timing estimation must be accomplished during the ranging period in order to align all active users with the BS time-slot and to avoid IBI over the data section of the frame.
FIG. 2 is a flowchart showing step 18 of FIG. 1 in greater detail. As mentioned above, step 18 of FIG. 1 allows for the calculation of the number of active codes in the considered (modeled) ranging subchannel. This information is utilized to derive CFO estimates, as well as estimates of timing delays and power levels of the corresponding active codes.
In step 30, the number of received codes in the considered subchannel is determined by evaluating the correlation matrix and applying an eigendecomposition (EVD) procedure. It is assumed that the CFOs are confined within a small fraction of the subcarrier spacing. In such a case, the demodulated signals incur negligible phase rotations over one OFDM block, and W(ω) in Equation 4 above can be approximated by IN. Substituting Equation 3 into Equation 2, and observing that FFH=IQV, yields:




  
    
      
        
          Y
          m
        
        ≈
        
          
            
              ∑
              
                k
                =
                J
              
              K
            
            ⁢
            
              
                
                  c
                  k
                
                ⁡
                
                  (
                  m
                  )
                
              
              ⁢
              
                e
                
                  j
                  ⁢
                  
                    
                  
                  ⁢
                  m
                  ⁢
                  
                    
                  
                  ⁢
                  
                    ω
                    k
                  
                  ⁢
                  
                    N
                    T
                  
                
              
              ⁢
              
                
                  S
                  k
                
                ⁡
                
                  (
                  
                    θ
                    k
                  
                  )
                
              
            
          
          +
          
            n
            m
          
        
      
    
    
      
        (
        9
        )
      
    
  




from which it follows that each CFO results only in a phase shift between contiguous OFDM blocks. The DFT output over the iq,vth subcarrier of the mth block is thus given by:





  
    
      
        
          
            Y
            m
          
          ⁡
          
            (
            
              i
              
                q
                ,
                v
              
            
            )
          
        
        ≈
        
          
            
              ∑
              
                k
                =
                1
              
              K
            
            ⁢
            
              
                
                  c
                  k
                
                ⁡
                
                  (
                  m
                  )
                
              
              ⁢
              
                e
                
                  j
                  ⁢
                  
                    
                  
                  ⁢
                  m
                  ⁢
                  
                    
                  
                  ⁢
                  
                    ω
                    k
                  
                  ⁢
                  
                    N
                    T
                  
                
              
              ⁢
              
                
                  S
                  k
                
                ⁡
                
                  (
                  
                    
                      θ
                      k
                    
                    ,
                    
                      i
                      
                        q
                        ,
                        v
                      
                    
                  
                  )
                
              
            
          
          +
          
            
              n
              m
            
            ⁡
            
              (
              
                i
                
                  q
                  ,
                  v
                
              
              )
            
          
        
      
    
    
      
        (
        10
        )
      
    
  




The vector Y(iq,v)=[Y0(iq,v),Y1(iq,v), . . . ,YM−1(iq,v)]T is then defined for collecting the iq,vth DFT output across all M ranging blocks. Then, from Equation 10, it follows that:





  
    
      
        
          Y
          ⁡
          
            (
            
              i
              
                q
                ,
                v
              
            
            )
          
        
        =
        
          
            
              ∑
              
                k
                =
                1
              
              K
            
            ⁢
            
              
                
                  S
                  k
                
                ⁡
                
                  (
                  
                    
                      θ
                      k
                    
                    ,
                    
                      i
                      
                        q
                        ,
                        v
                      
                    
                  
                  )
                
              
              ⁢
              
                Γ
                ⁡
                
                  (
                  
                    ω
                    k
                  
                  )
                
              
              ⁢
              
                c
                k
              
            
          
          +
          
            n
            ⁡
            
              (
              
                i
                
                  q
                  ,
                  v
                
              
              )
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




where Γ(ωk)=diag{ejmωkNT;m=0,1, . . . ,M−1} while n(iq,v)=[n0(iq,v),n1(iq,v), . . . ,nM−1(iq,v)]T is Gaussian distributed with zero-mean and covariance matrix σ2IM.

From Equation 11, it is seen that Y(iq,v) is a superposition of frequency-rotated codes {Γ(ωk)ck}k=1K embedded in white Gaussian noise and with random amplitudes {Sk(θk,iq,v)}k=1K. This model has the same structure as measurements of multiple uncorrelated sources from an array of sensors. EVD can then be performed on the correlation matrix of Y(iq,v), whereupon the active codes and their corresponding CFOs can subsequently be determined.
Estimation of the number of active codes over the ranging subchannel can be accomplished in step 30 through the EVD of the correlation matrix RY=E{Y(iq,v)YH(iq,v)}. Let λ1≧λ2≧ . . . ≧λM be the eigenvalues of RY arranged in non-increasing order. Then, from Equation 11, it follows that the M−K smallest eigenvalues are equal to the noise power σ2 while the largest K eigenvalues are related to σ2 and also to the average signal power. Mathematically, the following Equation applies:




  
    
      
        
          λ
          k
        
        =
        
          {
          
            
              
                
                  
                    
                      P
                      _
                    
                    k
                  
                  +
                  
                    σ
                    2
                  
                
              
              
                
                  
                    for
                    ⁢
                    
                      
                    
                    ⁢
                    1
                  
                  ≤
                  k
                  ≤
                  K
                
              
            
            
              
                
                  σ
                  2
                
              
              
                
                  
                    
                      for
                      ⁢
                      
                        
                      
                      ⁢
                      K
                    
                    +
                    1
                  
                  ≤
                  k
                  ≤
                  M
                
              
            
          
        
      
    
    
      
        (
        12
        )
      
    
  




where Pk are strictly positive quantities related to the signal power. If the quantities Pk are much greater than σ2, estimating K is not problematic since the eigenvalues exhibit a clear drop in passing from k=K to k−K+1. In practice, however, the drop cannot be easily detected since users starting a ranging procedure are required to keep the transmission power at a relatively low level in order to limit interference to DSSs. A further difficulty is that RY is not available at the receiver and must be replaced with some suitable estimate. One approach to solve these difficulties is to use the sample correlation matrix:





  
    
      
        
          
            R
            ~
          
          Y
        
        =
        
          
            1
            QV
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  0
                
                
                  V
                  -
                  1
                
              
              ⁢
              
                
                  Y
                  ⁡
                  
                    (
                    
                      i
                      
                        q
                        ,
                        v
                      
                    
                    )
                  
                
                ⁢
                
                  
                    Y
                    H
                  
                  ⁡
                  
                    (
                    
                      i
                      
                        q
                        ,
                        v
                      
                    
                    )
                  
                
              
            
          
        
      
    
    
      
        (
        13
        )
      
    
  




which provides an unbiased and consistent estimate of RY. An alternative solution is based on the forward-backward (FB) principle. In such a case RY is replaced by:





  
    
      
        
          
            R
            ^
          
          Y
        
        =
        
          
            1
            2
          
          ⁢
          
            (
            
              
                
                  R
                  ~
                
                Y
              
              +
              
                J
                ⁢
                
                  
                    R
                    ~
                  
                  Y
                  T
                
                ⁢
                J
              
            
            )
          
        
      
    
    
      
        (
        14
        )
      
    
  




where {tilde over (R)}Y is still given in Equation 13 while J is the exchange matrix with 1's on the anti-diagonal and 0's elsewhere.

In step 32, after performing EVD on {circumflex over (R)}Y and arranging the eigenvalues {circumflex over (λ)}1≧{circumflex over (λ)}2≧ . . . {circumflex over (λ)}M in non-increasing order, an estimate {circumflex over (K)} of the number of active codes is found through information-theoretic criteria. Two examples of such criteria are represented by the Akaike and Minimum Description Length (MDL) criteria. The MDL approach looks for the minimum of the following objective function:




  
    
      
        
          F
          ⁡
          
            (
            
              K
              ~
            
            )
          
        
        =
        
          
            
              1
              2
            
            ⁢
            
              
                K
                ~
              
              (
              
                
                  2
                  ⁢
                  M
                
                -
                
                  K
                  ~
                
              
              )
            
            ⁢
            
              ln
              ⁡
              
                (
                QV
                )
              
            
          
          -
          
            
              QV
              (
              
                M
                -
                
                  K
                  ~
                
              
              )
            
            ⁢
            ln
            ⁢
            
              
            
            ⁢
            
              ρ
              ⁡
              
                (
                
                  K
                  ~
                
                )
              
            
          
        
      
    
    
      
        (
        15
        )
      
    
  




where ρ({tilde over (K)}) is the ratio between the geometric and arithmetic mean of {{circumflex over (λ)}{tilde over (K)}+1,{circumflex over (λ)}{tilde over (K)}+2, . . . ,{circumflex over (λ)}M}, i.e.,





  
    
      
        
          ρ
          ⁡
          
            (
            
              K
              ~
            
            )
          
        
        =
        
          
            
              (
              
                
                  ∏
                  
                    i
                    =
                    
                      
                        K
                        ~
                      
                      +
                      1
                    
                  
                  M
                
                ⁢
                
                  
                
                ⁢
                
                  
                    λ
                    ^
                  
                  i
                
              
              )
            
            
              1
              
                M
                -
                
                  K
                  ~
                
              
            
          
          
            
              1
              
                M
                -
                
                  K
                  ~
                
              
            
            ⁢
            
              
                ∑
                
                  i
                  =
                  
                    
                      K
                      ~
                    
                    +
                    1
                  
                
                M
              
              ⁢
              
                
                  λ
                  ^
                
                i
              
            
          
        
      
    
    
      
        (
        16
        )
      
    
  




and {tilde over (K)} denotes a tentative value of K and varies from 1 to M−1.

In step 34, once {circumflex over (K)} has been computed, from Equation 12 an estimate of the noise power σ2 can be computed as:




  
    
      
        
          
            σ
            ^
          
          2
        
        =
        
          
            1
            
              M
              -
              
                K
                ^
              
            
          
          ⁢
          
            
              ∑
              
                i
                =
                
                  
                    K
                    ^
                  
                  +
                  1
                
              
              M
            
            ⁢
            
              
                
                  λ
                  ^
                
                i
              
              .
            
          
        
      
    
    
      
        (
        17
        )
      
    
  




FIG. 3 is a flowchart showing processing step 20 of FIG. 1 in greater detail. Processing step 20 allows for the estimation of CFOs for a given subchannel, as well as the identification of active codes in the subchannel. Beginning in step 36, the observation space is decomposed into a signal subspace Sσ which is spanned by the rotated codes {Γ(ωk)ck}k=1K plus a noise subspace Sv. As is known, Sv is the orthogonal complement of Sσ so that each vector of Sσ is orthogonal to any other vector in Sv. To proceed further, it is assumed that the number of active codes has been correctly estimated (i.e., {circumflex over (K)}≈K), and {û1,û2, . . . ,ûM} denotes the eigenvectors of {circumflex over (R)}y corresponding to the ordered eigenvalues {circumflex over (λ)}1≧{circumflex over (λ)}2≧ . . . ≧{circumflex over (λ)}M.
In step 38, CFOs are estimated using a MUSIC approach The MUSIC algorithm relies on the fact that the eigenvectors {ûK+1,ûK+2, . . . ,ûM} associated with the M−K smallest eigenvalues of {circumflex over (R)}Y form an estimated basis of the noise subspace Sv and, accordingly, they are approximately orthogonal to all vectors in the signal space Sσ. An estimate of ωk is thus obtained by minimizing the projection of Γ({tilde over (ω)})ck onto the subspace spanned by the columns of Ûn=[ûK+1ûK+2 . . . ûM]. This leads to the following estimation algorithm:




  
    
      
        
          
            
              ω
              ^
            
            k
          
          =
          
            
              
                arg
                ⁢
                
                  
                
                ⁢
                max
              
              
                ω
                ~
              
            
            ⁢
            
              {
              
                
                  Ψ
                  k
                
                ⁡
                
                  (
                  
                    ω
                    ~
                  
                  )
                
              
              }
            
          
        
        ⁢
        
          

        
        ⁢
        
          with
          ⁢
          
            :
          
        
      
    
    
      
        (
        18
        )
      
    
  
  
    
      
        
          
            Ψ
            k
          
          ⁡
          
            (
            
              ω
              ~
            
            )
          
        
        =
        
          1
          
            
              
              
                
                  
                    U
                    ^
                  
                  n
                  H
                
                ⁢
                
                  Γ
                  ⁡
                  
                    (
                    
                      ω
                      ~
                    
                    )
                  
                
                ⁢
                
                  c
                  k
                
              
              
            
            2
          
        
      
    
    
      
        (
        19
        )
      
    
  




In principle, CFO recovery must only be accomplished for the active codes. However, since at this stage the BS has no knowledge of which codes have been transmitted over the considered subchannel, frequency estimates {{circumflex over (ω)}1,{circumflex over (ω)}2, . . . ,{circumflex over (ω)}M−1} must be evaluated for the complete set {c1,c2, . . . ,cM−1}.
In step 40, active codes are identified. The present invention searches for the K largest values in the set {Ψk({circumflex over (ω)}k)}k=1M−1, say {Ψuk({circumflex over (ω)}uk)}k=1K, and declares as active the corresponding codes {cuk}k=1K. The CFO estimates are eventually found as {circumflex over (ω)}u=[{circumflex over (ω)}u1,{circumflex over (ω)}u2, . . . ,{circumflex over (ω)}uK]T. The CFO estimation based on the MUSIC algorithm is referred to herein as the MUSIC-based frequency estimator (MFE), while the active code identification algorithm is referred to herein as the MUSIC-based code detector (MCD).
A key issue is the maximum CFO value that can be handled by MFE. To simplify the analysis, it is assumed that the ranging codes belong to the Fourier basis of order M and are given by:

ck(m)=ej2πmkiM, 0≦m≦M−1   (20)

for k=1, 2, . . . ,M−1. Then, identifiability of the rotated codes {Γ(ωk)ck}k=1K is guaranteed as long as the CFOs are smaller than N/(2MNT) in magnitude. Let εk=N/(2MNT) and εk+1=−N/(2MNT). In such a case, the received codes Γ(ωk)ck and Γ(ωk+1)ck+1 are identical and there is no way to distinguish among them. This fact limits the acquisition range of MFE to the interval |εk|<N/(2MNT). Although this result has been derived for the Fourier basis specified in Equation 20, it may also apply to other code designs.

The computational load (complexity) of MFE can be assessed as follows. Assume that the number of active codes has already been found and is available at the BS. Then, computing {circumflex over (R)}Y in Equation 14 requires approximately O(M2QV) operations, while the complexity for evaluating the metric Ψk({tilde over (ω)}) in Equation 19 is O(M2(M−K)) for each {tilde over (ω)}. Denoting by Nω the number of candidate values {tilde over (ω)} and observing that the maximization in Equation 18 must be performed over the entire code set {c1,c2, . . . ,cM−1}, it follows that a total of O(NωM2(M−K)) operations are required to find {circumflex over (ω)}u. The overall complexity of MFE is thus on the order of O(M2QV+NωM2(M−K)). This figure does not consider the EVD on {circumflex over (R)}Y since this operation is comparatively easier than the maximization step in Equation 18.
FIG. 4 is a flowchart showing step 22 of FIG. 1 in greater detail. Step 22 allows for the estimation of timing delays of all active RSSs. After code detection and CFO recovery, the BS must acquire information about the timing delays and power levels of all ranging signals. Timing delay estimation is achieved in steps 42 and 44 through the use of ML and LS estimation algorithms. An initial assumption is made that {circumflex over (K)}=K. To simplify the notation, the indices {uk}k=1K of the detected codes are relabeled according to the map uk→k.
First, Equation 2 is reformulated in a more compact form. For this purpose, matrices Dm(ωk)=ck(m)ejmωkNTA(ωk) are defined and collect the CFOs and timing errors in two K-dimensional vectors ω=[ω1,ω2, . . . ,ωK]T and θ=[θ1,θ22, . . . ,θK]T. Then, Ym is rewritten as:

Ym=Dm(ω)S(θ)+nm   (21)

with S(θ)=[S1T(θ1)S2T(θ2) . . . SKT(θK)]T and Dm(ω)=[Dm(ω1)Dm(ω2) . . . Dm(ω2)]. Timing recovery is then accomplished by means of a two step procedure in which an estimate of S(θ) is computed first and then is exploited to obtain the timing errors.

In step 42, it is assumed that the accuracy of the CFO estimates obtained in the first stage of the ranging process is such that {tilde over (ω)}=ω. Then, from Equation 21, the ML estimate of S(θ) based on the observed vectors {Ym}m=0m−1 is found to be:




  
    
      
        
          
            S
            ^
          
          ⁡
          
            (
            θ
            )
          
        
        =
        
          
            
              T
              
                -
                1
              
            
            ⁡
            
              (
              
                ω
                ^
              
              )
            
          
          ⁢
          
            
              ∑
              
                m
                =
                0
              
              
                M
                -
                1
              
            
            ⁢
            
              
                
                  D
                  m
                  H
                
                ⁡
                
                  (
                  
                    ω
                    ^
                  
                  )
                
              
              ⁢
              
                Y
                m
              
            
          
        
      
    
    
      
        (
        22
        )
      
    
  




where T({circumflex over (ω)}) is the following matrix:





  
    
      
        
          T
          ⁡
          
            (
            
              ω
              ^
            
            )
          
        
        =
        
          
            ∑
            
              m
              =
              0
            
            
              M
              -
              1
            
          
          ⁢
          
            
              
                D
                m
                H
              
              ⁡
              
                (
                
                  ω
                  ^
                
                )
              
            
            ⁢
            
              
                
                  D
                  m
                
                ⁡
                
                  (
                  
                    ω
                    ^
                  
                  )
                
              
              .
            
          
        
      
    
    
      
        (
        23
        )
      
    
  




Substituting Equation 21 into Equation 22 yields:





  
    
      
        
          
            S
            ^
          
          ⁡
          
            (
            θ
            )
          
        
        =
        
          
            S
            ⁡
            
              (
              θ
              )
            
          
          +
          
            
              
                T
                
                  -
                  1
                
              
              ⁡
              
                (
                
                  ω
                  ^
                
                )
              
            
            ⁢
            
              
                ∑
                
                  m
                  =
                  0
                
                
                  M
                  -
                  1
                
              
              ⁢
              
                
                  
                    D
                    m
                    H
                  
                  ⁡
                  
                    (
                    
                      ω
                      ^
                    
                    )
                  
                
                ⁢
                
                  n
                  m
                
              
            
          
        
      
    
    
      
        (
        24
        )
      
    
  




from which it follows that Ŝ(θ) is an unbiased estimate of S(θ) and can be similarly partitioned as [Ŝ1(θ1)Ŝ2(θ2) . . . ŜK(θK)]T. The kth subvector Ŝk(θk) is composed of Q segments {Ŝk(θk,iq)}q=0Q−1 each corresponding to a specific tile in the considered ranging subchannel. From Equation 24, it can be seen that the vth entry of Ŝk(θk,iq), say Ŝk(θk,iq,v), is an unbiased estimate of the quantity Sk(θk,iq,v) defined in Equation 6. Hence, from Equations 6 and 7 it follows that:





  
    
      
        
          
            
              S
              ^
            
            k
          
          ⁡
          
            (
            
              i
              
                q
                ,
                v
              
            
            )
          
        
        =
        
          
            
              e
              
                
                  -
                  j
                
                ⁢
                
                  
                    2
                    ⁢
                    
                      
                    
                    ⁢
                    π
                    ⁢
                    
                      
                    
                    ⁢
                    
                      θ
                      k
                    
                  
                  N
                
                ⁢
                
                  i
                  
                    q
                    ,
                    v
                  
                
              
            
            ⁢
            
              
                ∑
                
                  l
                  =
                  0
                
                
                  L
                  -
                  1
                
              
              ⁢
              
                
                  
                    h
                    k
                  
                  ⁡
                  
                    (
                    l
                    )
                  
                
                ⁢
                
                  e
                  
                    
                      -
                      j
                    
                    ⁢
                    
                      
                        2
                        ⁢
                        
                          
                        
                        ⁢
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        n
                      
                      N
                    
                    ⁢
                    
                      i
                      
                        q
                        ,
                        v
                      
                    
                  
                
              
            
          
          +
          
            
              ξ
              k
            
            ⁡
            
              (
              
                i
                
                  q
                  ,
                  v
                
              
              )
            
          
        
      
    
    
      
        (
        25
        )
      
    
  




where ξk(iq,v) is a zero-mean disturbance term. The quantities {Ŝk(θk,iq,v)}q=0Q−1 are collected into a Q-dimensional vector Ŝk(θk,v) and Φ(θk,v)=diag{e−j2πθkiq,v/N;0≦q≦Q−1} is defined. Then, omitting for simplicity the functional dependence of Ŝk(θk,v) on θk, from Equation 25 it follows that

Ŝk(v)=Φ(θk,v)F(v)hk+ξk(v)   (26)

where ξk(v)=[ξk(i0,v),ξk(i1,v), . . . ,ξk(iQ−1,v)]T while F(v) is a matrix of dimensions Q×L with entries:





  
    
      
        
          
            
              [
              
                F
                ⁡
                
                  (
                  v
                  )
                
              
              ]
            
            
              q
              ,
              l
            
          
          =
          
            e
            
              
                -
                j
              
              ⁢
              
                
                  2
                  ⁢
                  
                    
                  
                  ⁢
                  π
                
                N
              
              ⁢
              
                li
                
                  q
                  ,
                  v
                
              
            
          
        
        ,
        
          0
          ≤
          q
          ≤
          
            Q
            -
            
              1
              ⁢
              
                
              
              ⁢
              and
              ⁢
              
                
              
              ⁢
              0
            
          
          ≤
          l
          ≤
          
            L
            -
            1.
          
        
      
    
    
      
        (
        27
        )
      
    
  




Equation 26 indicates that, apart from the disturbance vector ξk(v), Ŝk(v) is only contributed by the kth RSS, meaning that ranging signals have been successfully decoupled at the BS. Vectors {Ŝk(v);v=0,1, . . . ,V−1} can thus be used to obtain LS estimates of (θk,hk) individually for each RSS. This amounts to looking for the minimum of the following objective function:





  
    
      
        
          
            Λ
            k
          
          ⁡
          
            (
            
              
                
                  θ
                  ~
                
                k
              
              ,
              
                
                  h
                  ~
                
                k
              
            
            )
          
        
        =
        
          
            ∑
            
              v
              =
              0
            
            
              V
              -
              1
            
          
          ⁢
          
            
              
              
                
                  
                    
                      S
                      ~
                    
                    k
                  
                  ⁡
                  
                    (
                    v
                    )
                  
                
                -
                
                  
                    Φ
                    (
                    
                      
                        
                          θ
                          ~
                        
                        k
                      
                      ,
                      v
                    
                    )
                  
                  ⁢
                  
                    F
                    ⁡
                    
                      (
                      v
                      )
                    
                  
                  ⁢
                  
                    
                      h
                      ~
                    
                    k
                  
                
              
              
            
            2
          
        
      
    
    
      
        (
        28
        )
      
    
  




with respect to {tilde over (0)}k and {tilde over (h)}k. For a given {tilde over (θ)}k, the minimum occurs at;





  
    
      
        
          
            h
            ^
          
          k
        
        =
        
          
            1
            QV
          
          ⁢
          
            
              ∑
              
                v
                =
                0
              
              
                V
                -
                1
              
            
            ⁢
            
              
                
                  F
                  H
                
                ⁡
                
                  (
                  v
                  )
                
              
              ⁢
              
                
                  Φ
                  H
                
                (
                
                  
                    
                      θ
                      ^
                    
                    k
                  
                  ,
                  v
                
                )
              
              ⁢
              
                
                  
                    S
                    ^
                  
                  k
                
                ⁡
                
                  (
                  v
                  )
                
              
            
          
        
      
    
    
      
        (
        29
        )
      
    
  




where the identity FH(v)F(v)=Q·IL has been used, which holds true as long as L≦Q. Then, substituting Equation 29 into Equation 28 and minimizing with respect to {tilde over (θ)}k, yields the timing estimate used in step 44, in the form:





  
    
      
        
          
            θ
            ^
          
          k
        
        =
        
          
            
              arg
              ⁢
              
                
              
              ⁢
              max
            
            
              0
              ≤
              
                
                  θ
                  ~
                
                k
              
              ≤
              
                θ
                max
              
            
          
          ⁢
          
            {
            
              
                Υ
                k
              
              ⁡
              
                (
                
                  
                    θ
                    ~
                  
                  k
                
                )
              
            
            }
          
        
      
    
    
      
        (
        30
        )
      
    
  




where Υk({tilde over (θ)}k) is given by:





  
    
      
        
          
            Υ
            k
          
          ⁡
          
            (
            
              
                θ
                ~
              
              k
            
            )
          
        
        =
        
          
            
              
              
                
                  ∑
                  
                    v
                    =
                    0
                  
                  
                    V
                    -
                    1
                  
                
                ⁢
                
                  
                    
                      F
                      H
                    
                    ⁡
                    
                      (
                      v
                      )
                    
                  
                  ⁢
                  
                    
                      Φ
                      H
                    
                    (
                    
                      
                        
                          θ
                          ~
                        
                        k
                      
                      ,
                      v
                    
                    )
                  
                  ⁢
                  
                    
                      
                        S
                        ^
                      
                      k
                    
                    ⁡
                    
                      (
                      v
                      )
                    
                  
                
              
              
            
            2
          
          .
        
      
    
    
      
        (
        31
        )
      
    
  




The timing metric can be rewritten as:





  
    
      
        
          
            Υ
            k
          
          ⁡
          
            (
            
              
                θ
                ~
              
              k
            
            )
          
        
        =
        
          
            ∑
            
              l
              =
              
                
                  θ
                  ~
                
                k
              
            
            
              
                
                  θ
                  ~
                
                k
              
              +
              L
              -
              1
            
          
          ⁢
          
            
              
              
                
                  ∑
                  
                    v
                    =
                    0
                  
                  
                    V
                    -
                    1
                  
                
                ⁢
                
                  
                    
                      
                        s
                        ^
                      
                      k
                    
                    ⁡
                    
                      (
                      
                        v
                        ,
                        l
                      
                      )
                    
                  
                  ⁢
                  
                    e
                    
                      j
                      ⁢
                      
                        
                      
                      ⁢
                      
                        
                          2
                          ⁢
                          
                            
                          
                          ⁢
                          π
                          ⁢
                          
                            
                          
                          ⁢
                          l
                        
                        N
                      
                      ⁢
                      v
                    
                  
                
              
              
            
            2
          
        
      
    
    
      
        (
        32
        )
      
    
  




where ŝk(v,l) is the Q-point inverse DFT (IDFT) of the sequence {Ŝk(iq,v)}q=0Q−1, i.e.,





  
    
      
        
          
            
              s
              ^
            
            k
          
          ⁡
          
            (
            
              v
              ,
              l
            
            )
          
        
        =
        
          
            ∑
            
              q
              =
              0
            
            
              Q
              -
              1
            
          
          ⁢
          
            
              
                
                  S
                  ^
                
                k
              
              ⁡
              
                (
                
                  i
                  
                    q
                    ,
                    v
                  
                
                )
              
            
            ⁢
            
              
                e
                
                  j
                  ⁢
                  
                    
                      2
                      ⁢
                      
                        
                      
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      l
                    
                    Q
                  
                  ⁢
                  q
                
              
              .
            
          
        
      
    
    
      
        (
        33
        )
      
    
  




The algorithm of Equation 30 is referred to herein as the LS-based timing estimator (LS-TE). Note that for V=1 the timing metric reduces to:





  
    
      
        
          
            Υ
            k
          
          ⁡
          
            (
            
              
                θ
                ~
              
              k
            
            )
          
        
        =
        
          
            ∑
            
              l
              =
              
                
                  θ
                  ~
                
                k
              
            
            
              
                
                  θ
                  ~
                
                k
              
              +
              L
              -
              1
            
          
          ⁢
          
            
              
              
                
                  
                    s
                    ^
                  
                  k
                
                ⁡
                
                  (
                  
                    0
                    ,
                    l
                  
                  )
                
              
              
            
            2
          
        
      
    
    
      
        (
        34
        )
      
    
  




and becomes periodic in {tilde over (θ)}k with period Q. In such a case, the estimate {tilde over (θ)}k is affected by an ambiguity of multiples of Q. This ambiguity does not represent a serious problem as long as Q can be chosen greater than θmax. Unfortunately, in some applications this may not be the case. To avoid such an impairment, each tile must be composed by V>1 adjacent subcarriers.

FIG. 5 is a flowchart showing step 24 of FIG. 1 in greater detail. Step 24 aims at estimating the power levels of signals transmitted by the RSSs. In step 46, the ML estimation vector produced in step 42 and shown in Equation 24 is utilized, and power level estimation is accomplished on the basis of Equation 8 above. It is observed from Equation 24 that the entry Ŝk(iq,v) is an unbiased estimate of Sk(iq,v) with variance:

σk2(iq,v)=σ2[T−1({circumflex over (ω)})]μk(q,v),μk(q,v)   (35)

where μk(q,v)=(k−1)QV+qV+v+1. An unbiased estimate of |Sk(iq,v)|2 is thus given by |Ŝk(iq,v)|2−σk2(iq,v). This fact, together with Equations 8 and 35, leads to the following LS-based power estimation algorithm (referred to herein as “LS-PE”) utilized in step 48 to calculate the power levels:





  
    
      
        
          
            P
            ^
          
          k
        
        =
        
          
            1
            QV
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  0
                
                
                  V
                  -
                  1
                
              
              ⁢
              
                
              
              ⁢
              
                [
                
                  
                    
                      
                      
                        
                          
                            S
                            ^
                          
                          k
                        
                        ⁡
                        
                          (
                          
                            i
                            
                              q
                              ,
                              v
                            
                          
                          )
                        
                      
                      
                    
                    2
                  
                  -
                  
                    
                      
                        σ
                        ^
                      
                      k
                      2
                    
                    ⁡
                    
                      (
                      
                        i
                        
                          q
                          ,
                          v
                        
                      
                      )
                    
                  
                
                ]
              
            
          
        
      
    
    
      
        (
        36
        )
      
    
  




where {circumflex over (σ)}k2(iq,v)={circumflex over (σ)}2[T−1({circumflex over (ω)})]μk(q,v),μk(q,v) and {circumflex over (σ)}2 is given in Equation 17.

To assess the computational load (complexity) of LS-TE, it is useful to distinguish between a first stage, leading to the sequences {ŝk(v,l)}l=0Q−1 in Equation 33, for 0≦v≦V−1, and a second stage where these sequences are used to compute the timing metrics of the detected codes. First, it is observed that a complexity of order O(KN log2(N)) is required to determine the matrices Dm({circumflex over (ω)}k) for k=1, 2, . . . , K, while computing the matrix T({circumflex over (ω)}) and vector




  
    ∑
    
      m
      =
      0
    
    
      M
      -
      1
    
  
  ⁢
  
    
      
        D
        m
        H
      
      ⁡
      
        (
        
          ω
          ^
        
        )
      
    
    ⁢
    
      Y
      m
    
  




requires O(MK2Q3V3) and O(MKQ2V2) operations, respectively. Evaluation of Ŝ(θ) through Equation 22 is equivalent to solving a linear system of order KQV and needs O(K2Q2V2) operations. Finally, the complexity of the IDFT in Equation 33 is O(Q log2(Q)) for each value of v. Since most of the processing is devoted to the computation of T({circumflex over (ω)}), the overall complexity in the first stage can reasonably be approximated as O(MK2Q3V3). In the second stage of the estimation process, the timing metric of Equation 32 is computed for all candidate values {tilde over (θ)}k and all detected codes. Observing that O(VL) operations are required for each quantity Υk({tilde over (θ)}k), the computational burden of the second stage is O(KVLθmax). The overall complexity of LS-TE is summarized in the first row of Table 1 below. Reduced complexity timing delay estimates (“RC-TE”) are provided for comparison purposes. RC-TE is described below in connection with FIG. 6.












TABLE 1







Algorithm
First stage
Second stage









LS-TE
O(MK2Q3V3)
O(KVLθmax)



RC-TE
O(KQV(M + K))
O(KVL(Q + PV))



RC-TE (V = 2)
O(2KQ(M + K))
O(2LKQ)










FIG. 6 is a flowchart showing processing steps according to the present invention, indicated generally at 50, for estimating power levels and timing delays during initial ranging with reduced computational load (complexity). This approach relies on the assumption that the CFOs are small during the ranging process and result only in negligible phase rotations over one OFDM block. In such a case, an estimate of S(θ) can easily be obtained from the observations {Y(iq,v)} in Equation 11. In step 52, an ML estimation procedure is applied to {Y(iq,v)}. First, a K-dimensional vector S(iq,v)=[S1(iq,v),S2(iq,v), . . . ,SK(iq,v)]T is defined in order to rewrite Equation 11 in the equivalent form:

Y(iq,v)=C(ω)S(iq,v)+n(iq,v)   (37)

where C(ω)=[Γ(ω1)c1Γ(ω2)c2 . . . Γ(ωK)cK] is a matrix of dimensions M×K. Then, after letting {circumflex over (ω)}≈ω, Equation 37 is used to find the ML estimate in step 52 of S(iq,v) as:

Ŝ(iq,v)=[CH({circumflex over (ω)})C({circumflex over (ω)})]−1CH({circumflex over (ω)})Y(iq,v)   (38)

Substituting Equation 37 into Equation 38 yields:

Ŝ(iq,v)=S(iq,v)+ψ(iq,v)   (39)

where ψ(iq,v) is a zero-mean disturbance vector with covariance matrix σ2[CH({circumflex over (ω)})C({circumflex over (ω)})]−1. This means that, as long as the simplified model specified in Equation 37 is sufficiently accurate, Ŝ(iq,v) is an unbiased estimate of S(iq,v). Its entries {Ŝ(iq,v)}k=1K are then used to compute an estimate of the power level Pk as indicated in step 54:





  
    
      
        
          
            P
            ^
          
          k
        
        =
        
          
            1
            QV
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  0
                
                
                  V
                  -
                  1
                
              
              ⁢
              
                
              
              ⁢
              
                [
                
                  
                    
                      
                      
                        
                          
                            S
                            ^
                          
                          k
                        
                        ⁡
                        
                          (
                          
                            i
                            
                              q
                              ,
                              v
                            
                          
                          )
                        
                      
                      
                    
                    2
                  
                  -
                  
                    
                      
                        
                          σ
                          ^
                        
                        2
                      
                      ⁡
                      
                        [
                        
                          
                            
                              C
                              H
                            
                            ⁡
                            
                              (
                              
                                ω
                                ^
                              
                              )
                            
                          
                          ⁢
                          
                            C
                            ⁡
                            
                              (
                              
                                ω
                                ^
                              
                              )
                            
                          
                        
                        ]
                      
                    
                    
                      k
                      ,
                      k
                    
                    
                      -
                      1
                    
                  
                
                ]
              
            
          
        
      
    
    
      
        (
        40
        )
      
    
  




where {circumflex over (σ)}2 is still given in Equation 17.

To quantify the computational saving of the simplified ML estimation algorithm of Equation 38 utilized in step 52 with respect to Equation 22, it is observed that O(MK2) and O(MK) operations are needed to obtain CH({circumflex over (ω)})C({circumflex over (ω)}) and CH({circumflex over (ω)})Y(iq,v), respectively, while O(K2) operations are further required for evaluating the right hand side of Equation 38. Since CH({circumflex over (ω)})C({circumflex over (ω)}) is independent of q and v, the overall complexity involved in the computation of Ŝ(iq,v) for 0≦q≦Q−1 and 0≦v≦V−1 is approximately O(KQV(M+K)). Additional O(QV log2 Q) operations are eventually needed by the IDFTs in Equation 33. Compared to the first stage of LS-TE, it follows that the computational burden is reduced by a factor MKQ2V2/(M+K), which may be significant even for small values of the product QV.
Evaluating Υk({tilde over (θ)}k) for each value {tilde over (θ)}k=0,1, . . . ,θmax may still be computationally demanding, especially when θmax is large. A viable method to further reduce the system complexity is to decompose the timing error θk into a fractional part βk, less than Q, plus an integer part which is multiple of Q. This amounts to setting:

θk=βk+pkQ   (41)

where βk ε {0,1, . . . ,Q−1} while pk is an integer parameter taking values in the set {0,1, . . . ,P−1}, with P=└θmax/Q┘. In such a case, the timing metric of Equation 32 can be manipulated into the equivalent form:





  
    
      
        
          
            Υ
            k
          
          ⁡
          
            (
            
              
                θ
                ~
              
              k
            
            )
          
        
        =
        
          
            
              Υ
              k
              
                (
                1
                )
              
            
            ⁡
            
              (
              
                
                  β
                  ~
                
                k
              
              )
            
          
          +
          
            
              
                Υ
                k
                
                  (
                  2
                  )
                
              
              (
              
                
                  
                    β
                    ~
                  
                  k
                
                ,
                
                  
                    p
                    ~
                  
                  k
                
              
              )
            
            ⁢
            
              
            
            ⁢
            where
            ⁢
            
              :
            
          
        
      
    
    
      
        (
        42
        )
      
    
  
  
    
      
        
          
            Υ
            k
            
              (
              1
              )
            
          
          ⁡
          
            (
            
              
                θ
                ~
              
              k
            
            )
          
        
        =
        
          
            ∑
            
              l
              =
              
                
                  θ
                  ~
                
                k
              
            
            
              
                
                  θ
                  ~
                
                k
              
              +
              L
              -
              1
            
          
          ⁢
          
            
          
          ⁢
          
            
              ∑
              
                v
                =
                0
              
              
                V
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  
                  
                    
                      
                        s
                        ^
                      
                      k
                    
                    ⁡
                    
                      (
                      
                        v
                        ,
                        l
                      
                      )
                    
                  
                  
                
                2
              
              ⁢
              
                
              
              ⁢
              and
              ⁢
              
                :
              
            
          
        
      
    
    
      
        (
        43
        )
      
    
  
  
    
      
        
          
            Υ
            k
            
              (
              2
              )
            
          
          (
          
            
              
                β
                ~
              
              k
            
            ,
            
              
                p
                ~
              
              k
            
          
          )
        
        =
        
          2
          ⁢
          Re
          ⁢
          
            {
            
              
                ∑
                
                  l
                  =
                  
                    
                      β
                      ~
                    
                    k
                  
                
                
                  
                    
                      β
                      ~
                    
                    k
                  
                  +
                  L
                  -
                  1
                
              
              ⁢
              
                
              
              ⁢
              
                
                  ∑
                  
                    v
                    =
                    0
                  
                  
                    V
                    -
                    2
                  
                
                ⁢
                
                  
                
                ⁢
                
                  
                    ∑
                    
                      n
                      =
                      1
                    
                    
                      V
                      -
                      v
                      -
                      I
                    
                  
                  ⁢
                  
                    
                  
                  ⁢
                  
                    
                      
                        
                          s
                          ^
                        
                        k
                      
                      ⁡
                      
                        (
                        
                          v
                          ,
                          l
                        
                        )
                      
                    
                    ⁢
                    
                      
                        
                          s
                          ^
                        
                        k
                        *
                      
                      ⁡
                      
                        (
                        
                          
                            v
                            +
                            n
                          
                          ,
                          l
                        
                        )
                      
                    
                    ⁢
                    
                      e
                      
                        
                          -
                          j
                        
                        ⁢
                        
                          
                            2
                            ⁢
                            π
                            ⁢
                            
                              
                            
                            ⁢
                            n
                          
                          N
                        
                        ⁢
                        
                          (
                          
                            l
                            +
                            
                              
                                
                                  p
                                  ~
                                
                                k
                              
                              ⁢
                              Q
                            
                          
                          )
                        
                      
                    
                  
                
              
            
            }
          
        
      
    
    
      
        (
        44
        )
      
    
  




Clearly, maximizing Υk({tilde over (θ)}k) with respect to {tilde over (θ)}k requires a grid search over all candidate pairs ({tilde over (β)}k,{tilde over (p)}k). However, the decomposition of Equation 42 provides a reduced-complexity solution to approach the global maximum. This simplified algorithm computes an estimate in step 56 of timing delays βk by looking for the maximum of Υk(l)({tilde over (β)}k), i.e.,





  
    
      
        
          
            β
            ^
          
          k
        
        =
        
          
            
              arg
              ⁢
              
                
              
              ⁢
              max
            
            
              0
              ≤
              
                
                  β
                  ~
                
                k
              
              ≤
              
                Q
                -
                1
              
            
          
          ⁢
          
            {
            
              
                Υ
                k
                
                  (
                  1
                  )
                
              
              ⁡
              
                (
                
                  
                    β
                    ~
                  
                  k
                
                )
              
            
            }
          
        
      
    
    
      
        (
        45
        )
      
    
  




Then, after substituting {circumflex over (β)}k into the right-hand-side of Equation 42 and maximizing with respect to {tilde over (p)}k, an estimate of pk is obtained in step 56 in the form:





  
    
      
        
          
            p
            ^
          
          k
        
        =
        
          
            
              arg
              ⁢
              
                
              
              ⁢
              max
            
            
              0
              ≤
              
                
                  p
                  ~
                
                k
              
              ≤
              
                P
                -
                1
              
            
          
          ⁢
          
            {
            
              
                Υ
                k
                
                  (
                  2
                  )
                
              
              (
              
                
                  
                    β
                    ^
                  
                  k
                
                ,
                
                  
                    p
                    ~
                  
                  k
                
              
              )
            
            }
          
        
      
    
    
      
        (
        46
        )
      
    
  




The timing delay estimate is eventually obtained from Equation 41 as {circumflex over (θ)}k={circumflex over (β)}k+{circumflex over (p)}kQ.

Inspection of Equations 45 and 46 indicates that the simplified estimator maximizes Υk({tilde over (θ)}k) in a decoupled fashion, where an estimate of βk is obtained first and is then exploited to get {circumflex over (p)}k. In assessing the computational load of this scheme, it is assumed that the quantities {ŝk(v,l)}l=0Q−1 have already been computed for 0≦v≦V−1. Then, evaluating {circumflex over (β)}k through Equation 45 requires O(QVL) operations while O(QV2L) operations are needed to get {circumflex over (p)}k from Equation 46. Since estimates {circumflex over (β)}k and {circumflex over (p)}k must be obtained for each detected code, the overall complexity in the second stage of the estimation process is approximately O(KVL(Q+PV)). Compared to LS-TE, the computational requirement is thus reduced by a factor θmax/(Q+PV). This advantage is achieved at the price of some performance degradation since the described search procedure is not guaranteed to locate the global maximum of Υk({tilde over (θ)}k).
A further simplification is possible by setting V=2. In this case, it is found from Equation 44 that maximizing Υk(2)({circumflex over (β)}k,{tilde over (p)}k) with respect to {tilde over (p)}k is equivalent to maximizing cos(φk−2π{tilde over (p)}kQ/N), with.




  
    
      
        
          φ
          k
        
        =
        
          arg
          ⁢
          
            {
            
              
                ∑
                
                  l
                  =
                  
                    
                      β
                      ~
                    
                    k
                  
                
                
                  
                    
                      β
                      ~
                    
                    k
                  
                  +
                  L
                  -
                  1
                
              
              ⁢
              
                
              
              ⁢
              
                
                  
                    
                      s
                      ^
                    
                    k
                  
                  ⁡
                  
                    (
                    
                      0
                      ,
                      l
                    
                    )
                  
                
                ⁢
                
                  
                    
                      s
                      ^
                    
                    k
                    *
                  
                  ⁡
                  
                    (
                    
                      1
                      ,
                      l
                    
                    )
                  
                
                ⁢
                
                  e
                  
                    
                      -
                      j
                    
                    ⁢
                    
                      
                        2
                        ⁢
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        l
                      
                      N
                    
                  
                
              
            
            }
          
        
      
    
    
      
        (
        47
        )
      
    
  




The estimate {circumflex over (p)}k is thus obtained in closed-form as:





  
    
      
        
          
            p
            ^
          
          k
        
        =
        
          
            
              N
              ⁢
              
                
              
              ⁢
              
                φ
                k
              
            
            
              2
              ⁢
              π
              ⁢
              
                
              
              ⁢
              Q
            
          
          .
        
      
    
    
      
        (
        48
        )
      
    
  




The reduced-complexity power estimation provided by the algorithm of Equation 40, and implemented in step 54, is referred to herein as the reduced-complexity power estimator (RC-PE).

FIG. 7 is a block diagram showing a sample hardware and software/firmware arrangement for implementing the present invention, indicated generally at 90. The present invention can be implemented as a plurality of software and/or firmware modules 96 which are executed by a standard OFDMA base station transceiver 94 having an antenna system 92. The base station receiver 94 could be any suitable base station transceiver which operates using OFDM technology. The modules 96 could be coded in any suitable high-level or low-level programming language, such as C, C++, C#, Java, assembly language, or any suitable equivalent. The modules could be provided in an OFDMA transceiver in the form of a customized integrated circuit (IC), or distributed using any suitable computer-readable medium (e.g., diskette, tape, CD-ROM, DVD, portable flash memory device, etc.) for subsequent installation of the modules in an existing OFDMA transceiver. Still farther, the modules could be downloaded to an OFDMA transceiver by way of a network connection (e.g., via the Internet).
The modules 96 are programmed in accordance with the processing steps discussed above, and include an active code identification module 98, a CFO estimation module 100, a timing delay estimation module 102, and a power level estimation module 104. The active code identification module 98 processes the ranging OFDM blocks received by the transceiver 94 from RSSs in communication with the transceiver 94, to identify active codes. The active codes are then processed by the modules 100-104 to provide estimates of the corresponding CFOs, timing delays, and power levels. These estimates are then broadcast by the transceiver 94 to the RSSs so that each one can adjust its transmission parameters to the BS reference scales 94.
Software simulation testing was conducted to ascertain the performance characteristics of the present invention. Simulations of an IEEE 802.16 wireless MAN programmed in accordance with the present invention was conducted, which had a total of N=1024 subcarriers over an uplink bandwidth of 11.4 MHz. The sampling period is Ts=87.5 ns, corresponding to a subcarrier distance of 1/(NTs)=11.16 kHz. It was assumed that R=8 subchannels were available for initial ranging. Each subchannel was divided into Q=16 tiles uniformly spaced over the signal spectrum at a distance of N/Q=64 subcarriers. To reduce the system complexity, the number of subcarriers in any tile was set to V=2. Each ranging time-slot included M=4 OFDM blocks. The ranging codes were taken from a Fourier set of length 4 and are randomly selected by the RSSs at each simulation run (expect for the sequence [1,1,1,1]T). The discrete-time CIRs had order L=12. Their entries were modeled as independent and circularly symmetric Gaussian random variables with zero means and an exponential power delay profile, i.e.,

E{|hk(l)|2}=σh,k2exp(−l/12), l=0,1, . . . ,11   (49)

where σh,k2 is related to the average transmission power of the kth RSS. In all subsequent simulations, σh,k2 was chosen such that E{∥hk∥2}=1 for 1≦k≦K. In this way, the powers {Pk}k=1K of the ranging signals are random variables with unit means. Channels of different users were assumed to be statistically independent of each other. They were generated at each new simulation run and kept fixed over an entire time-slot. A cell radius of 10 km was considered, corresponding to a maximum propagation delay of θmax=204 sampling periods. Ranging blocks were preceded by a CP of length NG=256 in order to avoid IBI. The normalized CFOs were uniformly distributed over the interval [−Ω,Ω] and varied at each run. Recalling that the estimation range of MFE is |εk|<N/(2MNT), Ω≦0.1 was set throughout the simulations. The results of simulations are shown in FIGS. 8-14.

Simulations began by investigating the performance of MCD in terms of probability of making an incorrect detection, say Pf. The latter is defined as the probability that the set of detected codes does not coincide with the set of codes participating in the ranging process. FIG. 8 illustrates Pf as a function of the signal-to-noise ratio (SNR) defined as SNR=1/σ2. The number of active RSSs is either 2 or 3 while the maximum CFO is Ω=0.05. Comparisons were made with the ranging scheme proposed by Fu, Li and Minn (FLM), where the kth ranging code is declared active provided that the quantity:




  
    
      
        
          Z
          k
        
        =
        
          
            1
            
              M
              2
            
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  0
                
                
                  V
                  -
                  1
                
              
              ⁢
              
                
              
              ⁢
              
                
                  
                  
                    
                      c
                      k
                      H
                    
                    ⁢
                    
                      Y
                      ⁡
                      
                        (
                        
                          i
                          
                            q
                            ,
                            v
                          
                        
                        )
                      
                    
                  
                  
                
                2
              
            
          
        
      
    
    
      
        (
        50
        )
      
    
  




exceeds a suitable threshold η which is proportional to the estimated noise power {circumflex over (σ)}2. The results of FIG. 8 indicate that MCD performs remarkably better than FLM. As expected, the system performance deteriorates as K approaches M. The reason is that increasing K reduces the dimensionality of the noise subspace, which degrades the accuracy of MCD.

FIG. 9 illustrates the root mean-square error (RMSE) of the frequency estimates obtained with MFE versus SNR for K=2 or 3 and Ω=0.05. It can be seen that the system performance deteriorates as K increases. Nevertheless, the accuracy of MFE is satisfactory under all investigated conditions.
The impact of the CFOs on the performance of MFE was assessed in FIG. 10, where the RMSE is shown as a function of Ω. The SNR was fixed to 10 dB while K was still 2 or 3. As expected, the accuracy of the frequency estimates degrades with Ω due to increased amount of inter-channel interference (ICI) that affects the received signal. However, the loss is quite small over the full range Ω≦0.1.
The performance of the timing estimators was measured in terms of probability of making a timing error, say, P(ε). An error event was declared to occur whenever the estimate {circumflex over (θ)}k gave rise to IBI during the data section of the frame. In such a case, the quantity {circumflex over (θ)}k−θk+(L−NG,D)/2 is larger than zero or smaller than L−NG,D−1, where NG,D is the CP length during the data transmission period. In our simulations, we set NG,D=64. FIG. 11 illustrates P(ε) vs. SNR as obtained with RC-TE, LS-TE and FLM. The number of active codes was K=2 or 3 while Ω=0.05. It can be seen that for SNR larger than 6 dB RC-TE and LS-TE provided much better results than FLM. Also, the loss of RC-TE with respect to LS-TE is negligible in all considered situations.
FIG. 12 shows the timing estimators compared in terms of their sensitivity to CFOs. For this purpose, P(ε) is shown as a function of Ω for K=2 or 3 and SNR=10 dB. Again, it can be seen that LS-TE was slightly better than RC-TE and both schemes outperformed ELM. For K=2 and 3, RC-TE and LS-TE were quite robust to CFOs.
FIG. 13 illustrates the RMSE of the power estimates as a function of SNR. The number of active RSSs was either 2 or 3 while Ω=0.05. Comparisons were made with the known FLM algorithm, where an estimate of Pk is computed as {circumflex over (P)}k=Zk/(QV)−{circumflex over (σ)}2/M. As can be seen, in all considered scenarios, RC-PE and LS-PE outperformed FLM. The accuracy of RC-PE and LS-PE reduced as K grew large, while the performance of FLM was nearly the same with either K=2 or 3.
The impact of residual CFOs on the RMSE of the power estimates was assessed as shown in FIG. 14, where the RMSE is shown as a function of Ω. Again, the SNR was fixed to 10 dB and K was either 2 or 3. As can be seen, LS-PE performed slightly better than RC-PE. Both schemes were extremely robust to CFOs when K is either 2 or 3.
FIG. 15 is a flowchart showing another embodiment of the present invention, indicated generally at 120, wherein CFOs, timing errors, and power levels are estimated using the known EPSRIT algorithm. The ESPRIT algorithm is beneficial because it can be implemented with reduced computational complexity, and is suitable for use in systems where pilot symbols are spread in both the time and frequency domains using ranging codes belonging to a Fourier basis.
The embodiment shown in FIG. 15 can be employed in an OFDMA system employing N subcarriers with index set {0, 1, . . . , N−1}. It is assumed that a ranging time-slot is composed by M consecutive OFDMA blocks while the N available subcarriers are grouped into ranging subchannels and data subchannels. The former are used by the active RSSs to initiate their ranging processes, while the latter are assigned to DSSs for data transmission. R denotes the number of ranging subchannels, and it is assumed that each of them is divided into Q subbands. A given subband is composed of a set of V adjacent subcarriers, which is called a tile. The subcarrier indices of the qth tile (q=0, 1, . . . , Q−1) in the rth subchannel (r=0, 1, . . . , R−1) are collected into a set Jq(r)={iq(r)+v}v=0V−1, where the tile index iq(r) can be chosen adaptively according to the actual channel conditions. The only constraint in the selection of iq(r) is that different tiles must be disjoint, i.e., Jq1(r1)∩Jq2(r2)=ø for q1≠q2 or r1≠r2. The rth ranging subchannel is thus composed of QV subcarriers with indices in the set J(r)=∪q=0Q−1Jq(r), while a total of NR=QVR ranging subcarriers is available in each OFDMA block.
It is assumed that each subchannel can be accessed by a maximum number of Kmax=min{V,M}−1. RSSs, which are separated by means of orthogonal codes in both the time and frequency domains. The codes are selected in a pseudo-random fashion from a predefined set {C0,C1, . . . ,CKmax−1} with




  
    
      
        
          
            [
            
              C
              k
            
            ]
          
          
            v
            ,
            m
          
        
        =
        
          e
          
            j
            ⁢
            
              
            
            ⁢
            2
            ⁢
            π
            ⁢
            
              
            
            ⁢
            
              k
              ⁡
              
                (
                
                  
                    v
                    
                      V
                      -
                      1
                    
                  
                  +
                  
                    m
                    
                      M
                      -
                      1
                    
                  
                
                )
              
            
          
        
      
    
    
      
        (
        51
        )
      
    
  




where v=0, 1, . . . , V−1 counts the subcarriers within a tile and is used to perform spreading in the frequency domain, while m=0, 1, . . . , M−1 is the block index by which spreading is done in the time domain across the ranging time-slot. It is assumed that different RSSs select different codes. Also, it is assumed that a selected code is employed by the corresponding RSS over all tiles in the considered subchannel. Without loss of generality, we concentrate on the rth subchannel and denote by K≦Kmax the number of simultaneously active RSSs. To simplify the notation, the subchannel index (r) is dropped henceforth.

The signal transmitted by the kth RSS propagates through a multipath channel characterized by a channel impulse response (CIR) hk=[hk(0),hk(1), . . . ,hk(L−1)]T of length L (in sampling periods). We denote by θk the timing error of the kth RSS expressed in sampling intervals Ts, while εk is the frequency offset normalized to the subcarrier spacing. During IR, the CFOs are only due to Doppler shifts and/or to estimation errors and, in consequence, they are assumed to lie within a small fraction of the subcarrier spacing. Timing offsets depend on the distance of the RSSs from the BS and their maximum value is thus limited to the round trip delay from the cell boundary. In order to eliminate IBI, it is assumed that, during the ranging process, the CP length comprises NG≧θmax+L sampling periods, where θmax is the maximum expected timing error. This assumption is not restrictive as initialization blocks are usually preceded by long CPs in many standardized OFDM systems.
In step 126, a ranging subchannel is extracted from the DFT outputs corresponding to the received OFDM blocks. We denote by Xm(q)=[Xm(iq),Xm(iq+1), . . . ,Xm(iq+V−1)]T the DFT outputs corresponding to the qth tile in the mth OFDMA block. Since DSSs have successfully completed their IR processes, they are perfectly aligned to the BS references and their signals do not contribute to Xm(q). In contrast, the presence of uncompensated CFOs destroys orthogonality among ranging signals and gives rise to ICI. The latter results in a disturbance term plus an attenuation of the useful signal component. To simplify the analysis, in the ensuing discussion the disturbance term is treated as a zero-mean Gaussian random variable while the signal attenuation is considered as part of the channel impulse response. Under the above assumptions, we may write




  
    
      
        
          
            X
            m
          
          ⁡
          
            (
            
              
                i
                q
              
              +
              v
            
            )
          
        
        =
        
          
            
              ∑
              
                k
                =
                1
              
              K
            
            ⁢
            
              
            
            ⁢
            
              
                
                  [
                  
                    C
                    
                      l
                      k
                    
                  
                  ]
                
                
                  v
                  ,
                  m
                
              
              ⁢
              
                e
                
                  j
                  ⁢
                  
                    
                  
                  ⁢
                  m
                  ⁢
                  
                    
                  
                  ⁢
                  
                    ω
                    k
                  
                  ⁢
                  
                    N
                    T
                  
                
              
              ⁢
              
                
                  d
                  k
                
                ⁡
                
                  (
                  q
                  )
                
              
              ⁢
              
                
                  H
                  k
                
                ⁡
                
                  (
                  
                    
                      θ
                      k
                    
                    ,
                    
                      ɛ
                      k
                    
                    ,
                    
                      
                        i
                        q
                      
                      +
                      v
                    
                  
                  )
                
              
            
          
          +
          
            
              w
              m
            
            ⁡
            
              (
              
                
                  i
                  q
                
                +
                v
              
              )
            
          
        
      
    
    
      
        (
        52
        )
      
    
  




where ωk−2πεk/N, NT=N+NG denotes the duration of the cyclically extended block and Clk is the code matrix selected by the kth RSS. The quantity dk(q) is a unit-amplitude PSK symbol transmitted by the kth RSS on the qth tile while Hk(θk,εk,n) is the kth equivalent channel frequency response over the nth subcarrier and is given by

Hk(θk,εk,n)=γN(εk)Hk′(n)e−j2πnθk/N   (53)

where





  
    
      H
      k
      ′
    
    ⁡
    
      (
      n
      )
    
  
  =
  
    
      ∑
      
        l
        =
        0
      
      
        L
        -
        1
      
    
    ⁢
    
      
        
          h
          k
          ′
        
        ⁡
        
          (
          l
          )
        
      
      ⁢
      
        e
        
          
            -
            j2π
          
          ⁢
          
            
          
          ⁢
          
            nl
            /
            N
          
        
      
    
  




is the true channel frequency response, while





  
    
      
        
          
            γ
            N
          
          ⁡
          
            (
            ɛ
            )
          
        
        =
        
          
            
              sin
              ⁡
              
                (
                
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  ɛ
                
                )
              
            
            
              N
              ⁢
              
                
              
              ⁢
              
                sin
                ⁡
                
                  (
                  
                    π
                    ⁢
                    
                      
                    
                    ⁢
                    
                      ɛ
                      /
                      N
                    
                  
                  )
                
              
            
          
          ⁢
          
            e
            
              j
              ⁢
              
                
              
              ⁢
              π
              ⁢
              
                
              
              ⁢
              
                
                  ɛ
                  ⁡
                  
                    (
                    
                      N
                      -
                      1
                    
                    )
                  
                
                /
                N
              
            
          
        
      
    
    
      
        (
        54
        )
      
    
  




is the attenuation factor induced by the CFO. In many practical situations, the contribution of γN(ε) can be ignored since its magnitude is close to unity. As an example, for N=1024 and |ε|<0.1 we have |γN(ε)|>0.98. The last term in Equation 52 accounts for background noise and is modeled as a circularly symmetric complex Gaussian random variable with zero-mean and variance σw2=σn2+σICI2, where σn2 and σICI2 are the average noise and ICI powers, respectively. As mentioned earlier, the code matrix selected by each RSS is employed over all tiles in the considered subchannel. Furthermore, from Equation 53 we see that θk appears only as a phase shift across the DFT outputs. The reason is that the CP duration is longer than the maximum expected propagation delay, thereby making the ranging signals quasi-synchronous at the BS.

To proceed further, it is assumed that the tile width is much smaller than the channel coherence bandwidth. In this case, the channel response is nearly flat over each tile and we may reasonably replace the quantities {Hk′(iq+v)}v=0V−1 with an average frequency response




  
    
      
        
          
            H
            _
          
          q
          ′
        
        =
        
          
            1
            V
          
          ⁢
          
            
              ∑
              
                v
                =
                0
              
              
                V
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  H
                  k
                  ′
                
                ⁡
                
                  (
                  
                    
                      i
                      q
                    
                    +
                    v
                  
                  )
                
              
              .
            
          
        
      
    
    
      
        (
        55
        )
      
    
  




Substituting Equations 51 and 53 into Equation 52, and bearing in mind Equation 55, yields





  
    
      
        
          
            X
            m
          
          ⁡
          
            (
            
              
                i
                q
              
              +
              v
            
            )
          
        
        =
        
          
            
              ∑
              
                k
                =
                1
              
              K
            
            ⁢
            
              
            
            ⁢
            
              
                e
                
                  j
                  ⁢
                  
                    
                  
                  ⁢
                  2
                  ⁢
                  
                    π
                    ⁡
                    
                      (
                      
                        
                          m
                          ⁢
                          
                            
                          
                          ⁢
                          
                            ξ
                            k
                          
                        
                        +
                        
                          v
                          ⁢
                          
                            
                          
                          ⁢
                          
                            η
                            k
                          
                        
                      
                      )
                    
                  
                
              
              ⁢
              
                
                  S
                  k
                
                ⁡
                
                  (
                  q
                  )
                
              
            
          
          +
          
            
              w
              m
            
            ⁡
            
              (
              
                
                  i
                  q
                
                +
                v
              
              )
            
          
        
      
    
    
      
        (
        56
        )
      
    
  




where Sk(q)=γN(εk) Hk′(q)e−j2πiqθk/N and we have defined the quantities





  
    
      
        
          
            ξ
            k
          
          =
          
            
              
                1
                k
              
              
                M
                -
                1
              
            
            +
            
              
                
                  ɛ
                  k
                
                ⁢
                
                  N
                  T
                
              
              N
            
          
        
        ⁢
        
          

        
        ⁢
        and
      
    
    
      
        (
        57
        )
      
    
  
  
    
      
        
          η
          k
        
        =
        
          
            
              1
              k
            
            
              V
              -
              1
            
          
          -
          
            
              θ
              k
            
            N
          
        
      
    
    
      
        (
        58
        )
      
    
  




which are referred to as the effective CFOs and timing errors, respectively.

As discussed below, the DFT outputs {Xm(iq+v)} can be exploited to identify the active codes and to estimate the corresponding timing errors and CFOs.
In step 128, the number K of active codes over the considered ranging subchannel is determined. For this purpose, we collect the (iq+v) th DFT outputs across all ranging blocks into an M-dimensional vector Y(iq,v)=[X0(iq+v),X1(iq+v), . . . ,XM−1(iq+v)|T given by




  
    
      
        
          Y
          ⁡
          
            (
            
              i
              
                q
                ,
                v
              
            
            )
          
        
        =
        
          
            
              ∑
              
                k
                =
                1
              
              K
            
            ⁢
            
              
            
            ⁢
            
              
                e
                
                  j
                  ⁢
                  
                    
                  
                  ⁢
                  2
                  ⁢
                  π
                  ⁢
                  
                    
                  
                  ⁢
                  v
                  ⁢
                  
                    
                  
                  ⁢
                  
                    η
                    k
                  
                
              
              ⁢
              
                
                  S
                  k
                
                ⁡
                
                  (
                  q
                  )
                
              
              ⁢
              
                
                  e
                  M
                
                ⁡
                
                  (
                  
                    ξ
                    k
                  
                  )
                
              
            
          
          +
          
            w
            ⁡
            
              (
              
                i
                
                  q
                  ,
                  v
                
              
              )
            
          
        
      
    
    
      
        (
        59
        )
      
    
  




where w(iq,v)=[w0(iq+v),w1(iq+v), . . . ,wM−1(iq+v)]T is Gaussian distributed with zero mean and covariance matrix σw2IM while eM(ξ)=[1,ej2πξ,ej4πξ, . . . ,ej2π(M−1)ξ]T.

From Equation 59 above, it can be observed that Y(iq,v) has the same structure as measurements of multiple uncorrelated sources from an array of sensors. Hence, an estimate of K can be obtained by performing an EVD of the correlation matrix RY=E{Y(iq,v)YH(iq,v)}. In practice, however, RY is not available at the receiver and must be replaced by some suitable estimate. One popular strategy to get an estimate of RY is based on the FB principle. Following this approach, RY is replaced by




  
    
      
        
          
            R
            ^
          
          Y
        
        =
        
          
            1
            2
          
          ⁢
          
            (
            
              
                
                  R
                  ~
                
                Y
              
              +
              
                J
                ⁢
                
                  
                    R
                    ~
                  
                  Y
                
                ⁢
                
                  J
                  T
                
              
            
            )
          
        
      
    
    
      
        (
        60
        )
      
    
  




where {tilde over (R)}Y is the sample correlation matrix





  
    
      
        
          
            R
            ~
          
          Y
        
        =
        
          
            1
            QV
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  v
                  =
                  0
                
                
                  V
                  -
                  1
                
              
              ⁢
              
                
              
              ⁢
              
                
                  Y
                  ⁡
                  
                    (
                    
                      i
                      
                        q
                        ,
                        v
                      
                    
                    )
                  
                
                ⁢
                
                  
                    Y
                    H
                  
                  ⁡
                  
                    (
                    
                      i
                      
                        q
                        ,
                        v
                      
                    
                    )
                  
                
              
            
          
        
      
    
    
      
        (
        61
        )
      
    
  




while J is the exchange matrix with 1's on the anti-diagonal and 0's elsewhere. Arranging the eigenvalues {circumflex over (λ)}1≧{circumflex over (λ)}2≧ . . . ≧{circumflex over (λ)}M−1 of {circumflex over (R)}Y in non-increasing order, an estimate {circumflex over (K)} of the number of active codes can be found by applying the MDL information-theoretic criterion. This amounts to looking for the minimum of the following objective function





  
    
      
        
          F
          ⁡
          
            (
            
              K
              ~
            
            )
          
        
        =
        
          
            
              1
              2
            
            ⁢
            
              
                K
                ~
              
              (
              
                
                  2
                  ⁢
                  M
                
                -
                
                  K
                  ~
                
              
              )
            
            ⁢
            
              ln
              ⁡
              
                (
                QV
                )
              
            
          
          -
          
            
              QV
              (
              
                M
                -
                
                  K
                  ~
                
              
              )
            
            ⁢
            ln
            ⁢
            
              
            
            ⁢
            
              ρ
              ⁡
              
                (
                
                  K
                  ~
                
                )
              
            
          
        
      
    
    
      
        (
        62
        )
      
    
  




where ρ({tilde over (K)}) is the ratio between the geometric and arithmetic means of {{circumflex over (λ)}{tilde over (K)}−1,{circumflex over (λ)}{tilde over (K)}+2, . . . ,{circumflex over (λ)}M}.

In step 130, CFOs are estimated using the ESPRIT algorithm. For simplicity, it is assumed that the number of active codes has been perfectly estimated. Then, an estimate of ξ=[ξ1,ξ2, . . . ,ξK]T can be found by applying the ESPRIT algorithm to the model set forth in Equation 59. The eigenvectors of {circumflex over (R)}Y associated to the K largest eigenvalues {circumflex over (λ)}1≧{circumflex over (λ)}2≧ . . . ≧{circumflex over (λ)}K are arranged into an M×K matrix Z=[z1z2 . . . zK]. Next, we consider the matrices Z1 and Z2 that are obtained by collecting the first M−1 rows and the last M−1 rows of Z, respectively. The entries of ξ are finally estimated in a decoupled fashion as




  
    
      
        
          
            
              ξ
              ^
            
            k
          
          =
          
            
              1
              
                2
                ⁢
                π
              
            
            ⁢
            arg
            ⁢
            
              {
              
                
                  ρ
                  y
                
                ⁡
                
                  (
                  k
                  )
                
              
              }
            
          
        
        ,
        
          

        
        ⁢
        
          k
          =
          1
        
        ,
        2
        ,
        ...
        ⁢
        
          
        
        ,
        K
      
    
    
      
        (
        63
        )
      
    
  




where {ρy(1),ρy(2), . . . ,ρy(K)} are the eigenvalues of

ZY=(Z1HZ1H)−1Z1HZ2   (64)

and arg{ρy(k)} denotes the phase angle of ρy(k) taking values in the interval [−π,π).

After computing estimates of the effective CFOs through Equation 63, the problem arises of matching each {circumflex over (ξ)}k to the corresponding code Clk. This amounts to finding an estimate of lk starting from {circumflex over (ξ)}k. For this purpose, we denote by |εmax| the magnitude of the maximum expected CFO and observe from Equation 57 that (M−1)ξk belongs to the interval Jlk=[lk−β;1k+β], with β=|εmax|NT(M−1)/N. It follows that the effective CFOs can be univocally mapped to their corresponding codes as long as β<1/2 since only in that case intervals {Jlk}k=1K are disjoint. The acquisition range of the frequency estimator is thus limited to |εmax|<N/(2NT(M−1)) and an estimate of the pair (lk,εk) is computed as




  
    
      
        
          
            l
            ^
          
          k
        
        =
        
          
            round
            (
            
              
                (
                
                  M
                  -
                  1
                
                )
              
              ⁢
              
                
                  ξ
                  ^
                
                k
              
            
            )
          
          ⁢
          
            
          
          ⁢
          and
        
      
    
    
      
        (
        65
        )
      
    
  
  
    
      
        
          
            ɛ
            ^
          
          k
        
        =
        
          
            N
            
              N
              T
            
          
          ⁢
          
            (
            
              
                
                  ξ
                  ^
                
                k
              
              -
              
                
                  
                    l
                    ^
                  
                  k
                
                
                  M
                  -
                  1
                
              
            
            )
          
        
      
    
    
      
        (
        66
        )
      
    
  




It is worth noting that the arg{•} function in Equation 63 has an inherent ambiguity of multiples of 2π, which translates into a corresponding ambiguity of the quantity {circumflex over (l)}k by multiples of M−1. Hence, recalling that lk ε{0,1, . . . ,Kmax−1} with Kmax<M, a refined estimate of lk can be found as

{circumflex over (l)}k(F)=[{circumflex over (l)}k]M−1   (67)

where [x]M−1 is the value of x reduced to the interval [0, M−2]. Equation 66 is referred to below as the ESPRIT-based frequency estimator (EFE).

In step 132, timing errors are estimated next. We call Xm(q)=[Xm(iq),Xm(iq+1), . . . ,Xm(iq+V−1)]T the V-dimensional vector of the DFT outputs corresponding to the qth tile in the mth OFDM block. Then, from Equation 56, we have




  
    
      
        
          
            X
            m
          
          ⁡
          
            (
            q
            )
          
        
        =
        
          
            
              ∑
              
                k
                =
                1
              
              K
            
            ⁢
            
              
            
            ⁢
            
              
                e
                
                  
                    j
                    ⁢
                    
                      
                    
                    ⁢
                    2
                    ⁢
                    π
                    ⁢
                    
                      
                    
                    ⁢
                    m
                    ⁢
                    
                      
                    
                    ⁢
                    
                      ξ
                      k
                    
                  
                  ⁢
                  
                    
                  
                
              
              ⁢
              
                
                  S
                  k
                
                ⁡
                
                  (
                  q
                  )
                
              
              ⁢
              
                
                  e
                  V
                
                ⁡
                
                  (
                  
                    η
                    k
                  
                  )
                
              
            
          
          +
          
            
              w
              m
            
            ⁡
            
              (
              q
              )
            
          
        
      
    
    
      
        (
        68
        )
      
    
  




where wm(q)=[wm(iq),wm(iq+1), . . . ,wm(iq+V−1)]T is Gaussian distributed with zero mean and covariance matrix σw2Iv while eV(η)=[l,ej2πη,ej4πη, . . . ,eJ2π(V−1)η]T. Since Xm(q) is a superposition of complex sinusoidal signals with random amplitudes embedded in white Gaussian noise, an estimate of η=[η1,η2, . . . ,ηK]T can still be obtained by resorting to the ESPRIT algorithm. Following the previous steps, we first compute





  
    
      
        R
        ^
      
      X
    
    =
    
      
        1
        2
      
      ⁢
      
        (
        
          
            
              R
              ~
            
            X
          
          +
          
            J
            ⁢
            
              
                R
                ~
              
              X
              T
            
            ⁢
            J
          
        
        )
      
    
  
  ,




with





  
    
      
        
          
            R
            ~
          
          X
        
        =
        
          
            1
            MQ
          
          ⁢
          
            
              ∑
              
                m
                =
                0
              
              
                M
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                ∑
                
                  q
                  =
                  0
                
                
                  Q
                  -
                  1
                
              
              ⁢
              
                
              
              ⁢
              
                
                  
                    X
                    m
                  
                  ⁡
                  
                    (
                    q
                    )
                  
                
                ⁢
                
                  
                    
                      X
                      m
                      H
                    
                    ⁡
                    
                      (
                      q
                      )
                    
                  
                  .
                
              
            
          
        
      
    
    
      
        (
        69
        )
      
    
  




Then, we define a V×K matrix U=[u1u2 . . . uK] whose columns are the eigenvectors of {circumflex over (R)}X associated to the K largest eigenvalues. The effective timing errors are eventually estimated as





  
    
      
        
          
            
              η
              ^
            
            k
          
          =
          
            
              1
              
                2
                ⁢
                π
              
            
            ⁢
            arg
            ⁢
            
              {
              
                
                  ρ
                  x
                
                ⁡
                
                  (
                  k
                  )
                
              
              }
            
          
        
        ,
        
          

        
        ⁢
        
          k
          =
          1
        
        ,
        2
        ,
        ...
        ⁢
        
          
        
        ,
        K
      
    
    
      
        (
        70
        )
      
    
  




where {ρX(1),ρX(2), . . . ,ρX(K)} are the eigenvalues of

UX=(U1HU1)−1U1HU2   (71)

while the matrices U1 and U2 are obtained by collecting the first V−1 rows and the last V−1 rows of U, respectively.

The quantities {{circumflex over (η)}k}k−1K are eventually used to find estimates ({circumflex over (l)}k,{circumflex over (θ)}k) of the associated ranging code and timing error. To accomplish this task, we let α−θmax(V−1)/(2N). Then, recalling that 0≦θk≦θmax, from Equation 58, it can be seen that (V−1)ηk+α falls into the range Ilk=[lk−α;lk+α]. If θmax<N/(V−1), the quantity α is smaller than 1⁄2 and, in consequence, intervals {Ilk}k=1K are disjoint. In this case, there is only one pair (lk,θk) that results into a given ηk and an estimate of (lk,θk) is found as




  
    
      
        
          
            
              l
              ^
            
            k
          
          =
          
            round
            ⁡
            
              (
              
                
                  
                    (
                    
                      V
                      -
                      1
                    
                    )
                  
                  ⁢
                  
                    
                      η
                      ^
                    
                    k
                  
                
                +
                α
              
              )
            
          
        
        ⁢
        
          

        
        ⁢
        and
      
    
    
      
        (
        72
        )
      
    
  
  
    
      
        
          
            θ
            ^
          
          k
        
        =
        
          N
          (
          
            
              
                
                  l
                  ^
                
                k
              
              
                V
                -
                1
              
            
            -
            
              
                η
                ^
              
              k
            
          
          )
        
      
    
    
      
        (
        73
        )
      
    
  




A refined estimate of lk is obtained in the form

{circumflex over (l)}k(f)=[{circumflex over (l)}k]V−1.   (74)

Equation 73 is referred to below as the ESPRIT-based timing estimator (ETE).

In step 34, code detection is performed. From Equations 67 and 74, it can be seen that two distinct estimates {circumflex over (l)}k(F) and {circumflex over (l)}k(f) are available at the receiver for each code index lk. These estimates are now used to decide which codes are actually active in the considered ranging subchannel. For this purpose, we define two sets I(f)={{circumflex over (l)}k(f)}k=1K and I(F)={{circumflex over (l)}k(F)}k=1K and observe that, in the absence of any detection error, it should be I(f)=I(F)={lk}k=1K. Hence, for any code matrix CmεC the following detection strategy is utilized:




  
    
      
        
          C
          m
        
        ⁢
        
          
        
        ⁢
        is
        ⁢
        
          
        
        ⁢
        declared
        ⁢
        
          
        
        ⁢
        as
        ⁢
        
          {
          
            
              
                
                  detected
                
                
                  
                    
                      if
                      ⁢
                      
                        
                      
                      ⁢
                      m
                    
                    ∈
                    
                      
                        I
                        
                          (
                          f
                          )
                        
                      
                      ⋂
                      
                        I
                        
                          (
                          F
                          )
                        
                      
                    
                  
                
              
              
                
                  undetected
                
                
                  
                    
                      if
                      ⁢
                      
                        
                      
                      ⁢
                      m
                    
                    ∉
                    
                      
                        I
                        
                          (
                          f
                          )
                        
                      
                      ⋂
                      
                        I
                        
                          (
                          F
                          )
                        
                      
                    
                  
                
              
            
            .
          
        
      
    
    
      
        (
        75
        )
      
    
  




The BS notifies the detected codes through a downlink response message. Clearly, if an active RSS is using an undetected code, it must repeat its ranging process until success notification. Equation 75 is referred to below as the ESPRIT-based code detector (ECD).

In step 136, power levels are eventually estimated. Assume that the accuracy of the frequency and timing estimates provided by EFE and ETE is such that {circumflex over (ξ)}≈ξ and {circumflex over (η)}≈η. Then, the ML estimate of Ŝ(q)=[S1(q),S2(q), . . . ,SK(q)]T for q=0,1, . . . ,Q−1 based on the observation vectors Xm(q) is found to be




  
    
      
        
          
            
              S
              ^
            
            ⁡
            
              (
              q
              )
            
          
          =
          
            
              
                T
                
                  -
                  1
                
              
              (
              
                
                  ξ
                  ^
                
                ,
                
                  η
                  ^
                
              
              )
            
            ⁢
            
              
                X
                _
              
              (
              
                
                  q
                  ;
                  
                    ξ
                    ^
                  
                
                ,
                
                  η
                  ^
                
              
              )
            
          
        
        ⁢
        
          

        
        ⁢
        where
      
    
    
      
        (
        76
        )
      
    
  
  
    
      
        
          
            T
            (
            
              
                ξ
                ^
              
              ,
              
                η
                ^
              
            
            )
          
          =
          
            
              ∑
              
                m
                =
                0
              
              
                M
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  E
                  m
                  H
                
                (
                
                  
                    ξ
                    ^
                  
                  ,
                  
                    η
                    ^
                  
                
                )
              
              ⁢
              
                
                  E
                  m
                
                (
                
                  
                    ξ
                    ^
                  
                  ,
                  
                    η
                    ^
                  
                
                )
              
            
          
        
        ⁢
        
          

        
        ⁢
        and
      
    
    
      
        (
        77
        )
      
    
  
  
    
      
        
          
            X
            _
          
          (
          
            
              q
              ;
              
                ξ
                ^
              
            
            ,
            
              η
              ^
            
          
          )
        
        =
        
          
            ∑
            
              m
              =
              0
            
            
              M
              -
              1
            
          
          ⁢
          
            
          
          ⁢
          
            
              
                E
                m
                H
              
              (
              
                
                  ξ
                  ^
                
                ,
                
                  η
                  ^
                
              
              )
            
            ⁢
            
              
                X
                m
              
              ⁡
              
                (
                q
                )
              
            
          
        
      
    
    
      
        (
        78
        )
      
    
  




while Em(ξ,η) is a V×K matrix with entries

[Em({circumflex over (ξ)},{circumflex over (η)})]v+1,k=ej2π(mξk+νηk)0≦v≦V−1, 1≦k≦K.   (79)

From Equation 52 and letting Ŝ(q)=[S1(q),S2(q), . . . ,SK(q)]T yields

Ŝk(q)=Sk(q)+μk(q), k=1,2, . . . , K   (80)

where μk(q) is a zero-mean disturbance term with variance σμk2=σw2[T−1({circumflex over (ξ)},{circumflex over (η)})]k,k. The above equation indicates that an unbiased estimator of |Sk(q)|2 is given by |Sk(q)|2−σμk2. Finally, we define the average power of the kth RSS over the employed ranging subcarriers as





  
    
      
        
          P
          k
        
        =
        
          
            1
            Q
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                
                  
                  
                    
                      S
                      k
                    
                    ⁡
                    
                      (
                      q
                      )
                    
                  
                  
                
                2
              
              .
            
          
        
      
    
    
      
        (
        81
        )
      
    
  




This leads to the following ad hoc power estimator (AHPE):





  
    
      
        
          
            P
            ^
          
          k
        
        =
        
          
            1
            Q
          
          ⁢
          
            
              ∑
              
                q
                =
                0
              
              
                Q
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              [
              
                
                  
                    
                    
                      
                        
                          S
                          ^
                        
                        k
                      
                      ⁡
                      
                        (
                        q
                        )
                      
                    
                    
                  
                  2
                
                -
                
                  
                    σ
                    ^
                  
                  
                    μ
                    k
                  
                  2
                
              
              ]
            
          
        
      
    
    
      
        (
        82
        )
      
    
  




where {circumflex over (σ)}μk2={circumflex over (σ)}w2[T−1({circumflex over (ξ)},{circumflex over (η)})]k,k and {circumflex over (σ)}w2 is given by





  
    
      
        
          
            σ
            ^
          
          w
          2
        
        =
        
          
            1
            
              M
              -
              
                K
                ^
              
            
          
          ⁢
          
            
              ∑
              
                i
                =
                
                  
                    K
                    ^
                  
                  +
                  1
                
              
              M
            
            ⁢
            
              
            
            ⁢
            
              
                
                  λ
                  ^
                
                i
              
              .
            
          
        
      
    
    
      
        (
        83
        )
      
    
  




In step 138, a determination is made as to whether additional ranging subchannels exist. If a positive determination is made, step 126 occurs again, and the processing steps shown in FIG. 15 are repeated so that the additional ranging subchannels can be processed. Otherwise, step 140 occurs, wherein the estimated CFOs, timing errors, and power levels are broadcast to remote users. In the downlink response message, the BS will indicate only the detected codes while undetected RSSs must restart their ranging process.
FIGS. 16A-16C are graphs showing performance results of the ESPRIT-based embodiment of the present invention discussed above in connection with FIG. 15, obtained in software simulations. The investigated system had a total of N=1024 subcarriers over an uplink bandwidth of 3 MHz. The sampling period was Ts=0.33 μs, corresponding to a subcarrier distance of 1/(NTs)=2960 Hz. It was assumed that R=4 subchannels were available for IR. Each subchannel was divided into Q=16 tiles uniformly spaced over the signal spectrum at a distance of N/Q=64 subcarriers. The number of subcarriers in any tile was V=4 while M=4. The discrete-time CIRs had L=12 taps which were modeled as independent and circularly symmetric Gaussian random variables with zero means and an exponential power delay profile, i.e., E{|hk(1)|2}=σh2×exp(−1/12) for l=0,1, . . . , 11, where σh2 is chosen such that E{∥hk∥2}=1. Channels of different users were statistically independent of each other and were kept fixed over an entire time-slot. We considered a maximum propagation delay of θmax=204 sampling periods. Ranging blocks were preceded by a CP of length NG=256. The normalized CFOs are uniformly distributed over the interval [−Ω,Ω] and vary at each run. Recalling that the estimation range of EFE is |εk|<N/(2NT(M−1)), we set Ω≦0.1.
Performance testing began by investigating the performance of ECD in terms of probability of making an incorrect detection. FIG. 16A illustrates Pf as a function of SNR=1/σn2. The number of active RSSs was 2 and 3 while the maximum CFO was either Ω=0.05 or 0.1. Comparisons were made with the ranging scheme proposed by Fu, Li and Minn (FLM). The results of FIG. 16A indicate that ECD performs remarkably better than FLM.
FIG. 16B illustrates the RMSE of the frequency estimates obtained with EFE vs. SNR for K=2 or 3 and Ω=0.05 or 0.1. It can be seen that the accuracy of EFE is satisfactory at SNR values of practical interest. Moreover, EFE has virtually the same performance as K or Ω increase.
The performance of the timing estimators is measured in terms of probability of making a timing error, say P(ε). An error event is declared to occur whenever the estimate {circumflex over (θ)}k gives rise to IBI during the data section of the frame. This is tantamount to saying that the timing error {circumflex over (θ)}k−θk+(L−NG,D)/2 is larger than zero or smaller than L−NG,D−1, where NG,D is the CP length during the data transmission phase. For testing purposes, we set NG,D=32. FIG. 16C illustrates P(ε) vs. SNR as obtained with ETE and FLM. The number of active codes is K=2 or 3 while Ω=0.05 or 0.1. It can be seen that ETE provides much better results than FLM.
Having thus described the invention in detail it is to be understood that the foregoing description is not intended to limit the spirit or scope thereof. What is desired to be protected is set forth in the following claims.

## Claims

1. A method for initial ranging in a wireless communication system, comprising the steps of:
receiving at a base station a plurality of blocks transmitted to the base station from a plurality of remote user devices in wireless communication with the base station;
processing the plurality of blocks to extract a ranging subchannel from the plurality of blocks;
calculating a total number of active codes in the ranging subchannel by computing a correlation matrix and applying an eigendecomposition procedure to the correlation matrix, and using a minimum description length algorithm;
estimating carrier frequency offsets for all active devices of the plurality of remote user devices in a first processing step using the active codes;
estimating timing delays for all active devices of the plurality of remote user devices in a second processing step using the active codes;
estimating power levels for all active devices of the plurality of remote user devices in a third processing step using the active codes;
broadcasting the estimated carrier frequency offsets, timing delays, and power levels to the plurality of remote user devices; and
adjusting transmission parameters of the plurality of remote user devices in response to the broadcasted estimates to synchronize communication with the base station.

2. The method of claim 1, wherein the step of processing the plurality of blocks further comprises applying a discrete Fourier transformation algorithm to the plurality of blocks to extract the ranging subchannel.

3. The method of claim 1, wherein the step of identifying active codes in the ranging subchannel comprises decomposing an observation space into a signal subspace and a noise subspace.

4. The method of claim 3, wherein the step of identifying active codes in the ranging subchannel further comprises identifying active codes by applying a multiple signal classification algorithm to the signal subspace and the noise subspace.

5. The method of claim 1, wherein the step of estimating carrier frequency offsets further comprises estimating carrier frequency offsets by applying a multiple signal classification (MUSIC) algorithm to a signal subspace and a noise subspace.

6. The method of claim 1, wherein the step of estimating timing delays further comprises applying a maximum likelihood (ML) estimation algorithm to the ranging subchannel to produce an ML estimation vector.

7. The method of claim 6, wherein the step of estimating timing delays further comprises applying a least-squares estimation algorithm to the ML estimation vector to estimate timing delays.

8. The method of claim 1, wherein the step of estimating power levels comprises applying a least-squares estimation algorithm to a maximum likelihood (ML) estimation vector to estimate power levels.

9. The method of claim 1, wherein the step of estimating timing delays further comprises applying a simplified ML estimation algorithm having reduced computational complexity to the model to produce an ML estimation vector.

10. The method of claim 9, wherein the step of estimating timing delays further comprises applying a simplified least-squares estimation algorithm having reduced computational complexity to the ML estimation vector to estimate timing delays.

11. The method of claim 1, wherein the step of estimating power levels comprises applying a simplified least-squares estimation algorithm having reduced computational complexity to a maximum likelihood (ML) estimation vector to estimate power levels.

12. The method of claim 1, further comprising the step of processing all subcarriers and subchannels of the plurality of blocks to provide updated estimates of carrier frequency offsets, timing delays, and power levels for the plurality of remote user devices.

13. The method of claim 1, wherein the step of estimating carrier frequency offsets further comprises estimating carrier frequency offsets using an Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) algorithm.

14. A system for initial ranging in a wireless communication system, comprising:
an orthogonal frequency division multiple access (OFDMA) transceiver for receiving a plurality of blocks from a plurality of remote user devices in wireless communication with the OFDMA receiver;
an active code identification module programmed into the OFDMA transceiver for identifying active codes in a ranging subchannel of the plurality of; blocks, the active code identification module calculating a total number of active codes in the ranging subchannel by computing a correlation matrix and applying an eigendecomposition procedure to the correlation matrix, and using a minimum description length algorithm;
a carrier frequency offset estimation module for estimating carrier frequency offsets for all active devices of the plurality of remote user devices in a first processing step using the active codes;
a timing delay estimation module for estimating timing delays for all active devices of the plurality of remote user devices in a second processing step using the active codes; and
a power level estimation module for estimating power levels for all active devices of the plurality of remote user devices in a third processing step using the active codes,
wherein the OFDMA transceiver broadcasts the estimated carrier frequency offsets, timing delays, and power levels to the plurality of remote user devices so that the plurality of remote user devices can adjust transmission parameters in response to the broadcasted estimates to synchronize communication with the base station OFDMA transceiver.

15. The system of claim 14, wherein the modules process all subcarriers and subchannels of the plurality of blocks to provide updated estimates of carrier frequency offsets, timing delays, and power levels for the plurality of remote user devices, and the OFDMA transceiver broadcasts the updated estimates to the plurality of remote user devices.

16. A non-transitory computer-readable medium having computer-executable instructions for performing a method comprising the steps of:
extracting a ranging subchannel from a plurality of blocks received by a base station from a plurality of remote user devices in wireless communication with the base station;
calculating a total number of active codes in the ranging subchannel by determining a total number of received codes by creating a correlation matrix and applying an eigendecomposition procedure to the correlation matrix, and by applying a minimum description length algorithm to the received codes to calculate the total number of active codes;
estimating carrier frequency offsets for all active devices of the plurality of remote user devices in a first processing step using the active codes;
estimating timing delays for all active devices of the plurality of remote user devices in a second processing step using the active codes;
estimating power levels for all active devices of the plurality of remote user devices in a third processing step using the active codes; and
instructing the base station to broadcast the estimated carrier frequency offsets, timing delays, and power levels to the plurality of remote user devices, wherein the plurality of remote user device adjust transmission parameters in response to the broadcasted estimates to synchronize communication with the base station.

17. The computer-readable medium of claim 16, wherein the step of identifying active codes in the ranging subchannel comprises decomposing received codes in the ranging subchannel into a signal subspace and a noise subspace.

18. The computer-readable medium of claim 17, wherein the step of identifying active codes in the ranging subchannel further comprises identifying active codes by applying a multiple signal classification algorithm to the signal subspace and the noise subspace.

19. The computer-readable medium of claim 16, wherein the step of estimating carrier frequency offsets further comprises estimating carrier frequency offsets by applying a multiple signal classification algorithm to a signal subspace and a noise subspace.

20. The computer-readable medium of claim 16, wherein the step of estimating timing delays further comprises applying a maximum likelihood (ML) estimation algorithm to the ranging subchannel to produce an ML estimation vector.

21. The computer-readable medium of claim 20, wherein the step of estimating timing delays further comprises applying a least-squares estimation algorithm to the ML estimation vector to estimate timing delays.

22. The computer-readable medium of claim 16, wherein the step of estimating power levels comprises applying a least-squares estimation algorithm to a maximum likelihood (ML) estimation vector to estimate power levels.

23. The computer-readable medium of claim 16, wherein the step of estimating timing delays further comprises applying a simplified maximum likelihood (ML) estimation algorithm having reduced computational complexity to the model to produce an ML estimation vector.

24. The computer-readable medium of claim 23, wherein the step of estimating timing delays further comprises applying a simplified least-squares estimation algorithm having reduced computational complexity to the ML estimation vector to estimate timing delays.

25. The computer-readable medium of claim 16, wherein the step of estimating power levels comprises applying a simplified least-squares estimation algorithm having reduced computational complexity to a maximum likelihood (ML) estimation vector to estimate power levels.

26. The computer-readable medium of claim 16, further comprising the step of processing all subcarriers and subchannels of the plurality of blocks to provide updated estimates of carrier frequency offsets, timing delays, and power levels for the plurality of remote user devices.

27. The computer-readable medium of claim 16, wherein the step of estimating carrier frequency offsets further comprises estimating carrier frequency offsets using an Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) algorithm.