Doppler Spread Estimation in MIMO Frequency-selective Fading Channels

1 Doppler Spread Estimation in MIMO Frequenc y-selecti v e F ading Channels Mostafa Mohammadka rimi, Student Member , IEEE, Ebrahim Karami, Stude nt Me mber , IEEE, Octavia A. Dob re, Senior Member , IEEE, and Moe Z. W in F ellow , IEEE Abstract One of the main challenges in high-speed mobile communication s is the presence of large Doppler spreads. Thus, accurate estimation of maximum Doppler spread (MDS) plays an important role in improving the perform ance of the communication link. In this paper , we deriv e the data-aided (DA) and non-data-aided (ND A) Cramer-Rao lower bounds (CR L Bs) and maximum likelihood estimators (MLEs) for the MDS in multiple-input multi ple-output (MIMO) frequency-selecti ve fading channel. Moreov er , a low- complex ity NDA-mom ent-based estimator (MBE) is proposed. The proposed ND A-MBE relies on the second- and fourth-order moments of the receiv ed si gnal, which are employ ed to estimate the normalized squared autocorrelation function of the fading channel. T hen, the problem of MDS estimation is formulated as a non-linear regression problem, and the least-squares curve- fitting optimization technique is applied t o determine the estimate of the MDS. This is the fir st time in the l iterature when DA- and NDA -MDS esti mation is in vestigated for MIMO frequenc y-selectiv e fading channel. S imulation results sho w that there is no significant performance gap between the deri ved ND A-MLE and ND A-CRLB e ven w hen the observ ation windo w is relativ ely small. Furthermore, the significant reduced-complex ity in the NDA-M BE leads to lo w root-mean-squa re error (NRMSE) ov er a wide range of MDSs when the observ ation windo w is selected large enough. Index T erms Maximum Doppler spread, data-aided, non-data-aided, multiple-input multiple-output, frequency-selecti ve, Cramer-Rao lower bound (CRLB), fourth-order moment, autocorrelation, non-linear regression, maximum likelihood estimator (MLE). I . I N T R O D U C T I O N M AXIMUM DOPPLER SPREAD m easures the coh erence time, related to the rate of ch ange, of wireless com mu- nication ch annels. I ts k n owledge is imp ortant to d e sign efficient wireless commun ication systems for high-spee d vehicles [ 1 ]–[ 3 ]. In particular, accur ate estimation of the maximum Do ppler spread (MDS) is requ ired f o r the design of adaptive tra nsceiv ers, as well as in c e llular and smart an tenna systems [ 3 ]–[ 12 ]. For example, in the context of ad a p tiv e transceivers, system p arameters such as coding, mod ulation, and power a re ada p ted to the ch anges in the chann e l [ 4 ]–[ 7 ]. In cellular systems, h andoff is dictated b y the velocity of the mobile station, which is also dire c tly obtained from the Dop pler informa tio n. Knowledge of the rate of the ch a nnel change is also em ployed to reduce unnecessary hando ff; the hando ff is initiated based on the receiv ed power at the mobile station, and the optimu m window size for power estimation dep ends o n the MDS [ 7 ]–[ 10 ]. In the c o ntext of smart antenna systems, the MDS is used in the design of th e maximum likelihood (ML) M. Mohammadkari mi, E. Karami, and O. A. Dobre are with the Department of Electrical and Computer Engineeri ng, Memorial Univ ersity , St. John’ s, NL, Canada (e-mail: { m.mohamm adkarimi, ekarami, odobre } @mun.ca ). M. Z. W in is with the Laboratory for Information and Decision Systems (LIDS), Massachusett s Institute of T echn ology , Cambridge, MA, USA (e-mail: moewin@mit.edu ). 2 space-time tran scei vers [ 11 ], [ 1 2 ]. In addition , k nowledge o f MDS is required fo r chann el tracking and equalizatio n, as well as for the selection of the op timal interleaving length in wireless commu nication systems [ 13 ]. In general, parameter estimators can be categorized as: i) data-aided (D A), where the estimation relies on a pilo t or prea m ble sequence [ 14 ]–[ 18 ], ii) no n-data-aid ed (NDA ), where the estimation is perf ormed with no a priori kn owledge ab out the transmitted symbo ls [ 19 ]–[ 23 ], and iii) code- aided ( CA), where th e decod in g ga in is used via iterative feed back to enha nce the estimation perfo rmance of the desired pa r ameters [ 24 ]–[ 29 ]. W ith regard to the MDS estimatio n, the DA app roach o f ten pr ovides accurate estimates fo r slowly-varying ch annels b y employing a red uced nu mber o f pilo t symb o ls, wher eas this do es not hold fo r fast-varying chann els. In the latter case, the details of the channel variations cann ot be captured accura tely , and mor e pilots ar e requir ed, which r e su lts in increa sed overhead and reduced system capacity . There are five majo r classes of MDS estimato rs: ML-based, power spectral density (PSD) - based, lev el-cro ssing-rate (LCR)- based, covariance-based, and cyclostationar ity -based estimator s. T he ML-based estimator maxim izes th e likelihood fun ction, and, in gen eral, is asympto tically unbiased, achieving the Cramer-Rao lower bound (CRLB) [ 30 ]–[ 32 ]. Howe ver, maximum likelihood e stimator (MLE) for MDS suf fers from significant computationa l co mplexity . He n ce, dif fere nt mod ified lo w- complexity MLEs for M DS in single - input single-o utput (SISO) flat-fading chan nel were developed [ 33 ], [ 34 ]. W ith the PSD-based estimator s, some u nique features from the Doppler spectrum are obtain ed through th e sample per io dogram of the received sign al [ 35 ]. Cov ariance- based estimators extract the Doppler inform ation which exists in the sample auto-covariance of the received signa l [ 36 ]–[ 38 ]. LCR-based estimators rely on the numbe r of le vel crossings of the received sig n al statistics, which is p r oportio n al to the MDS [ 39 ]. The cyclostationarity-b ased estimators exploit the cyclostationarity o f the received signal [ 40 ]. Com paring with oth er MDS estimators, the advantage of the cyclostationarity- b ased estimato r s is the robustness to stationary noise and interferen ce. While the pro b lem of MDS estimatio n in SISO flat-fading chann e l h a s been extensively investigated in the liter ature [ 30 ]– [ 40 ], the MDS estimation in multiple-in put m u ltiple-outp ut (MIM O) freque ncy-selectiv e or in MIMO flat-fadin g ch annel has not been con siderably exp lo red. Fu r thermor e , DA -MDS estimation ha s mainly b een studied in the literature . T o th e best of ou r knowledge, only a few works have ad dressed MDS etimation in co njunction with multiple antenna sy stems. In [ 32 ], the au thors derived an asymptotic DA-MLE and D A-CRLB for joint MDS and n oise variance estimation in MIMO flat-fading ch a n nel. I n [ 40 ], the cyclic correlatio n (CC) of linearly modu lated sign als is exploited fo r the MDS estimatio n for single tran sm it an tenna scenarios. While bo th DA and ND A estimators ar e stud ied in [ 40 ], only f requency-flat fading an d sing le transmit antenn a are considered . In this paper, we inv estigate the pro blem o f MDS estimation in MI M O f requency-selective fading chann el for bo th DA an d ND A scenarios. Th e D A-CRLB, ND A-CRLB, DA-MLE, and NDA -MLE in MIMO frequ ency-selectiv e fading channe l are derived. In addition , a low-comp lexity ND A-moment- based estima to r (M BE) is propo sed. Th e pr oposed MBE relies on the second- an d fourth -order moments of the r eceiv ed signal a long w ith the lea st-sq uare (LS) curve-fitting optimization techn ique to estimate the normalized squared autocorr elation function (AF) and MDS o f the fading cha n nel. Sin ce the p roposed MBE is NDA , it removes th e need of pilots an d preambles used for D A-MDS estimation , and th us, it results in in c reased system 3 capacity . The N DA-MBE outperf orms the derived DA -MLE in the p resence of imp erfect time-freq uency synchroniza tion. Also, the MBE outperfo rms the ND A-CC estimator (CCE) in [ 40 ] and the DA low-complexity MLE in [ 33 ], [ 34 ] in SISO systems and under flat fading channels and in the presence of perfec t time-fre quency sy nchron ization. A. Contributions This paper brings the following or ig inal con tributions: • The D A- and ND A-CRLBs for MDS estimation in MIMO freq uency-selective fadin g chan nel are deriv ed; • The D A- and ND A-MLEs for MDS in MIMO fre quency-selective fading cha n nel are der i ved; • A low-complexity NDA-MBE is p r oposed. T he proposed estimator exhibits the f ollowing advantages: – lo wer comp utational complexity comp ared to the M LEs; – does not requir e time synchr onization; – is robust to th e c a rrier frequen cy offset; – increases system capacity; – does not requir e a priori knowledge o f no ise p ower , signal power , and cha n nel delay profile; – does not requir e a priori knowledge o f the nu m ber o f transmit antenna s; – removes the n e ed of joint par ameter estimation , such as carr ie r frequen cy offset, signal power , n o ise p ower , and channel delay pro file estimation; • The optimal combin ing me thod f or the ND A-MBE in case of multiple receiv e antennas is derived throu gh the bo otstrap technique . B. Notations Notation. Rand om variables are displayed in sans serif, uprig ht fo nts; their realizatio ns in serif, italic fon ts. V ectors and matrices ar e denoted by bold lowercase an d upper c ase letter s, respectively . For example, a rando m variable and its realiza tion are deno ted by x an d x ; a random vector and its realization are d enoted by x and x ; a rando m matrix and its realiza tio n are denoted by X and X , respectively . Th r ougho ut the paper, ( · ) ∗ is used f o r the complex co n jugate, ( · ) † is u sed for tran spose, | · | represents the absolute value operato r , ⌊·⌋ is the floor f unction, δ i,j denotes the Kron ecker delta f unction, n ! is th e factorial of n , E {·} is the statistical expectation , ˆ x is an estimate of x , an d det( A ) denotes th e determina nt of the matrix A . The re st of the paper is organized as follows: Section II de scr ibes the system mod el; Section III obtains th e DA - and ND A-CRLBs for MDS estimation in MI M O frequen cy-selecti ve fading; Section IV de r i ves the DA- and ND A-MLEs f or MDS in MIM O freq uency-selective fading ch a n nel; Section V introduces the proposed NDA-MBE for M DS; Section VI evaluates the co mputation al c o mplexity o f the der i ved estimators; Section VII pr e sents num erical results; and Section VII I conclu d es the paper . 4 I I . S Y S T E M M O D E L Let us consider a MIMO wireless co mmunicatio n system with n t transmit anten n as and n r receive antenn as, wher e the received signals are affected b y time-varying freq uency-selective Rayleigh fading and are corr upted b y additive white Gaussian noise. The discrete-time comp lex-valued b aseband signal at the n th receive antenna is expressed a s [ 41 ] r ( n ) k = n t X m =1 L X l =1 h ( mn ) k,l s ( m ) k − l + w ( n ) k k = 1 , ..., N , (1) where N is the numb er of observation symbols, L is the length of the channel impulse resp onse, s ( m ) k is the sym bol tran smitted from the m th antenna at time k , satisfy ing E  s ( m 1 ) k 1 ( s ( m 2 ) k 2 ) ∗  = σ 2 s m 1 δ m 1 ,m 2 δ k 1 ,k 2 , with σ 2 s m 1 being the tr a nsmit power of the m 1 th anten na, w ( n ) k is the complex-valued additive white Gaussian noise at the n th receive antenna a t time k , whose variance is σ 2 w n , and h ( mn ) k,l denotes the ze r o-mean complex-valued Gaussian fading process between the m th transmit and n th receiv e antennas fo r the l th tap of the fading channe l and at time k . I t is con sidered th at the chan nels for different ante n nas are indepen d ent, with the cross-co rrelation of the l 1 and l 2 taps gi ven by 1 E n h ( mn ) k,l 1  h ( mn ) k + u,l 2  ∗ o = σ 2 h ( mn ) ,l 1 J 0 (2 π f D T s u ) δ l 1 ,l 2 , (2) where J 0 ( · ) is the zero-or der Bessel func tion of the first kind, σ 2 h ( mn ) ,l 1 is the variance of the l 1 th tap between the m th transmit and n th receive antennas, T s denotes the symbol p eriod, and f D = v / λ = f c v / c represents th e MDS in H z , with v as the relativ e spee d between the tra nsmitter and receiver , λ as the wavelength, f c as the car rier freq uency , and c as the speed of light. I I I . C R L B F O R M D S E S T I M AT I O N In this section, the D A- and NDA-CR LB for MDS estimation in MIMO f requency-selective fadin g channel are d erived. A. D A-CRLB Let us co n sider s ( m ) k = s ( m ) k , m = 1 , 2 , · · · , n t , k = 1 , 2 , · · · , N − L + 1 , as employed pilots for DA-MDS estimation. The received sign al a t n th receive an tenna in ( 1 ) can be written as r ( n ) k = ¯ r ( n ) k + j ˘ r ( n ) k = n t X m =1 L X l =1 ¯ h ( mn ) k,l ¯ s ( m ) k − l − ˘ h ( mn ) k,l ˘ s ( m ) k − l + ¯ w ( n ) k + j n t X m =1 L X l =1 ¯ h ( mn ) k,l ˘ s ( m ) k − l + ˘ h ( mn ) k,l ¯ s ( m ) k − l + ˘ w ( n ) k ! , (3) where ¯ r ( n ) k , Re  r ( n ) k  , ˘ r ( n ) k , I m  r ( n ) k  , ¯ h ( mn ) k,l , Re  h ( mn ) k,l  , ˘ h ( m,n ) k,l , I m  h ( mn ) k,l  , ¯ s ( mn ) k − l , Re  s ( mn ) k − l  , and ˘ s ( mn ) k − l , Im  s ( mn ) k − l  . Let us define r ( n ) , h ¯ r ( n ) 1 ¯ r ( n ) 2 · · · ¯ r ( n ) N ˘ r ( n ) 1 ˘ r ( n ) 2 · · · ˘ r ( n ) N i † (4) 1 Here we consider the Jakes channel; it is worth noting that dif ferent parametric channel models can be also considere d. 5 and r , h r (1) † r (2) † · · · r ( n r ) † i † . (5) The e le m ents of the vector r ( n ) , n = 1 , 2 , · · · , n t , are linear combination s of the cor related Gaussian ra ndom variables as in ( 3 ). Thus, r , is a Gau ssian r andom vector with p robability d ensity function (PDF) given by p ( r | s ; θ ) = exp  − 1 2 r † Σ − 1 ( s , θ ) r  (2 π ) N n r det 1 2  Σ ( s , θ )  , (6) where Σ ( s , θ ) , E { rr † } , s ,  s (1) † s (2) † · · · s ( n t ) †  † , s ( m ) ,  ¯ s ( m ) 1 ¯ s ( m ) 2 · · · ¯ s ( m ) N − L +1 ˘ s ( m ) 1 ˘ s ( m ) 2 · · · ˘ s ( m ) N − L +1  † , and θ , [ ξ ϑ f D ] † is the parameter vector, with ξ , [ σ 2 w 1 · · · σ 2 w n r ] † (7a) ϑ , [ ϑ † 1 ϑ † 2 · · · ϑ † L ] † (7b) ϑ l , h σ 2 h (11) ,l · · · σ 2 h (1 n r ) ,l σ 2 h (21) ,l · · · (7c) σ 2 h (2 n r ) ,l · · · σ 2 h ( n t 1) ,l · · · σ 2 h ( n t n r ) ,l i † . Since r ( n 1 ) and r ( n 2 ) , n 1 6 = n 2 , a re unc orrelated rando m vectors, i.e . E  r ( n 1 ) r ( n 2 ) †  = 0 , the covariance m atrix of r , Σ ( s , θ ) , is b lock diagona l as Σ ( s , θ ) , E { rr † } =          Σ (1) Σ (2) . . . Σ ( n r )          , (8) where Σ ( n ) , E  r ( n ) r ( n ) †  . By employing ( 2 ), ( 3 ), and ( 4 ), using the fact the rea l an d imag inary pa rt of the fading tap are independ ent random v ariab les with E  | ¯ h ( mn ) k,l | 2  =  | ˘ h ( mn ) k,l | 2  = σ 2 h ( mn ) ,l / 2 , and after some alg ebra, the elements of the covariance matr ix Σ ( n ) , n ∈ { 1 , 2 , · · · , n r } , are obtaine d as E n ¯ r ( n ) k ¯ r ( n ) k + u o = E n ˘ r ( n ) k ˘ r ( n ) k + u o (9a) = 1 2 n t X m =1 L X l =1 σ 2 h ( mn ) ,l  ¯ s ( m ) k − l ¯ s ( m ) k + u − l + ˘ s ( m ) k − l ˘ s ( m ) k + u − l  J 0 (2 π f D T s u ) + σ 2 w n 2 δ u, 0 E n ¯ r ( n ) k ˘ r ( n ) k + u o = − E n ˘ r ( n ) k ¯ r ( n ) k + u o (9b) 1 2 n t X m =1 L X l =1 σ 2 h ( mn ) ,l  ¯ s ( m ) k − l ˘ s ( m ) k + u − l − ˘ s ( m ) k − l ¯ s ( m ) k + u − l  6 J 0 (2 π f D T s u ) . The Fisher inform ation matrix of the paramete r vector θ , I ( θ ) , for the zero-mea n Gau ssian observation vector in ( 6 ) is obtained as [ I ( θ )] ij , − E ( ∂ 2 ln p ( r | s ; θ ) ∂ θ i ∂ θ j ) (10) = 1 2 tr " Σ − 1 ( s , θ ) ∂ Σ ( s , θ ) ∂ θ i Σ − 1 ( s , θ ) ∂ Σ ( s , θ ) ∂ θ j # . For the MDS, f D , I ( f D ) , [ I ( θ )] xx , x = n t n r L + n r + 1 , and one obtains I ( f D ) = − E ( ∂ 2 ln p ( r | s ; θ ) ∂ f 2 D ) (11) = 1 2 tr "  Σ − 1 ( s , θ ) ∂ Σ ( s , θ ) ∂ f D  2 # , where ∂ Σ ( s , θ ) ∂ f D is ob ta in ed by r e placing J 0 (2 π f D T s u ) with − 2 π uT s J 1 (2 π f D T s u ) in Σ ( s , θ ) , where J 1 ( · ) is the Bessel function of the first kind. Finally , by employing ( 11 ), the D A-CRLB f or MDS estimation in MIMO f requency-selective fading ch annel is obtained as V ar( ˆ f D ) ≥ I − 1 ( f D ) = 1 1 2 tr "  Σ − 1 ( s , θ ) ∂ Σ ( s , θ ) ∂ f D  2 # . (12) I ( f D ) = − E ( ∂ 2 ln p ( r ; ϕ ) ∂ f 2 D ) = − 1 | M | N ′ n t Z x ∂ 2 ∂ f 2 D ln | M | N ′ n t X i =1 exp  − 1 2 x † Σ − 1 ( c h i i , ϕ ) x  det 1 2  Σ ( c h i i , ϕ )  ! | M | N ′ n t X q =1 exp  − 1 2 x † Σ − 1 ( c h q i , ϕ ) x  (2 π ) N n r det 1 2  Σ ( c h q i , ϕ )  d x . (19) B. ND A-CRLB Let u s conside r that the symb ols transmitted by each antenna are selected from a constellation with eleme n ts { c 1 c 2 · · · c | M | } , where 1 | M | P | M | i =1 | c i | 2 = 1 . The PDF of the received vector r for ND A-MDS estimation is expressed as p ( r ; ϕ ) = X c p ( r , c ; ϕ ) , (13) where c is the c onstellation vector as c ,  c (1) † c (2) † · · · c ( n t ) †  † , c ( m ) ,  ¯ c ( m ) 1 − L ¯ c ( m ) 2 − L · · · ¯ c ( m ) N − 1 ˘ c ( m ) 1 − L ˘ c ( m ) 2 − L · · · ˘ c ( m ) N − 1  † , c ( m ) k = ¯ c ( m ) k + j ˘ c ( m ) k is the co nstellation po int of the m th transmit antenn a at time k , and ϕ , [ β † ξ † ϑ † f D ] † with β , [ σ 2 s 1 σ 2 s 2 · · · σ 2 s n t ] † , and ξ and ϑ are given in ( 7 ). 7 By employing th e chain rule of proba bility and usin g p ( c = c h i i ) = | M | − N ′ n t , N ′ , N + L − 1 , one can write ( 13 ) as p ( r ; ϕ ) = X c p ( r , c ; ϕ ) = X c p ( c = c ) p ( r | c = c ; ϕ ) = 1 | M | N ′ n t | M | N ′ n t X i =1 p ( r | c = c h i i ; ϕ ) , (14) where c h i i represents the i th possible constellatio n vector a t the transmit-side. Similar to the D A-CRLB, p ( r | c = c h i i ; ϕ ) is Gaussian and p  r | c = c h i i ; ϕ  = exp  − 1 2 r † Σ − 1 ( c h i i , ϕ ) r  (2 π ) N n r det 1 2  Σ ( c h i i , ϕ )  , (15) where Σ ( c h i i , ϕ ) , E  r h i i r † h i i  is the covariance matrix of th e received vecto r r h i i giv en the co nstellation vector is c = c h i i , i = 1 , 2 , · · · , | M | N ′ n t . The 2 N n r × 2 N n r covariance matrix Σ ( c h i i , ϕ ) is block diag o nal a s in ( 8 ), where its diagon al elem e n ts, i.e., Σ ( n ) h i i , E n r ( n ) h i i r ( n ) h i i † o , n ∈ { 1 , 2 , · · · , n r } , are obtaine d as E n ¯ r ( n ) k, h i i ¯ r ( n ) k + u, h i i o = E n ˘ r ( n ) k, h i i ˘ r ( n ) k, h i i o (16a) = 1 2 n t X m =1 L X l =1 σ 2 h ( mn ) ,l σ 2 s m  ¯ c ( m ) k − l, h i i ¯ c ( m ) k + u − l, h i i + ˘ c ( m ) k − l, h i i ˘ c ( m ) k + u − l, h i i  J 0 (2 π f D T s u ) + σ 2 w n 2 δ u, 0 E n ¯ r ( n ) k, h i i ˘ r ( n ) k + u, h i i o = − E n ˘ r ( n ) k, h i i ¯ r ( n ) k, h i i o (16b) = 1 2 n t X m =1 L X l =1 σ 2 h ( mn ) ,l σ 2 s m  ¯ c ( m ) k − l, h i i ˘ c ( m ) k + u − l, h i i − ˘ c ( m ) k − l, h i i ¯ c ( m ) k + u − l, h i i  J 0 (2 π f D T s u ) . By substituting ( 15 ) into ( 14 ), one obtains p ( r ; ϕ ) = 1 | M | N ′ n t | M | N ′ n t X i =1 exp  − 1 2 r † Σ − 1 ( c h i i , ϕ ) r  (2 π ) N n r det 1 2  Σ ( c h i i , ϕ )  . (17) Finally , by employing ( 17 ), the ND A-CRLB for MDS estimation in MIMO frequency-selective fading channel is expressed as V ar( ˆ f D ) ≥ I − 1 ( f D ) = 1 − E n ∂ 2 ln p ( r ; ϕ ) ∂ f 2 D o , (18) where I ( f D ) is given in ( 19 ) on the top of this page, and R x , R x 1 R x 2 · · · R x (2 N n r ) . As seen, the r e is no an explicit expression for ( 19 ), and thus, for the CRLB in ( 18 ). Therefo re, num e rical meth o ds are u sed to solve ( 1 9 ) and ( 18 ). 8 κ ( n ) u = E n   r ( n ) k   2 | r ( n ) k + u   2 o = n t X m =1 L X l =1 E n   h ( mn ) k,l   2   h ( mn ) k + u,l   2 o σ 4 s m + n t X m 1 =1 n t X m 2 6 = m 1 L X l =1 E n   h ( m 1 n ) k,l   2   h ( m 2 n ) k + u,l   2 o σ 2 s m 1 σ 2 s m 2 + n t X m =1 L X l 1 =1 L X l 2 6 = l 1 E n   h ( mn ) k,l 1   2   h ( mn ) k + u,l 2   2 o σ 4 s m + n t X m 1 =1 n t X m 2 6 = m 1 L X l 1 =1 L X l 2 6 = l 1 E n   h ( m 1 n ) k,l 1   2   h ( m 2 n ) k + u,l 2   2 o σ 2 s m 1 σ 2 s m 2 + 2 σ 2 w n n t X m =1 L X l =1 E n   h ( mn ) k,l   2 o σ 2 s m + σ 4 w n , u ≥ L. (28) I V . M L E S T I M AT I O N F O R M D S In this section, we derive th e D A- and ND A-MLEs for MDS in MIMO fr equency-selective fading chann el. A. D A-MLE for MDS The D A-MLE for f D is obtain e d as ˆ f D = arg max f D p ( r | s ; θ ) , (20) where p ( r | s ; θ ) is giv en in ( 6 ). Since p ( r | s ; θ ) is a differentiable fu n ction, the DA-MLE for f D is obtained from ∂ ln p ( r | s ; θ ) ∂ f D = 0 . (21) By substituting ( 6 ) in to ( 21 ) and after some mathematical manipu la tio ns, one obtains ∂ ln p ( r | s ; θ ) ∂ f D = − 1 2 tr " Σ − 1 ( s , θ ) ∂ Σ ( s , θ ) ∂ f D # (22) + 1 2 r † Σ − 1 ( s , θ ) ∂ Σ ( s , θ ) ∂ f D Σ − 1 ( s , θ ) r . As seen in ( 22 ), there is n o closed-for m solutio n for ( 21 ). Thu s, numerical methods need to be used to obtain solution . By employing th e Fisher -scor ing m e thod [ 42 ], 2 the solution of ( 22 ) can be iteratively obtaine d as ˆ f [ t +1] D = ˆ f [ t ] D + I − 1 ( f D ) ∂ ln p ( r | s ; θ ) ∂ f D      f D = ˆ f [ t ] D , (23) where I ( f D ) and ∂ ln p ( r | s ; θ ) ∂ f D are giv en in ( 11 ) and ( 22 ), respectively . B. ND A-MLE for MDS Similar to the D A-MLE, the NDA -MLE for MDS is ob tain ed fr om ˆ f D = arg max f D p ( r ; ϕ ) , (24) 2 The Fisher-scorin g method replac es the Hessian matrix in the Ne wtown -Raphson method with the negat iv e of the Fisher information matrix [ 43 ]. 9 where p ( r ; ϕ ) is g i ven in ( 17 ). Since p ( r ; ϕ ) is a linear combin ation of differentiable fu nctions, the NDA-MLE for f D is obtained from ∂ ln p ( r ; ϕ ) ∂ f D = 0 . (25) By substituting ( 17 ) into ( 25 ) and a fter some algebra , one obtains | M | N ′ n t X i =1 ( r † Σ − 1 ( c h i i , ϕ ) ∂ Σ ( c h i i , ϕ ) ∂ f D Σ − 1 ( c h i i , ϕ ) r det 1 2 Σ ( c h i i , ϕ ) − tr h Σ − 1 ( c h i i , ϕ ) ∂ Σ ( c h i i , ϕ ) ∂ f D i det 1 2 Σ ( c h i i , ϕ ) ) = 0 (26) Similar to the D A-MLE, there is no closed -form solu tion fo r ( 26 ); thus, numer ic a l methods are used to solve ( 26 ). V . N DA - M O M E N T - B A S E D ( M B ) E S T I M A T I O N O F M D S In th is section, we propo se an ND A-MB MDS estimator for multiple inp ut single o utput (MISO) systems under frequency- selecti ve Rayleigh fading chan nel by e m ploying the fourth- order moment of the rec ei ved sign al. Then, an extension of the propo sed estimator to the MIMO systems is provid ed. A. ND A-MBE for MDS in MIS O S ystems Let us assume that th e par ameter vector ϕ = [ β † ξ † ϑ † f D ] † is u nknown at the receive-side. The statistical MB approach enables u s to pr opose an NDA-MBE to estimate f D without any pr iori knowledge of β , ξ , and ϑ . Let u s consider th e fourth - order two-conjug ate momen t of the received signa l a t the n th receive ante n na, d efined as κ ( n ) u ∆ = E n   r ( n ) k   2   r ( n ) k + u   2 o . (27) W ith the transmitted symbols, s ( m ) k , m = 1 , ..., n t being ind ependen t, drawn from sym metric comp lex-valued co nstellation points, 3 and with u ≥ L , κ ( n ) u is expressed as in ( 28 ) at the top of this page (see Append ix I fo r pro of). By employing th e first-order autoregressive mod e l of the Rayleigh fading chann el, on e can wr ite [ 45 ], [ 46 ] h ( mn ) k,l = Ψ u h ( mn ) k + u,l + v ( mn ) k,l , (29) where Ψ u , J 0 (2 π f D T s u ) an d v ( mn ) k,l is a zero -mean complex-valued Gaussian white p r ocess with variance E {| v ( mn ) k,l | 2 } = (1 − | Ψ u | 2 ) σ 2 h ( mn ) ,l , which is indepen dent of h ( mn ) k,l . By using ( 29 ) and exploiting the p roperty of a complex-valued Gau ssian rando m variable z ∼ N c  0 , σ 2 z  that E {| z | 2 n } = n ! σ 2 n z [ 47 ], one obtains E n   h ( mn ) k,l   2 | h ( mn ) k + u,l   2 o (30) = E     ( J 0 (2 π f d T s u ) h ( mn ) k + u,l + v ( mn ) k,l )    2   h ( mn ) k + u,l   2  3 E  ( s ( m ) k ) 2  = 0 for M -ary phase-shift -keyi ng (PSK) and quadrature amplit ude modulation (QAM), M > 2 [ 44 ]. 10 = J 2 0 (2 π f d T s u ) E n   h ( mn ) k + u,l   4 o + E n   v ( mn ) k,l   2 | h ( mn ) k + u,l   2 o + J 0 (2 π f d T s u ) E n h ( mn ) k + u,l   h ( mn ) k + u,l   2 ( v ( mn ) k,l ) ∗ o + J 0 (2 π f d T s u ) E n ( h ( mn ) k + u,l ) ∗   h ( mn ) k + u,l   2 ( v ( mn ) k,l ) o = 2 J 2 0 (2 π f d T s u ) σ 4 h ( mn ) ,l +  1 − J 2 0 (2 π f d T s u )  σ 4 h ( mn ) ,l =  1 + J 2 0 (2 π f d T s u )  σ 4 h ( mn ) ,l m = 1 , ...n t , l = 1 , ..., L. W ith the channel taps l 1 and l 2 being uncorrela ted for eac h tran smit a n tenna, i.e., E  h ( mn ) k,l 1 ( h ( mn ) k,l 2 ) ∗  = σ 2 h ( mn ) ,l 1 δ l 1 ,l 2 and employing E n   h ( m 1 n ) k,l 1   2   h ( m 2 n ) k + u,l 2   2 o (31) = σ 2 h ( m 1 n ) ,l 1 σ 2 h ( m 2 n ) ,l 2 h (1 − δ l 1 ,l 2 )(1 − δ m 1 ,m 2 ) + δ l 1 ,l 2 (1 − δ m 1 ,m 2 ) + (1 − δ l 1 ,l 2 ) δ m 1 ,m 2 i + σ 4 h ( m 1 n ) ,l 1  1 + J 2 0 (2 π f d T s u )  δ m 1 ,m 2 δ l 1 ,l 2 , one can write ( 28 ) as κ ( n ) u = n t X m =1 L X l =1 σ 4 h ( mn ) ,l σ 4 s m  1 + J 2 0 (2 π f d T s u )  (32) + n t X m 1 =1 n t X m 2 6 = m 1 L X l =1 σ 2 h ( m 1 n ) ,l σ 2 h ( m 2 n ) ,l σ 2 s m 1 σ 2 s m 2 + n t X m =1 L X l 1 =1 L X l 2 6 = l 1 σ 2 h ( mn ) ,l 1 σ 2 h ( mn ) ,l 2 σ 4 s m + n t X m 1 =1 n t X m 2 6 = m 1 L X l 1 =1 L X l 2 6 = l 1 σ 2 h ( m 1 n ) ,l 1 σ 2 h ( m 2 n ) ,l 2 σ 2 s m 1 σ 2 s m 2 + 2 σ 2 w n n t X m =1 L X l =1 σ 2 h ( mn ) ,l σ 2 s m + σ 4 w n . ˆ f [ t +1] D = ˆ f [ t ] D − M max P u = U min 8 π T s u  ˆ Ψ ( n ) u − J 2 0 (2 π f [ t ] D T s u )  J 0  2 π f [ t ] D T s u  J 1  2 π f [ t ] D T s u  ∂ 2 ∂ f 2 D U max P u = U min  ˆ Ψ ( n ) u − J 2 0 (2 π f D T s u )  2   f D = f [ t ] D (42) ∂ 2 ∂ f 2 D U max X u = U min  ˆ Ψ u − J 2 0 (2 π f D T s u )  2   f D = f [ t ] D = M max X u = U min ( 32 π 2 T 2 s u 2 J 2 0  2 π f [ t ] D T s u  J 2 1  2 π f [ t ] D T s u  + 8 π T s u 2 π T s u  J 2 0  2 π f [ t ] D T s u  − J 2 1  2 π f [ t ] D T s u   − J 0  2 π f [ t ] D T s u  J 1  2 π f [ t ] D T s u  f [ t ] D !  ˆ Ψ ( n ) u − J 2 0 (2 π f [ t ] D T s u )  ) 11 Further, let us consid er the second-o rder momen t of th e received sig nal, i.e., µ ( n ) 2 ∆ = E {| r ( n ) k | 2 } . By using ( 1 ), it can b e easily shown that µ ( n ) 2 = n t X m =1 L X l =1 σ 2 h ( mn ) ,l σ 2 s m + σ 2 w n . (33) By employing ( 32 ) and ( 33 ), one obtain s the nor malized sq uared AF of the fading ch annel a s (see Appendix II for proof ) Ψ u ∆ = J 2 0 (2 π f d T s u ) = η ( n )  κ ( n ) u −  µ ( n ) 2  2  , (34) where η ( n ) = 1 / n t P m =1 L P l =1 σ 4 h ( mn ) ,l σ 4 s m . For n on-con stant modulu s constellation s, η ( n ) is expressed in terms o f µ ( n ) 4 ∆ = E  | r ( n ) k | 4  and µ ( n ) 2 as (see Appe n dix III for proof ) η ( n ) = 2(Ω s − 1) µ ( n ) 4 − 2  µ ( n ) 2  2 , (35) where Ω s = 1 | M | P | M | i =1 | c i | 4 is a constan t, and 1 < Ω s ≤ 2 . 4 Finally , substitutin g ( 35 ) in to ( 34 ) y ields Ψ u ∆ = 2(Ω s − 1) κ ( n ) u −  µ ( n ) 2  2 µ ( n ) 4 − 2  µ ( n ) 2  2 . (36) As seen, the normalized squared AF o f the fading channe l is expressed as a n on-linear functio n of the µ ( n ) 2 , µ ( n ) 4 , and κ ( n ) u . In practice, statistical momen ts are estimated by time averages of the r e c ei ved signa l. For ( 36 ), the following e stima to rs of the moments are employed ˆ µ ( n ) 2 = 1 N N X k =1   r ( n ) k   2 (37) ˆ µ ( n ) 4 = 1 N N X k =1 | r ( n ) k   4 ˆ κ ( n ) u = 1 N − u N − u X k =1   r ( n ) k   2   r ( n ) k + u   2 , where u ≥ L > 0 . By substituting the correspo nding estimators in ( 36 ), the estimate of the normalize d squared AF is given as ˆ Ψ ( n ) u ∆ = 2 (Ω s − 1) ˆ κ ( n ) u −  ˆ µ ( n ) 2  2 ˆ µ ( n ) 4 − 2  ˆ µ ( n ) 2  2 . (38) Now , based on ( 34 ) and ( 38 ), the problem of MDS estimation ca n be formulated as a non-line a r regression pro blem. G iven the estimated n o rmalized squ ared AF , ˆ Ψ ( n ) u , the non-linear regression m odel assumes th a t the relationship between ˆ Ψ ( n ) u and 4 For 16-QAM, 64-QAM, and complex-v alue d zero-mean Gaussian signals, Ω s is 1.32, 1.38, and 2, respecti vely [ 44 ]. 12 Ψ u is modeled through a disturb ance term or err or variable ǫ ( n ) u as [ 48 ], [ 49 ] ˆ Ψ ( n ) u = Ψ u + ǫ ( n ) u (39) = J 2 0 (2 π f D T s u ) + ǫ ( n ) u , u = U min , . . . , U max , where U min and U max are the maximu m and minim um d elay lags, r espectiv ely . T o so lve the no n-linear regression prob lem in ( 39 ), the LS cu rve-fitting o p timization tech n ique is em ployed. Based on the LS cur ve-fitting optimization, the estimate of f D , i.e., ˆ f D , is ob tained thro ugh minimizin g the sum of the squared residuals (SSR) as [ 49 ] minimize f D U max X u = U min  ˆ Ψ ( n ) u − J 2 0 (2 π f D T s u )  2 subject to f l ≤ f D ≤ f h , (40) where f l and f h are the m inimum an d maximum possible MDSs, respectively . T o o btain ˆ f D , we consider the de r i vati ve o f the SSR with respect to f D and set it equa l to zero as follows: M max X u = U min 8 π T s u  ˆ Ψ ( n ) u − J 2 0 (2 π f D T s u )  (41) J 0 (2 π f D T s u ) J 1 (2 π f D T s u ) = 0 . As seen , for the non- linear regression, the deriv ative in ( 41 ) is a function o f f D . Thus, an explicit solu tion f o r ˆ f D cannot be obtained. Howe ver , num erical me th ods [ 50 ] can b e employed to solve the LS cu rve-fitting optim iz a tion pr oblem in ( 4 0 ). By employing the Newton-Raphson method, ˆ f D can be iterativ ely obtained as it is shown in ( 42 ) at th e top o f next p age. The main pro b lem with the Newton-Raphso n m ethod is that it suffers f rom the co n vergence pro b lem [ 43 ]. Since the p arameter space for th e MDS estimation is on e-dimensio n al, the g rid search m ethod can be employed, which en sures the glob al optima lity of the solution . W ith the grid search me thod, the parame te r space, i.e. , [ f l , f h ] is discretized as a grid with step size δ , a n d the value which minimizes SSR is considered as the e stima ted f D . This pro cedure can be perfo rmed in two steps, includ ing a r ough estimate of the MDS, ˆ f (r) D , by cho osing a larger step size ∆ fo llowed b y a fine estimate, ˆ f (s) D , through small grid step size δ aroun d the ro u gh estimate, i.e.,  ˆ f (r) D − ∆ , ˆ f (r) D + ∆  . A fo rmal description of th e propo sed ND A-MBE f or MDS in MISO frequen cy-selecti ve c h annel is p resented in Algorithm 1 . It is worth noting that f D can be estimated by using a d ownsampled version of ˆ Ψ ( n ) u . For the case of unifo rm downsampling, i.e., u = ℓu s , the SSR is given as N la − 1 X ℓ =0  ˆ Ψ ( n ) U min + ℓu s − Ψ U min + ℓu s  2 , (43) where u s is the downsampling per iod expr e ssed in delay lags, N la is the numb er of delay lag, Ψ ℓu s = J 2 0  2 π f D T s ( U min + ℓu s )  , (44) 13 Algorithm 1 : ND A-MBE for MDS in MISO systems 1: Set f l , f h , ∆ , and δ 2: Acquire the measurem ents  r ( n ) k  N k =1 3: Estimate the statistics ˆ µ ( n ) 2 , ˆ µ ( n ) 4 , and ˆ κ ( n ) u , by employing ( 37 ) 4: Compute ˆ Ψ ( n ) u , ∀ u ∈  U min , . . . , U max  by using ( 38 ) 5: Obtain ˆ f (r) D by solving the minimization prob lem in ( 40 ) thr ough the grid search method with grid step size ∆ 6: Obtain ˆ f (s) D by solving the m inimization problem in ( 40 ) through the gr id search method over  ˆ f (r) D − ∆ , ˆ f (r) D + ∆  with grid step size δ 7: ˆ f D = ˆ f (s) D and ˆ Ψ ( n ) U min + ℓu s = 2 (Ω s − 1) ˆ κ ( n ) U min + ℓu s −  ˆ µ ( n ) 2  2 ˆ µ ( n ) 4 − 2  ˆ µ ( n ) 2  2 . (45) The downsampled version of ˆ Ψ ( n ) u is usually employed for the roug h MDS estimation, where ∆ is a large value. B. ND A-MBE for MDS in MIMO S ystems The perf o rmance of the propo sed ND A-MBE fo r MDS in MISO system can be imp r oved when employin g multiple rec e i ve antennas due to the spatial diversity , b y combining the estimated normalized squared AFs, ˆ Ψ ( n ) u , n = 1 , ..., n r as e Ψ u = n r X n =1 λ ( n ) u ˆ Ψ ( n ) u , (46) where Λ u ,  λ (1) u λ (2) u · · · λ ( n r ) u  † , with P n r n =1 λ ( n ) u = 1 , is the weighting vector . Let u s defin e ˆ Ψ u ,  ˆ Ψ (1) u ˆ Ψ (2) u · · · ˆ Ψ ( n r ) u  † . The mean square error (MSE) of the comb ined nor malized sq u ared AF in ( 46 ) is expressed as E n  e Ψ u − Ψ u  2 o = Λ † u C u Λ u +  Λ † u µ u − Ψ u  2 , (47) where C u , E n  ˆ Ψ u − µ u  ˆ Ψ u − µ u  † o and µ u , E  ˆ Ψ u  . By e m ploying the method o f Lagran ge mu ltipliers, the optimal weighting vector Λ op u in ( 47 ) in terms of minim um MSE is obtained as Λ op u =  1 † y u  − 1 y u , (48) where y u ,  C u + ( µ u − Ψ u 1 )( µ u − Ψ u 1 ) †  − 1 1 and 1 , [1 1 · · · 1] † is an n r -dimension al vector of ones. As seen , the o ptimal weightin g vector, Λ op u , in ( 48 ) dep ends on th e true value o f MDS, i.e. , f D , thro ugh the true normalize d squared AF , Ψ u , in y u . T o obtain the optimal weightin g vector , the mean vector µ u and covariance m atrix C u are requir ed to be estimated fr om the received symbo ls. One approac h is bo otstrapping [ 51 ]–[ 53 ]. The b ootstrap method suggests to re-sample the empirica l joint cumulative d istribution functio n (CDF) o f ˆ Ψ u to estimate µ u and C u as summar ized in Algor ithm 2 . 5 5 Since ˆ Ψ ( n ) u , n = 1 , ..., n r are uncorrelate d random va riables, C u is a diagonal matrix. Thus, only the diagonal elements of ˆ C u are employed to obtain the optimal weighting vector . 14 Algorithm 2 : Bootstrap Algorithm for Optimal Combinin g 1: Set N B 2: for t = 1 , 2 , · · · , N B do 3: Draw a ra ndom sample of size N , with r eplacement, from X , { 1 , 2 , · · · , N } and nam e it X ⋆ 4: for n = 1 , 2 , · · · , n r do ˆ Ψ ( n ) ⋆ u [ t ] = 1 N − u P k ∈ X ⋆   r ( n ) k   2   r ( n ) k + u   2 −  1 N P k ∈ X ⋆ | r ( n ) k   2  2 1 N P k ∈ X ⋆ | r ( n ) k   4 − 2  1 N P k ∈ X ⋆ | r ( n ) k   2  2 5: end for 6: ˆ Ψ ⋆ u [ t ] , 2(Ω s − 1)  ˆ Ψ (1) ⋆ u ˆ Ψ (2) ⋆ u · · · ˆ Ψ ( n r ) ⋆ u  † 7: end for 8: Γ u = h ˆ Ψ ⋆ u [1] ˆ Ψ ⋆ u [2] · · · ˆ Ψ ⋆ u [ N B ] i 9: ˆ µ u = 1 N B P N B t =1 ˆ Ψ ⋆ u [ t ] 10: ˆ C u = 1 N B − 1 ( Γ u − ˆ µ u 1 † )( Γ u − ˆ µ u 1 † ) † 1 100 200 300 400 500 600 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Fig. 1. Illustra tion of the non-linear LS regression for the uniformly sampled normalize d squared AF for f D T s = 0 . 02 and f D T s = 0 . 005 , with n t = 1 , n r = 2 , L = 1 , u s = 2 , and at γ = 10 dB. As seen in Alg o rithm 2 , the optimal weig hting vector for e ach d e la y lag u is derived at the expense o f higher computa tio nal complexity . In order to avoid this com putational co mplexity , the sub optimal equ al weight co mbining meth od ca n be em ployed as e Ψ u = 1 n r n r X n =1 ˆ Ψ ( n ) u . (49) Fig. 1 shows h ow e Ψ u fits Ψ u throug h the equal weight combinin g in ( 49 ) for f d T s = 0 . 02 a n d f d T s = 0 . 00 5 with n t = 1 , n r = 2 , L = 1 , u s = 2 , an d at γ = 10 dB. Finally , similar to the MISO scena r io, the prob le m of MDS estimation for multip le receive antennas is fo rmulated as non - linear regression pr oblem in ( 39 ) for e Ψ u . A formal descrip tion of the prop osed NDA-MBE for MDS in MI MO frqu e ncy-selectiv e channel is presented in Algor ithm 3 . 15 Algorithm 3 : ND A-MBE for MDS in MIMO systems 1: Set f l , f h , ∆ , and δ 2: Acquire the measurem ents { r ( n ) k } N k =1 , ∀ n ∈  1 , . . . , n r  3: Estimate the statistics ˆ µ ( n ) 2 , ˆ µ ( n ) 4 , and ˆ κ ( n ) u , by employing ( 37 ) for { r ( n ) k } N k =1 , ∀ n ∈  1 , . . . , n r  4: Compute ˆ Ψ ( n ) u , ∀ u ∈  U min , . . . , U max  , ∀ n ∈  1 , . . . , n r  , by using ( 38 ) 5: Compute e Ψ u , ∀ u ∈  U min , . . . , U max  , by using ( 49 ) 6: Obtain ˆ f (r) D by solving the minimization prob lem in ( 40 ) for e Ψ u throug h the grid search meth o d with step size ∆ 7: Obtain ˆ f (s) D by solv ing the minimization prob lem in ( 40 ) for e Ψ u via the grid search method over  ˆ f (r) D − ∆ , ˆ f (r) D + ∆  with step size δ 8: ˆ f D = ˆ f (s) D T abl e I N U M B E R O F R E A L A D D I T I O N S , R E A L M U LT I P L I C A T I O N S , A N D C O M P L E X I T Y O R D E R O F T H E P RO P O S E D N DA - M B E. Algorithm Real additions Real multi plications Order MISO  N + 2 N g − ( U max + U min ) 2  N la + 3 N − N g − 1  N + N g − ( U max + U min ) 2 + 2  N la + 3 N + 4 O ( N ) MIMO n r   N − ( U max + U min ) 2  N la + 3 N − 1  + (2 N la − 1) N g n r   N − ( U max + U min ) 2 + 2  N la + 3 N + 4  + N g N la + 1 O ( N ) C. Semi-b lind ND A- MBE The p roposed ND A-MBE for MISO and MIMO systems do n ot req u ire knowledge of the p a rameter vector ϕ = [ β † ξ † ϑ † f D ] † . In other words, th e propo sed NDA-MBE in section V -A and V -B are blind. For the scenarios in which the variance of the ad ditiv e noise can be a ccurately estimated at the receive an tennas, i.e., ξ is k nown, a semi-b lind NDA -MBE for the case of SISO transmission and flat- fading channel, i.e., n t = 1 and L = 1 , can be pro posed. In this case, for the n th receive antennas, one can easily obtain 6 µ ( n ) 2 = σ 2 h n σ 2 s + σ 2 w n (50) and η ( n ) = ( σ 4 h n σ 4 s ) − 1 = 1  µ ( n ) 2 − σ 2 w n  2 . (51) By using ( 34 ), ( 50 ) and ( 51 ), and by rep lacing the statistical momen ts and the noise variance with their corresp onding estimates, one obtains ˆ Ψ ( n ) u = ˆ κ ( n ) u −  ˆ µ ( n ) 2  2  ˆ µ ( n ) 2 − ˆ σ 2 w  2 , (52) where ˆ σ 2 w n is th e estimate of the noise variance, and ˆ κ ( n ) u and ˆ µ ( n ) 2 are given in ( 37 ). Clearly , similar to the SISO transmission, the optimal and suboptima l comb in ing m ethods f or the multiple receiv e antenn a s can be employed , as well. 6 The index of transmit antenna , i.e., m = 1 and the inde x of channel tap, i.e., l = 1 is dropped. 16 256 512 1024 2048 10 4 10 6 10 8 10 10 10 12 Fig. 2. Computat ional complex ity comparison of the proposed NDA-MBE, the low-compl exity D A-ML E in [ 33 ], [ 34 ], and the D A-COMA T estimator in [ 38 ]. V I . C O M P L E X I T Y A NA L Y S I S By employing the two steps grid search method to solve the optimizatio n prob lem in ( 40 ), th e number o f real additio ns and multiplications emp loyed in th e pro posed ND A-MBE is shown in T able I , wh ere N la is th e numb er of delay lag, N g , ( N g 1 + N g 2 ) , and N g 1 and N g 2 are th e num ber of grid points used for the ro ugh and fine estimation, respe cti vely . As seen, th e propo sed ND A-MBE exhibits a complexity orde r of O ( N ) . It sho uld be mentio ned that the complexity order o f the der i ved D A-MLE and ND A-MLE are O ( N 3 ) and O ( | M | N ′ n t ) , respectively . Fig. 2 co mpares the total numb er of operatio ns used b y the pro posed NDA -MBE with the low-complexity DA-MLE in [ 33 ], [ 3 4 ] and the D A-COMA T estimator in [ 38 ]. As seen, the propo sed NDA-MBE exhibits sign ificantly lower comp utational complexity comp ared to the DA-COMA T [ 38 ] and the low-comp lexity D A-MLE in [ 33 ], [ 34 ]. This substantial reduced - complexity en a b les th e propo sed MBE to exhib it goo d per forman c e in the ND A scenario s, wh ere th e o bservation wind ow c a n be selected large eno ugh. V I I . S I M U L AT I O N R E S U LT S In th is section, we examine th e per f ormance of the propo sed N DA-MBE, as well as th e der i ved DA -MLE a n d DA -CRLB for MDS in MIMO freq u ency-selective fading chann el thro ugh several simulation experim ents. A. Simulation Setup W e con sider a MIMO system e m ploying spatial m ultiplexing, with carr ie r f requency f c = 2 . 4 GHz. Unless oth erwise mentioned , n t = 2 , n r = 2 , T s = 10 µs , N = 1 0 5 , a nd the m o dulation is 64- QAM. The d elay pr ofile of the Rayleigh fading ch annel is σ 2 h ( mn ) ,l = β exp ( − l rms l / L ) , wh ere β is a norm alization factor, i.e., β P l ( − l rms l / L ) = 1 , with L = 5 and l rms = L/ 4 as the max im um and RMS delay spread o f the chan nel, respectively . The parame ters for the downsampled LS 17 curve-fitting optimization are U min = L , U max = ⌊ N 10 ⌋ , and u s = 10 . The add iti ve white noise was mod eled as a comp lex- valued Gaussian ra n dom variable with zero- mean and variance σ 2 w for each receive antenn as. Wit ho ut loss of gen erality , it was assumed that σ 2 s m = 1 / ( n t n r ) , m = 1 , 2 , .., n t , and thu s, the average SNR was defined as γ = 10 log  1  n r σ 2 w  . Un less otherwise men tioned, the perf ormance o f the MDS estimators was p resented in terms of n ormalized root-m ean-squar e error (NRMSE), i.e., E { ( ˆ f D T s − f D T s ) 2 } 1 / 2 /f D T s , ob tained from 100 0 Mo nte Carlo trials for each f D T s ∈ [10 − 3 , 1 8 × 10 − 3 ] , with the search step size ∆ = 10 Hz and δ = 0 . 5 Hz, resp ectiv ely . B. Simulation Results Fig. 3 shows the distributions of the estimated ˆ f D by the p roposed ND A-MBE for different MDSs, f D = 1 000 Hz an d f D = 10 0 Hz, with n t = 2 , n r = 2 , and a t γ = 10 d B. As seen, the distributions are not symmetric arou nd th eir me a n values; h ence, this leads to bias in MDS estimation . Furtherm ore, Fig . 4 illustrates E { ˆ f D /f D } versus f D for γ = 10 dB and γ = 20 dB. As seen , the pro posed ND A-MBE is near ly unbiased, i.e., E { ˆ f D } ≈ f D over a wide range of MDS. This can be (a) f D = 1000 Hz (b) f D = 100 Hz Fig. 3. Distributi on of the estimated ˆ f D for the proposed ND A -MBE for n t = 2 , n r = 2 , and at γ = 10 dB. 18 200 600 1000 1400 1800 0.95 0.96 0.97 0.98 0.99 1 Fig. 4. The mean value of the estimated MDS by NDA-MBE for va rious SNR va lues for n t = 2 and n r = 2 . 1 2 4 6 8 10 12 14 16 18 10 -3 10 -3 10 -2 10 -1 10 0 Fig. 5. The NRMSE of the proposed ND A-MBE versus f D T s for dif ferent SNR values. explained, as while the distribution of the estima ted f D is n ot symm etric, the estimated values are accum ulated aroun d their mean value. It should be mentioned that b y increasing the length of the observation wind ow , N , the bias o f the p roposed estimator approa c hes zero . In Fig. 5 , the NRMSE of th e N DA-MBE versus f d T s is illustrated fo r γ = 0 dB, γ = 10 dB, and γ = 20 dB. As seen , th e propo sed estimator exhibits a goo d perfo rmance over a wide range of Doppler rates, f D T s . As observed, the NRMSE decre ases as f D T s increases. This perform a n ce imp rovement can be expla in ed, as for lower Dopp ler rates, a larger observation window 19 1 2 4 6 8 10 12 14 16 18 10 -3 10 -3 10 -2 10 -1 10 0 Fig. 6. The ef fect of n t on the performance of the proposed ND A-MBE for n r = 2 and at γ = 10 dB. 1 2 4 6 8 10 12 14 16 18 10 -3 10 -3 10 -2 10 -1 10 0 Fig. 7. The ef fect of n r on the performance of the proposed ND A-MBE for n t = 2 and at γ = 10 dB. is re quired to capture the variation of the fading chann el. Also, as expected, the NRMSE decreases as γ inc reases. This can be easily explained, as an in crease in γ leads to more accura te estimates of the statistics in ( 3 8 ). Fig. 6 p resents the NRMSE of the p roposed NDA -MBE versus f D T s for different num bers of transmit an tennas, n t , fo r n r = 2 and at γ = 1 0 dB. As expected, the NRMSE increases as the num ber of transmit an te n nas increases. Th is inc r ease can be explained, as the variance of the statistics employed in ( 38 ) increases with the numb er of transmit antenna s, thus, lea ding to higher estimation error in the LS curve-fitting . In Fig. 7 , the NRMSE of the prop osed NDA-MBE is sh own versus f D T s for different nu m bers of rece ive antenn a s, n r , for n t = 2 , and at γ = 10 dB. I t can be seen that an incr ement in n r leads to a r e d uced NRMSE. T his decrease can be easily 20 1 2 3 4 5 6 7 8 9 10 -3 10 -2 10 -1 Fig. 8. The ef fect of the paramete r U min on the performance of the proposed NDA-MBE for differen t val ues of f d T s and at γ = 10 dB. explained, as averaging at the receive-side yields more ac c urate estima tio n of Ψ u , thus, leading to a mo re accu rate resu lt in the LS curve-fitting. In Fig. 8 , th e effect of the para meter U min on the perfo rmance of the pr oposed NDA -MBE is illu stra ted fo r L = 5 and u s = 1 0 . As ob served, the pr oposed estimato r exhibits a low sen siti vity to the value of U min . This can be explained , as a large number of lags, U min ≤ u ≤ U max , are e m ployed f or fitting e Ψ u to J 2 0 (2 π f D T s u ) in the LS estimatio n; thus, the estimator is nearly robust to a few m issing delay lags, L ≤ u < U min , or nuisan c e delay lags, 1 ≤ u < L . As such, b asically the estimator does not require an accurate estimate of L . Fig. 9 shows the effect of the ob servation window size, N , on the perform ance o f the pr o posed ND A-MBE. As expected, the perfo rmance of the prop osed e stima to r improves as the length of the observation window increases. Th is perfo rmance improvement can be explained, a s the variance of the estimated statistics employed in ( 38 ) decreases when N in creases. In Fig. 10 , the NRMSE is p lo tted versus f D T s for the pro posed NDA-MBE , the low-comp lexity DA-MLE ( D A-LMLE) in [ 33 ], [ 3 4 ], the ND A-CCE in [ 4 0 ], the D A-MLE in [ 30 ], and the D A-CRLB in [ 32 ] f o r MDS estimation in SI SO fr equency-flat fading chann el for N = 1 000 and at γ = 10 dB. As seen, the propo sed NDA-MBE o utperfo rms the NDA -CCE, and p rovides a similar perfo rmance as the D A-LMLE for f D T s ≥ 0 . 012 . The perform ance degradatio n of the DA-LMLE at hig h values of f D T s is related to the secon d -order T aylor expansion employed to ap proxima te the cov arian c e matrix; this is less a ccurate at higher MDSs. Fig. 11 illu stra te s the NRMSE versus f D T s for the prop osed NDA-MBE, th e derived DA-MLE, an d the derived D A-CRLB in MI MO fre quency-selective fadin g chann el for N = 1000 an d at γ = 10 dB. In ord er to show the co n vergence pr oblem in the deriv ed D A-ML E caused by th e Fisher-scoring numer ical m ethod employed to solve the ML [ 4 3 ], the p erforma nce of 21 1 2 4 6 8 10 12 14 16 18 10 -3 10 -3 10 -2 10 -1 10 0 10 1 10 2 Fig. 9. The ef fect of the observ ation windo w size, N , on the performanc e of the proposed ND A -MBE for n t = 2 and n r = 2 in frequency- selecti ve channel, and at γ = 10 dB. 1 2 4 6 8 10 12 14 16 18 10 -3 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Fig. 10. Performance comparison of the proposed NDA-MBE, D A-CRLB [ 32 ], D A-MLE [ 30 ], the low-c omplexity D A-MLE in [ 33 ], [ 34 ], and the NDA-C CE in [ 40 ] in SISO frequenc y-flat fading channel for N = 10 3 at γ = 10 dB. the d e riv ed DA -MLE for the cases of sin g le initial value (SIV ) and multiple initial values (MIV) is plotted , respec ti vely . 7 As seen, b y choosing MIV , the convergence p roblem of th e Fisher-scoring metho d employed in th e derived D A-MLE is solved. Moreover, as o bserved, the perfor mance of the derived DA-MLE with MIV is close to th e D A-CRLB. This h igh perfor mance is o btained at the expense of significant co m putational comp lexity in the o rder of O ( N 3 ) . On the other han d, the propo sed NDA-MBE canno t reach the D A-CRLB. T h is beh aviour can be explain ed, as th e NDA -MBE requ ir es a larger num ber of observation symbo ls to accurately estimate the secon d- and fou rth-ord er statistics in time-varying ch annel. Howev er, th e 7 W ith the MIV method, se veral initi al val ues are considered and at conv erge nce the one that yields the m aximum is chosen. 22 1 2 4 6 8 10 12 14 16 18 10 -3 10 -3 10 -2 10 -1 10 0 10 1 10 2 Fig. 11. Performance comparison of the proposed ND A-MBE, the deri ve d D A-CRLB, and the deri ved D A-MLE (SIV and MIV) in MIMO frequenc y-selecti ve Raylei gh fading channel for n t = 2 , n r = 2 , L = 5 , N = 10 3 , and at γ = 10 dB. substantial redu ced-com p lexity enables the p roposed MBE to exhibit significantly low NRMSE in the NDA scenar ios, where the observation wind ow can be selecte d large enou gh. In Fig . 12 , the NRMSE is shown versus f D T s for the p r oposed semi-blin d ND A-MBE in ( 5 2 ), the derived NDA -MLE, and the NDA-CRLB in SISO flat-fadin g chann el for BPSK signal, N = 10 , f D T s ∈ [5 × 10 − 3 , 4 5 × 10 − 3 ] , and at γ = 20 dB. 8 As expe c te d , th e ND A-MBE d oes no t exhibit go od per f ormance for a sho rt obser vation win dow size beca use the second- and fou rth-ord er statistics em p loyed in ( 52 ) are not accurately estimated. On the other hand, the derived NDA -MLE exhib its low NRMSE even for a shor t observation win dow . Mor eover , there is no sig n ificant perfo rmance g ap between the derived ND A-MLE and ND A-CRLB, as well as the DA-MLE in [ 30 ] and the D A-CRLB in [ 32 ]. V I I I . C O N C L U S I O N In this pa p er , we derived the DA - and ND A-CRLBs an d DA- and NDA-MLE s for MDS in MIM O fre quency-selective fadin g channel. Mor e over , a low-complexity ND A-MBE fo r MISO and MIMO systems was pr oposed. The NDA-MBE emp loys the statistical momen t-based appro ach and relies o n the secon d- a nd f ourth-o rder statistics of the received signal, as well as the LS curve-fitting optim ization techn ique. Compared to the existing DA estimators, the pro posed ND A-MBE provide s h igher system capacity due to absence of p ilot. Also, the substantial reduced-c o mplexity enab les the p roposed MBE to exhibit good perfor mance in th e NDA scenarios, where the observation window can b e selected large enoug h. The NDA-MBE does not require a priori knowledge of other pa rameters, such as the nu m ber of transmit antennas; furth ermore, the propo sed ND A- MBE is robust to the time-freq uency asynch ronization . When co mpared to th e NDA-CC E, the ND A-MBE exhibits better 8 The complexity order of the deri ved ND A-MLE and NDA-CRLB are in the order of O ( | M | N ′ n t ) ; for large value s of N ′ ( N ′ = N + L − 1 ), the correspond ing curve s are not obtainabl e ev en for | M | = 2 or n t = 1 . Hence, N = 10 and SISO flat-fadi ng channel are considere d. 23 0.005 0.015 0.025 0.035 0.045 10 -1 10 0 10 1 Fig. 12. Performance comparison of the proposed semi-blind ND A-MBE in ( 52 ), the deri ved NDA-MLE, the deri ved ND A-CRLB, the DA-MLE in [ 30 ], and the DA-CRLB in [ 32 ] in SISO frequenc y-flat-fa ding channe l for N = 10 and at γ = 20 dB SNR.   r ( n ) k   2 = n t X m =1 L X l =1   h ( mn ) k,l   2   s ( m ) k − l   2 + n t X m 1 =1 n t X m 2 6 = m 1 L X l =1 h ( m 1 n ) k,l  h ( m 2 n ) k,l  ∗ s ( m 1 ) k − l  s ( m 2 ) k − l  ∗ (53) + n t X m =1 L X l 1 =1 L X l 2 6 = l 1 h ( mn ) k,l 1  h ( mn ) k,l 2  ∗ s ( m ) k − l 1  s ( m ) k − l 2  ∗ + n t X m 1 =1 n t X m 2 6 = m 1 L X l 1 =1 L X l 2 6 = l 1 h ( m 1 n ) k,l 1  h ( m 2 n ) k,l 2  ∗ s ( m 1 ) k − l 1  s ( m 2 ) k − l 2  ∗ +  w ( n ) k  ∗ n t X m =1 L X l =1 h ( mn ) k,l s ( m ) k − l +  w ( n ) k  n t X m =1 L X l =1  h ( mn ) k,l  ∗  s ( m ) k − l  ∗ +   w ( n ) k   2 . κ ( n ) u = n t X m =1 L X l =1 E n   h ( mn ) k,l   2 | h ( mn ) k + u,l   2 o E n   s ( m ) k − l   2   s ( m ) k + u − l   2 o + n t X m 1 =1 n t X m 2 6 = m 1 L X l =1 E n   h ( m 1 n ) k,l   2 | h ( m 2 n ) k + u,l   2 o E n   s ( m 1 ) k − l   2   s ( m 2 ) k + u − l   2 o + n t X m =1 L X l 1 =1 L X l 2 6 = l 1 E n   h ( mn ) k,l 1   2   h ( mn ) k + u,l 2   2 o E n   s ( m ) k − l 1   2   s ( m ) k + u − l 2   2 o + n t X m 1 =1 n t X m 2 6 = m 1 L X l 1 =1 L X l 2 6 = l 1 E n   h ( m 1 n ) k,l 1   2   h ( m 2 n ) k + u,l 2   2 o E n   s ( m 1 ) k − l 1   2   s ( m 2 ) k + u − l 2   2 o + E n   w ( n ) k + u   2 o n t X m =1 L X l =1 E n   h ( mn ) k,l   2 o E n   s ( m ) k − l   2 o + E n   w ( n ) k   2 o n t X m =1 L X l =1 E n   h ( mn ) k + u,l   2 o E n   s ( m ) k + u − l   2 o + E n   w ( n ) k   2 o E n   w ( n ) k + u   2 o . (54) perfor mance, and when comp ared to th e low-complexity DA-MLE, it exhibits similar perf o rmance for hig h MDSs. On the other han d, the d erived D A-MLE’ s perfo rmance is very close to the deriv ed D A-CRLB in MIMO frequen cy-selecti ve channel ev en when the observation window is relati vely small. Similarly , there is n o significant perform ance g ap between the derived ND A-MLE and the ND A-CRLB. 24 A P P E N D I X I T o obtain an explicit closed-fo rm expression for κ ( n ) u ∆ = E  | r ( n ) k | 2 | r ( n ) k + u | 2  , we first write   r ( n ) k   2 = r ( n ) k  r ( n ) k  ∗ by emp loyin g ( 1 ), as in ( 53 ) at the top of next p a g e. Th en,   r ( n ) k + u   2 is straightforwardly calcu lated by r eplacing k with k + u in ( 53 ), and   r ( n ) k   2   r ( n ) k + u   2 can be easily expressed in a summation form, which is omitted due to space constra in ts. As fading is indepen dent o f the signal an d noise, the statistical expectatio n in κ ( n ) u can be de composed into statistical expec- tations over the signal, fading, and no ise distributions, r espectively . With in depend e n t and iden tically distributed tra n smitted symbols, s ( m ) k , k = 1 , . . . , N , and u ≥ L , the sym bols fr om the m th an tenna c ontributed in r ( n ) k , i.e., { s ( m ) k − l } L l =1 , are d ifferent from those con tributed in r ( n ) k + u , i.e., { s ( m ) k + u − l } L l =1 ; furth ermore, by using th at E  ( s ( m ) k ) 2  = 0 , E  s ( m ) k  = 0 , and the linearity proper ty of the statistical expectation, one obtain s κ ( n ) u as in ( 54 ). Finally , with E  | s ( m ) k | 2  = σ 2 s m and E  | w ( n ) k | 2  = σ 2 w n , ( 28 ) is obtained .  A P P E N D I X I I By employing ( 33 ), one can write  µ ( n ) 2  2 = n t X m =1 L X l =1 σ 2 h ( mn ) ,l σ 2 s m + σ 2 w n ! 2 = n t X m =1 L X l =1 σ 4 h ( mn ) ,l σ 4 s m + n t X m 1 =1 n t X m 2 6 = m 1 L X l =1 σ 2 h ( m 1 n ) ,l σ 2 h ( m 2 n ) ,l σ 2 s m 1 σ 2 s m 2 + n t X m =1 L X l 1 =1 L X l 2 6 = l 1 σ 2 h ( mn ) ,l 1 σ 2 h ( mn ) ,l 2 σ 4 s m + n t X m 1 =1 n t X m 2 6 = m 1 L X l 1 =1 L X l 2 6 = l 1 σ 2 h ( m 1 n ) ,l 1 σ 2 h ( m 2 n ) ,l 2 σ 2 s m 1 σ 2 s m 2 + 2 σ 2 w n n t X m =1 L X l =1 σ 2 h ( mn ) ,l σ 2 s m + σ 4 w n . (55) Then, by subtractin g  µ ( n ) 2  2 in ( 55 ) from κ ( n ) u in ( 32 ), one obtain s κ ( n ) u −  µ ( n ) 2  2 = J 2 0 (2 π f D T s u ) η ( n ) , (56) where η ( n ) = 1 / n t P m =1 L P l =1 σ 4 h ( mn ) ,l σ 4 s m . A P P E N D I X I I I W ith indepen d ent fadin g , no ise, and signal proce sses, by using ( 53 ) and th en ( 3 1 ) and that E { ( s ( m ) k ) 2 } = E { s ( m ) k } = E  h ( mn ) k,l  2  = E  w ( n ) k  = 0 , similar to Appendix I , one obtains µ ( n ) 4 = E n | r ( n ) k | 4 o = 2 n t X m =1 L X l =1 σ 4 h ( mn ) ,l σ 4 s m (57) + 2 n t X m =1 L X l =1 σ 4 h ( mn ) ,l (Ω s − 1) σ 4 s m 25 + 2 n t X m 1 =1 n t X m 2 6 = m 1 L X l =1 σ 2 h ( m 1 n ) ,l σ 2 h ( m 2 n ) ,l σ 2 s m 1 σ 2 s m 2 + 2 n t X m =1 L X l 1 =1 L X l 2 6 = l 1 σ 2 h ( mn ) ,l 1 σ 2 h ( mn ) ,l 2 σ 4 s m + 2 n t X m 1 =1 n t X m 2 6 = m 1 L X l 1 =1 L X l 2 6 = l 1 σ 2 h ( m 1 n ) ,l 1 σ 2 h ( m 2 n ) ,l 2 σ 2 s m 1 σ 2 s m 2 + 4 σ 2 w n n t X m =1 L X l =1 σ 2 h ( mn ) ,l σ 2 s m + 2 σ 4 w n , where Ω s is th e fourth- order two-conju gate statistic for un it variance sign al, which represen ts the effect of the modu lation format. Finally , by emp loying ( 33 ) and ( 57 ), ( 35 ) is easily obtained .  R E F E R E N C E S [1] A. F . Molisch, W ir eless Communicati ons . John W ile y & Sons, 2007. [2] G. L. St ¨ uber , Principles of Mobile Communication . Springer Science & Business Media, 2011. [3] J. Wu and P . Fan, “ A survey on high mobility wireless communicati ons: Challenges, opportunit ies and solutions, ” IEEE Access , vol. 4, pp. 450–476, Jan. 2016. [4] M. J. Chu and W . E. Stark, “Effect of mobile velocit y on communications in fading channels, ” IEEE T rans. V eh. T echnol. , vol. 49, no. 1, pp. 202–210, Jan. 2000. [5] A. J. Goldsmith and S.-G. Chua, “V ariab le-rate varia ble-po wer M-QAM for fading channel s, ” IEEE T rans. Commun. , vol. 45, no. 10, pp. 1218–1230, Oct. 1997. [6] K. Balacha ndran, S. R. Kadaba, and S. Nanda, “Channel quality estimation and rate adapt ation for cellular m obile radio, ” IEEE J . Sel. Area s Commun. , vol. 17, no. 7, pp. 1244–1256, J uly 1999. [7] C. T ep edelenlio ˘ glu, A. Abdi, G. B. Giannakis, and M. Kave h, “Estimation of Doppler spread and signal strength in mobile communications with applic ations to handof f and adapti ve transmission, ” W ir el. Commun. Mob . Comput , vol. 1, no. 2, pp. 221–242, Mar . 2001. [8] M. D . Austin and G. L. St ¨ uber , “V elocity adapti ve handoff algorithms for microcel lular systems, ” IEEE T rans. V eh. T echnol. , vol. 43, no. 3, pp. 549–561, Aug. 1994. [9] B. L. Mark and A. E. Leu, “Local avera ging for fast handof fs in cellular networks, ” IEEE T rans. W ire less Commun. , vol. 6, no. 3, pp. 866–874, Mar . 2007. [10] A. Duel-Hal len, S. Hu, and H. Hallen, “Long-range predicti on of fading signals, ” IEE E Signal Pr ocess. Mag. , vol. 17, no. 3, pp. 62–75, May 2000. [11] D.-S. Shiu, G. J . Foschini, M. J. Gans, and J. M. Kahn, “Fading correlatio n and its effec t on the capaci ty of multi-el ement antenna systems, ” IEEE T rans. Commun. , vol. 48, no. 3, pp. 502–513, Mar . 2000. [12] C.-N. Chuah et al. , “Capaci ty scali ng in MIMO wireless systems under correlate d fadi ng, ” IEEE T rans. Inf. Theory , vol. 48, no. 3, pp. 637–650, Mar . 2002. [13] T . Y ucek, R. T annious, and H. Arslan, “Doppl er spread estimation for wireless OFDM systems, ” in Proc. IEEE Advances in W ir ed and W ire less Communicat ion , Princeton , USA, Apr . 2005, pp. 233–236. [14] S. Coleri, M. Ergen, A. Puri, and A. Bahai, “Channe l estimat ion tec hniques based on pilot arrangement in OFDM systems, ” IEEE T rans. Br oadcast. , vol. 48, no. 3, pp. 223–229, Sept. 2002. [15] F . Bellili, R. Meftehi, S. Affes, and A. St ´ ephenne, “Maximum like lihood SNR estimation of linearl y-modulated signals ov er time-vary ing flat-fading SIMO channels, ” IE EE T rans. Signal Proc ess. , vol. 63, no. 2, pp. 441–456, Jan. 2015. [16] M. Marey , M. Samir , and O. A. Dobre, “EM-based joint channe l estimation and IQ imbalances for OFDM systems, ” IEE E Tr ans. Broad cast. , vol. 58, no. 1, pp. 106–113, Mar . 2012. [17] A. Khansefid and H. Minn, “On channel estimation for massiv e MIMO with pilot contaminati on, ” IEEE Commun. Lett , vol. 19, no. 9, pp. 1660–1663, Sept. 2015. 26 [18] Y . Li, H. Minn, and M. Z. W in, “Frequ ency offse t estimation for MB-OFDM-based UW B systems, ” IEE E T rans. Commun. , vol. 56, no. 6, pp. 968–979, Jun. 2008. [19] M. Mohammadkarimi, O. A. Dobre, and M. Z. W in, “Non-Data-Aided SNR Estimation for Multiple Antenna S ystems, ” in Proc. GLOBECOM , W ashi ngton, USA, Dec. 2016, pp. 1–5. [20] H. W ang, O. A. Dobre, C. Li, and D . C. Popescu, “Blind cyc lostationa rity-based symbol period estimation for FSK signals, ” IEEE Commun. Lett , vol. 19, no. 7, pp. 1149–1152, J uly 2015. [21] A. Masmoudi, F . Belli li, S. Aff es, and A. Stephenne, “ A non-data-a ided maximum likeliho od time delay estimator using importance s ampling, ” IEE E T rans. Signal Pr ocess. , vol . 59, no. 10, pp. 4505–4515, Oct. 2011. [22] M. Mohammadkarimi, E. Karami, O. A. Dobre, and M. Z. W in, “Number of transmit antenna s detect ion using time-di versi ty of the fading channel, ” IEEE T rans. Signal Pr ocess. , vol. 65, no. 15, pp. 4031–4046, Aug. 2017. [23] A. St ´ ephenne, F . Belli li, and S. Affes, “Moment-based SNR estimation ove r linearl y-modulated wirele ss SIMO channel s, ” IE EE T rans. W ireless Commun. , vol. 9, no. 2, pp. 714—-722, Feb . 2010. [24] F . Belli li, A. Methenni, S. B. Am or , S. Affes, and A. St ´ ephenne, “Time synchroniza tion of turbo-coded square-QAM-modulate d transmissions: code-aide d ML estimator and closed-form Cram ´ er-Ra o lower bounds, ” IEEE T rans. V eh. T echnol. , vol. PP , no. 91, pp. 1–1, 2017. [25] F . Bellili, A. Methen ni, and S. Affes, “Closed-form CRLBs for SNR estimat ion from turbo-coded BPSK-, MSK-, and square-QAM-modula ted signals, ” IEEE T rans. Signal Pr ocess. , vol. 62, no. 15, pp. 4018–4033, Aug. 2014. [26] N. W u, H. W ang, J. Kuang, and C. Y an, “Performance anal ysis of code-aided sym bol timing reco very on A WGN channels, ” IEE E T rans. Commun. , vol. 59, no. 7, pp. 1975–1984, Mar . 2011. [27] C. Herzet, N. Noels, V . Lottici, H. W ymeersch, M. Luise, M. Moenecla ey , and L. V andendorpe , “Code-aide d turbo synchroniz ation, ” Proc. IEEE , vol. 95, no. 6, pp. 1255–1271, July 2007. [28] F . Simoens and M. Moenecla ey , “Reduced complexi ty data-ai ded and code-aided frequenc y of fset estimation for flat-fadi ng MIMO channel s, ” vol. 5, no. 6, pp. 1558–1567, Jun 2006. [29] J. Sun and M. C. V alenti, “Joint synchronizati on and SNR estimation for turbo codes in A WGN channe ls, ” IEEE T rans. Commun. , vol. 53, no. 7, pp. 1136–1144, J uly 2005. [30] L. Krasny , H. Arslan, D. Koil pillai, and S. Chenna keshu, “Doppler spread estimation in mobile radio systems, ” IEE E Commun. Lett. , vol. 5, no. 5, pp. 197–199, May 2001. [31] Y .-R. Tsai and K.-J. Y ang, “ Approximate ML Doppler spread estimatio n ov er flat Rayle igh fading channels, ” IE EE Signal Proc ess. Lett. , vol. 16, no. 11, pp. 1007–1010, Nov . 2009. [32] A. Dogand ˇ zi ´ c and B. Zhang, “Estimating Jakes’ Doppler powe r s pectrum paramet ers using the Whittle approximat ion, ” IEEE T rans. Signal Proc ess. , vol. 53, no. 3, pp. 987–1005, Mar . 2005. [33] F . Bell ili, Y . Selmi, S. Affes, and A. Ghrayeb, “ A lo w-cost and robust maximum lik elihood joint estimator for the Doppler spread and CFO parameters ov er flat-fa ding Rayleigh channel s, ” IEE E T rans. Commun. , pp. 1–1, 2017. [34] F . Bellil i and S. Affes, “ A low-cost and robust maximum likeliho od Doppler spread estimator , ” in Pr oc. GLOBECOM , Atlanta, USA, Dec. 2013, pp. 4325–4330. [35] K. E. Baddour and N. C. Beauli eu, “Nonparametri c Doppler spread estimation for flat fading channels, ” in Proc. WCNC , New Orleans, USA, 2003, pp. 953–958. [36] K. D. Anim-Appiah, “On genera lized cov ariance -based velocit y estimation, ” IEEE T rans. V eh. T echnol . , vol. 48, no. 5, pp. 1546–1557, Sept. 1999. [37] A. Abdi, H. Zhang, and C. T epe delenliogl u, “ A unified approac h to the performance analysis of speed estimation techniques in mobile communication, ” IEEE T rans. W ir eless Commun. , vol . 56, no. 1, pp. 126–135, Jan. 2008. [38] M. Souden, S. Af fes, J. Benest y , and R. Bahroun, “Robust Doppler spread estimat ion in the presence of a residual carrier frequen cy offset , ” IEEE T rans. Signal P r ocess. , vol. 57, no. 10, pp. 4148–4153, Oct. 2009. [39] G. Park, D. Hong, and C. Kang, “Lev el crossing rate estimation with Doppler adapt iv e noise suppression techniqu e in frequency domain, ” in Proc. IEEE VTC , Orlando, USA, Oct. 2003, pp. 1192–1195. [40] H. Zhang and A. Abdi, “Cycl ostationari ty-based Doppler spread estima tion in mobile fa ding channel s, ” IEEE T ra ns. Commun. , vol. 57, no. 4, pp. 1061–1067, A pr . 2009. 27 [41] E. Biglieri, R. Calderbank, A. Constantinid es, A. Goldsmith, A. Paulraj, and H. V . Poor , MIMO W ire less Communicat ions . Cambridge Uni versity press, 2007. [42] N. T . L ongford, “ A fast s coring algorit hm for maximum lik elihood estimat ion in unbalanced m ixed models with nested random ef fects, ” Biometrika , vol. 74, no. 4, pp. 817–827, 1987. [43] S. M. Kay , Fundamental s of Statistical Signal Proc essing, V ol. I: Estimation Theory . Prent ice Hall, 1993. [44] A. Swami and B. M. Sadler , “Hierarch ical digital m odulat ion classificatio n using cumulant s, ” IEEE T rans. Commun. , vol. 48, no. 3, pp. 416–429, Mar . 2000. [45] K. E. Baddour and N. C. Beauli eu, “ Autoregre ssiv e modeli ng for fading chann el simulation, ” IEE E T ran s. W ir eless Commun. , vol. 4, no. 4, pp. 1650–1662, July 2005. [46] P . Sadeghi , R. A. Kenne dy , P . B. Rapajic, and R. Shams, “Finite -state Markov modeling of fading channels-a surve y of principle s and appl ications, ” IEEE Signal Proc ess. Mag . , vol. 25, no. 5, pp. 57–80, Sept. 2008. [47] I. Reed, “On a moment theorem for complex Gaussian processes, ” IEEE T rans. Inf. Theory , vol. 3, no. 8, pp. 194–195, Apr . 1962. [48] A. R. Gallant, Nonlinear Statist ical Models . John Wi ley & Sons, 2009, vol. 310. [49] N. R. Draper , H. Smith, and E . Pownell, A pplied Re gressio n Analysis . W ile y New Y ork, 1966, vol. 3. [50] S. S. Rao and S. Rao, Engineering Optimizati on: Theory and Practi ce . John W iley & Sons, 2009. [51] A. M. Zoubir and B. Boashash, “The bootstrap and its applica tion in signal proce ssing, ” IE EE Signal P r ocess. Mag. , vol. 15, no. 1, pp. 56–76, Jan. 1998. [52] A. M. Zoubir and D. R. Iskandler , “Bootstrap methods and application s, ” IEE E Signal Pr ocess. Mag . , vol. 24, no. 4, pp. 10–19, July 2007. [53] A. M. Zoubir and D. R. Iskander , “Bootstrap modeling of a class of nonstation ary signal s, ” IE EE T rans. Signal Pr ocess. , vol. 48, no. 2, pp. 399–408, Feb . 2000.