6 Deterministic channel parameter estimation inthishaprer we inroduce theoand reiques foretmhanelpreters rom mesuem (.thecap-fmethd()d the MUSIC Schmidt (c the maximum likelihood estimation method is elaborated we introduce the EM and the SAGE is described.In Section6.5,some new techniques relying on the Bayesian estimation methods are introduced.In section subspace-bas ebased techniquesare desribed which are applicable for tracking the variation of channels with relatively ow complexity. 6.1 Spectral-based method either maximiznor minimin the spectral height over the parameter(s),it is possible to obtain certain estimates o cteristi of the componentsin the propagation ett (1948).the Ca beamformer Capon (199),and the MUSIC Schmidt (16),as well as numerous variants of these method ic methods and the n-Tukey In the sequel,we briefly introduce these techniques,and illustrate the results when they are applied to processing the measurement data collected in real environment. 6.1.1 Beamforming method 了p In order tounderstand what is the refinements of the Bartlett beamformer,it is necessary to know the periodogram d new the the equation l=它 (6.1) thathas very regular formulation,which does not include the intrinsic system responses for thethsample. O础y
6 Deterministic channel parameter estimation In this chapter, we introduce the algorithms and techniques used for estimating channel parameters from measurement data. In Section 6.1, we introduce the spectral-based method, such as the Bartlett beamforming method Bartlett (1948), the Capon-beamforming method Capon (1969), and the MUSIC algorithm Schmidt (1986). In Section 6.2, the maximum likelihood estimation method is elaborated. In Section 6.3, we introduce the EM and the SAGE algorithm used for channel parameter estimation. In Section 6.4, the Richter’s Maximum likelihood estimation method is described. In Section 6.5, some new techniques relying on the Bayesian estimation methods are introduced. In Section 6.6, typical subspace-based method is elaborated. In Section 6.7, the Kalman filtering algorithms are described which are used to track the variation of the parameters of the channels. Finally, in Section 6.8, the modified particle filtering based techniques are described, which are applicable for tracking the variation of channels with relatively low complexity. 6.1 Spectral-based method The methods belonging to the category of spectral-based method have a common feature that a smooth spectrum with respect to a specific parameter, such as delay, direction of arrival, etc., is computed based on the observations. By either maximizing or minimizing the spectral height over the parameter(s), it is possible to obtain certain estimates of the parameters that represent the characteristics of the components in the propagation channel. The spectral-based method includes the periodogram, correlogram, Bartlett beamformer Bartlett (1948), the Capon beamformer Capon (1969), and the MUSIC Schmidt (1986), as well as numerous variants of these methods. Generally speaking, the spectral-based methods can be grouped into two classes, i.e. the nonparametric methods and the parametric methods. Typical nonparametric methods are the periodogram, correlogram, the Blackman-Tukey method, and the refined Blackman-Tukey method based on windows, and other refined periodogram methods, such as the Bartlett method, the Welch method, Daniell method, etc. To be added some new methods. Include some new citations after 2005 for examples. In the sequel, we briefly introduce these techniques, and illustrate the results when they are applied to processing the measurement data collected in real environment. 6.1.1 Beamforming method The reason that we focus on the refined periodogram method is that for channel parameter estimation, we cope with the samples collected not only in time, in frequency, but also in space. Sometimes, we also need to consider polarization domains. In order to understand what is the refinements of the Bartlett beamformer, it is necessary to know the periodogram and correlogram spectral estimators first. The periodogram power spectral estimation method can be represented by the computations in the following equation: pˆp(ω) = 1 N X N n=1 y(n)e −jωn 2 , (6.1) where y(n) is the received signal observation for the nth sampling instance in time or location in space, e −jωn represents the system response at the nth sample when the signal component has the frequency of ω. Its obvious that e −jωn has very regular formulation, which does not include the intrinsic system responses for the nth sample. This is a Book Title Name of the Author/Editor c XXXX John Wiley & Sons, Ltd
152 Deterministic channel parameter estimation The correlogram spectral estimation method is c(o)- (6.2) k=(-N-1) where hatr()represents the autocorrelation function of u(n). 1.For a given direction of arrival (DoA)0,the filter passed undistorted signals. 2.The filter attenuates all the other DoAs different from. hould the cmdu tha)h ond to the c s”of the spatial filte c8) h=c(0)We(0) (6.3) Thus,when this optimal weights are applied to the spatial filter,we obtain for the output signal of the spatial filter: Eflu(t))=h"Rh =hE(y(t)y(t)"] _c(0)Eiy(t)(t))c(0) (6.4) c()Hc(6) This equation shows that when the DoA of the impinging signal isonly when the weights of the fiter satisfy (6.3). c(0)Rc(e) pB()= (6.5) c(0Hc(0) whereis the covariance matrix of the received signal,which can be calculated according to R=E(y(t)u(t)) (6.6) nshebembgmhmhehenhngerecetedsgalsnom likelihood e thod Inmthe following an example of using the boinmethod the dband MIMO underOSou-builin ofu were not used in the Rx during the measurements.Thus,totally50 1000 spatial subchannels are measured in one measurement cycle. polarization and horizontal polarization.In this example,both vertical and horiz nt al polarization array responses are ctra ct on the
152 Deterministic channel parameter estimation The correlogram spectral estimation method is pˆc(ω) = N X−1 k=(−N−1) rˆ(k)e −jωk (6.2) where hatr(k) represents the autocorrelation function of y(n). Beamforming method is based on the assumption that the array response, which is also called steering vector c(θ) is know. The beamforming method is proposed for designing a spatial filter, which satisfies the following two conditions: 1. For a given direction of arrival (DoA) θ, the filter passed undistorted signals. 2. The filter attenuates all the other DoAs different from θ. These two conditions correspond to the constraint that the weights h of the multiple spatial “taps” of the spatial filter should satisfy the condition that min h h Hh subject to h Hc(θ) = 1 ??. It can be shown that in such a case, h = c(θ) c(θ)Hc(θ) . (6.3) Thus, when this optimal weights are applied to the spatial filter, we obtain for the output signal of the spatial filter: E{|y(t)| 2 } = h HRh = h HE{y(t)y(t) H} = c(θ) HE{y(t)y(t) H}c(θ) c(θ)Hc(θ) (6.4) This equation shows that when the DoA of the impinging signal is θ, only when the weights of the filter satisfy (6.3), the output signal has the maximum power. Thus, when the DoA of the impinging signal is unknown, we can vary the weights of the filter by letting θ being selected within a certain range. Therefore, a pseudo-power spectrum p(θ) with respect to θ can be computed. By selecting the specific value of θ which leads to the peak of p(θ), the estimate of the DoA is obtained. Based on this rationale, the power spectrum computed can be expressed as pB(θ) = c(θ) HRcb (θ) c(θ)Hc(θ) , (6.5) where Rb is the covariance matrix of the received signal, which can be calculated according to Rb = E{y(t)y(t) H} = 1 N X N n=1 yn(t)yn(t) H (6.6) discussing about the rank... It is worth mentioning that the maximum likelihood estimation method when assuming the received signal is normal distributed has the same expression as the beamforming method in the single-path scenario. In the following an example of using the beamforming method to estimate the direction of arrival is illustrated. In this example, measurement data is used which was collected in a campaign jointly conducted by Elektrobit and Technology University of Vienna in 2005 using the wideband MIMO sounder -PROPSound- in a building of Oulu University. Figure 6.1 (a) depicts the photograph of the 50-element antenna array used in the Tx and Rx during the measurements. Figure 6.1 (b) illustrates the indices of the antennas in the array. Notice that 18 antennas from No. 19 to 36 were not used in the Rx during the measurements. Thus, totally 50 × 32 = 1600 spatial subchannels are measured in one measurement cycle. Fig. 6.3 depicts the power spectrum estimated by using the beamforming method. The signals emitted by 50 Tx antennas and received from the first Rx antenna are used. The responses of the antennas are separated into vertical polarization and horizontal polarization. In this example, both vertical and horizontal polarization array responses are considered. It can be observed from Fig. 6.3 that the power spectra are not exactly the same. The spectrum estimated by using the vertical polarization has more fluctuations than that observed when horizontal polarization is considered. Furthermore, it can be seen that the power spectrum estimated by using the horizontal polarization has higher spectral height than that computed by using the vertical polarization. It is obvious that the array response may have significant impact on the estimated power spectrum
Deterministic channel parameter estimation 153 (a)50-element array (b)Antenna indices in the array 网 Figure 6.1 (a)The photograph of the antenna array used in measurements and (b)the indices of the array element (a)Vertical polarization (b)Horizontal polarizatio pulNNIaNTpnes 10 -6Azimuth 100 heiiof depecalaed byteebynand )tHortontalpobrizatios 30 3 rture calculated by using the beamforming method by and
Deterministic channel parameter estimation 153 (a) 50-element array (b) Antenna indices in the array 49 50 / 37 38 / 39 40 / 41 42 / 43 44 / 45 46 / 47 48 / 1 2/ 3 4/ 5 6/ 7 8/ /9 10 11 12 / 13 14 / 15 16 / 17 18 / 21 22 / 23 24 / 25 26 / / 27 28 29 30 / 31 32 / 33 34 / 35/36 19 20 / Figure 6.1 (a) The photograph of the antenna array used in measurements and (b) the indices of the array element −150 −100 −50 0 50 100 150 −80 −60 −40 −20 0 20 40 60 80 2 4 6 8 10 12 14 16 x 10−6 Azimuth [ ◦ ] Elevation [ ◦ ] Spectral height [linear] (a) Vertical polarization −150 −100 −50 0 50 100 150 −80 −60 −40 −20 0 20 40 60 80 1 2 3 4 5 6 7 8 9 10 x 10−5 Azimuth [ ◦ ] Elevation [ ◦ ] Spectral height [linear] (b) Horizontal polarization Figure 6.2 Power spectrum in direction of departure calculated by using the beamforming method by using 50 Tx antennas and the no. 1 Rx antenna −150 −100 −50 0 50 100 150 −80 −60 −40 −20 0 20 40 60 80 0.5 1 1.5 2 2.5 3 x 10−5 Azimuth [ ◦ ] Elevation [ ◦ ] Spectral height [linear] (a) Vertical polarization −150 −100 −50 0 50 100 150 −80 −60 −40 −20 0 20 40 60 80 2 4 6 8 10 12 14 16 x 10−5 Azimuth [ ◦ ] Elevation [ ◦ ] Spectral height [linear] (b) Horizontal polarization Figure 6.3 Power spectrum in direction of departure calculated by using the beamforming method by using 50 Tx antennas and the no. 2 Rx antenna
154 Deterministic channel parameter estimation 6.2 Maximum-likelihood method The maximum-likelihood estimation is considered as a standard method for estimating channel parameters based or certain predefined models.In this subsec ion,this method is briefly introduced by using the channel mode where the A narro (6.7刀 with spectral height of N The DoA is a unit vector uniquely determined by azimuth of arrival (AoA) and elevation of arrival 6 0,as =[sin()cos(),sin(@)sin(),cos()]. (6.8) 日=(a1.21,a2.22...a.2w (6.9) vtheohe ators of the model pa arameters can be calculated by maximizing the log-likelihood function of the A(e)=logp[eh. (6.10) npyin thethat the obe twhe The estimators can be obtained as 合L=argm8xA(⊙) (6.11) oblem to be solved in(6.11)is practically prohibitive due to the complexity of exhaustive searching mum-likelihood estimator can be evaluated by using the Cramer-Rao bound.In Fleury eta(1999).the derivation of the Cramer-Rao bound is provided.The readers are referred to this paper for the details 6.3 The EM and SAGE algorithms arrival of indi propagat paths.In Fleury et al.(2002a,b)the SAGE algorithm was e tended to include the initialization-and mproved SAGE (SAGB)algorithm.The acronym ISI stresses the fact that the ntiaizti effect of polariz
154 Deterministic channel parameter estimation 6.2 Maximum-likelihood method The maximum-likelihood estimation is considered as a standard method for estimating channel parameters based on certain predefined models. In this subsection, this method is briefly introduced by using the channel model where the propagation paths are only characterized by their directions of arrival (DoAs). Notice that this simplified model can be replaced by other more complicated channel models when deriving the ML estimators of the parameters of the model. A narrowband 1 × N single-input multiple-output (SIMO) scenario is considered where the propagation paths between the Tx and Rx have different DoA Ω. Following the nomenclature in Yin et al. (2003a), the narrowband representation of the impulse responses h ∈ CN of the SIMO channel can be written as h = X L ℓ=1 αℓc(Ωℓ) + w, (6.7) where ℓ is the index of specular propagation paths, L represents the number of paths, αℓ and Ωℓ are respectively, the complex attenuation and the direction of arrival (DoA) of the ℓth path, w represents the standard white Gaussian noise with spectral height of No. The DoA Ωℓ is a unit vector uniquely determined by azimuth of arrival (AoA) φ ∈ [−π, π] and elevation of arrival θ ∈ [0, π] as Ω = [sin(θ) cos(φ),sin(θ) sin(φ), cos(θ)]. (6.8) In (6.75), c(Ω) represents the array response at a given DoA. The parameters of interest for estimation in the above mentioned generic channel model (6.75) are Θ = (α1, Ω1, α2, Ω2, . . . , αN , ΩN ). (6.9) The ML estimators Θb ML of the model parameters can be calculated by maximizing the log-likelihood function of the Θ given the observations of h. Λ(Θ) = log p[Θ|h, (6.10) where p(θ) denotes the likelihood function of Θ. Applying the assumption that the noise component w is white Gaussian, it can be shown that Λ(Θ) = −N log(2πσw) − 1 2Nσ2 w xˆ [i+1] ℓ − X L ℓ=1 αℓc(Ωℓ) 2 . The estimators Θb ML can be obtained as Θb ML = arg max Θ Λ(Θ). (6.11) The optimization problem to be solved in (6.11) is practically prohibitive due to the complexity of exhaustive searching in multiple dimensions. The performance of the maximum-likelihood estimator can be evaluated by using the Cramer-Rao bound. In Fleury et al. (1999), the derivation of the Cramer-Rao bound is provided. The readers are referred to this paper for the details. 6.3 The EM and SAGE algorithms The space-alternating generalized expectation-maximization (SAGE) algorithm was applied for channel parameter estimation in Fleury et al. (1999) for estimating the delay, Doppler frequency, direction (azimuth and elevation) of arrival of individual propagation paths. In Fleury et al. (2002a,b) the SAGE algorithm was extended to include the estimation of the directions of departure of paths and used to estimate path parameters with the measurement data collected with multi-antenna array installed in both the transmitter and the receiver. This SAGE algorithm is called initialization-and-search-improved SAGE (ISI-SAGE) algorithm. The acronym ISI stresses the fact that the initialization and search procedures of the SAGE algorithm are optimized to speed up its convergence and enhance its capability of detecting weak paths. In order to accurately model the dispersion characteristics of the channels that incorporate the effect of polarizations, the ISI-SAGE algorithm has been extended in Fleury et al. (2003); Yin et al. (2003b) to include estimation of the polarization matrices of waves propagating from the Tx site to the Rx site in a MIMO system
Deterministic channel parameter estimation 155 ☒ exp(j2u,） -) A ) Figure 6.4 Contribution of the (th wave to the received signal in a MIMO system incorporating dual-antenna arrays In this chapter,we describe the signal model which uses 14 parameters to characterize a single specular propagation a设 ions to s ecial cases,an analysis addressir n specific scenarios such as complexity. 6.3.1 Signal model As shown in Figure 6.4.Rx and Tx antenna arrays have dual antenna configuration.Each of the dual antenna ransmits/rec on the potarbai s in the underlying pattern.The other ced ofth weighs for the propagation paths.The the of to the output of the MIMO system reads st)=ep2wd[e）c2282[a）c2ut-小. (6.12 从二兰 s(t:0)=exp(j2wvet)C2(S2)AC(S.i)u(t-T), (6.13)
Deterministic channel parameter estimation 155 Polarization Matrix of Path " " A (ȍ ) ,11,1,1 " c ( ) 1 " tu W ( ) 1 " u t W M ( ) " u t W (ȍ ) 1,1,2 ,2 " c ),( 1 ș" ts );( ș" exp( j2 t) s t SX " (ȍ ) ,12,1,1 " c (ȍ ) M1 2,,1 ,1 " c ),( 2 ș" s t M (ȍ ) 2,1,2 ,2 " c (ȍ ) M1 ,11,,1 " c (ȍ ) M2 1,,2 ,2 " c (ȍ ) M2 2,,2 ,2 " c Figure 6.4 Contribution of the ℓth wave to the received signal in a MIMO system incorporating dual-antenna arrays In this chapter, we describe the signal model which uses 14 parameters to characterize a single specular propagation path. Then the general Maximum-Likelihood estimates of channel parameters associated with the novel signal model is presented. The main steps in the derivation of the SAGE algorithm are explained. Furthermore, to generalize the conclusions to special cases, an analysis addressing the relevant adaptations for certain specific scenarios such as SIMO, SISO and MISO systems is articulated. In addition, the incoherent initialization technique is introduced which can be included into the estimation scheme based on the derived SAGE algorithm for parameter initialization with low complexity. 6.3.1 Signal model As shown in Figure 6.4, Rx and Tx antenna arrays have dual antenna configuration. Each of the dual antenna transmits/receives signal within two polarizations at the same time. This assumption is based on the realistic experience that antennas can not transmit or receive signals in one polarization only. To differentiate the two polarizations in the underlying model, one of them is called main polarization, specifying the dominant direction of the signal field pattern. The other is correspondingly referred to the complementary polarization. In order to describe the wave polarization, a polarization matrix Aℓ is introduced, which is composed of the complex weights for the attenuations along the propagation paths. The signal model describing the contribution of the ℓth wave to the output of the MIMO system reads s(t; θℓ) = exp(j2πνℓt) c2,1(Ω2,ℓ) c2,2(Ω2,ℓ) αℓ,1,1 αℓ,1,2 αℓ,2,1 αℓ,2,2 · c1,1(Ω1,ℓ) c1,2(Ω1,ℓ) T u(t − τℓ), (6.12) where ci,pi (Ω) denotes the steering vector of the transmitter array (i = 1) with totally M1 entries or receiver array (i = 2) with totally M2 entries, where M1 and M2 are the amount of antennas in the Tx and Rx respectively. Here pi ,(pi = 1, 2) denotes the polarization index. Written in matrix form, (6.12) is reformulated as s(t; θℓ) = exp(j2πνℓt)C2(Ω2,ℓ)AℓC T 1 (Ω1,ℓ)u(t − τℓ), (6.13)
156 Deterministic channel parameter estimation Cyele 1 Cycle 2 2 Mi-l …M网g口…g···☐2…M网a2…网a· Rx Figure 6.5 Timing structure of sounding and sensing windows with C2(2）=【e212c2.22j (6.14 C1(1.）=[c11(1,e）c1.21,] (6.15) (6.16 (9=a(),w(tj7 (6.17刀 The equation(6.13)can be recast as following s(t0)=exp(j2 avit0·{[ae1,1c21(2.)c1(1.d)+a1,2c2,1(2jcH2z(1) +a.21c2.2(2jc.(1)+au2.2c2.2(2.c.2(1.】ut-} (6.18) (6.19y 6.3.2 Channel Sounding Technique both and R.The timing structures used in these sounder systems are quire simlar and are depicred in Figure The mth antennae of Array 1 is active during the sounding windows (6.20) i=1 where i denotes the cycle index and t4m=(位-1)Tg+(m1-1)T. no theorn te m The so-called sensing window 9r(t-t.m.m)2m2=1,,M2,m1=1,,M corresponds to the case where ng is used.the s ndowsand the senin windows)eed merely
156 Deterministic channel parameter estimation 1 Tt 2 M1-1 M1 1 2 M2-1 M2 Cycle 1 Cycle 2 1 2 1 2 1 2 Tg Tsc Tcy t Tr Array 1 Switch 1 Tx Array 2 Switch 2 Rx M2-1 M2 M2-1 M2 M2-1 M2 Figure 6.5 Timing structure of sounding and sensing windows with C2(Ω2,ℓ) = c2,1(Ω2,ℓ) c2,2(Ω2,ℓ) (6.14) C1(Ω1,ℓ) = c1,1(Ω1,ℓ) c1,2(Ω1,ℓ) (6.15) Aℓ = αℓ,1,1 αℓ,1,2 αℓ,2,1 αℓ,2,2 = [αℓ,p2,p1 ] (6.16) u(t) = [u1(t), . . . , uM(t)]T . (6.17) The equation (6.13) can be recast as following s(t; θℓ) = exp(j2πνℓt) · αℓ,1,1c2,1(Ω2,ℓ)c T 1,1 (Ω1,ℓ) + αℓ,1,2c2,1(Ω2,ℓ)c T 1,2 (Ω1,ℓ) +αℓ,2,1c2,2(Ω2,ℓ)c T 1,1 (Ω1,ℓ) + αℓ,2,2c2,2(Ω2,ℓ)c T 1,2 (Ω1,ℓ) u(t − τℓ) (6.18) = exp(j2πνℓt) · X 2 p2=1 X 2 p1=1 αℓ,p2,p1 c2,p2 (Ω2,ℓ)c T 1,p1 (Ω1,ℓ) u(t − τℓ). (6.19) 6.3.2 Channel Sounding Technique Most of the widely-used channel sounders, such as the Propsound Czink (2007), the Medav sounder Richter and Thoma (2005), and the rBECS sounder Yin et al. (2012) are equipped with the radio frequency (RF) antenna switches at both Tx and Rx. The timing structures used in these sounder systems are quite similar and are depicted in Figure 6.5. The m1th antenna element of Array 1 is active during the sounding windows1 q1,m1 (t) = X I i=1 qTt (t − ti,m1 + Tg), m1 = 1, ..., M1, (6.20) where i denotes the cycle index and ti,m1 = (i − 1)Tcy + (m1 − 1)Tt. Here q1,m1 (t) is a real function, with value of 1 or 0 corresponding to the active or inactive moments of the m1th window. Let us define the sounding window vector q1(t) .= [q1,1(t), ..., q1,M1 (t)]T . The so-called sensing window qTsc (t − ti,m1,m2 ), m2 = 1, . . . , M2, m1 = 1, . . . , M1 corresponds to the case where 1Remarks: If another ordering of polarization sounding/sensing is used, the sounding windows q1(t) and the sensing windows q2(t) need merely to be appropriately redefined.
Deterministic channel parameter estimation 157 .The mth Tx antenna is active; The m2th Rx antenna is sensing. where 4m2m1=(-1)+(m1-1+(m2-10, The sensing window for the math Rx dual antenna is given by the real function (6.21) We can define the sensing window vector q2()=21),92.M2 as well as (6.22) Transmitted signal Making use of the sounding window vector (t),we have the explicit transmitted signal u(t)by concatenating the inputs of the M elements of Array 1 u)）=g()u( (6.23) Received signal The signal at the output of Switch 2 can be written as (6.24) with s(t:0)=exp(j2)()C2(D2)A.C(SL)q(t-T)u(t-T). 6.25) Implementing (6.18),we can rewrite s(t:0)=exp(j2mvt)q(t)s(t:0:) 2=11=】 (6.26) (2002a)to incorporate polarization. ng the same approa U(传：T）=q2()q(④'u(t-Ti. (6.27 With this definition,(6.26)can be further written as (6.28) We can also expresss(:)as (6.29 where 8)兰exp(j2mwut)cn(n2.t)Uncp(1） (6.30) an expression similar to (7)in Fleury et al.(2002a)
Deterministic channel parameter estimation 157 • The m1th Tx antenna is active; • The m2th Rx antenna is sensing, where ti,m2,m1 = (i − 1)Tcy + (m1 − 1)Tt + (m2 − 1)Tr. The sensing window for the m2th Rx dual antenna is given by the real function q1,m2 (t) = X I i X M1 m1=1 qTsc (t − ti,m2,m1 ). (6.21) We can define the sensing window vector q2(t) .= [q2,1(t), ..., q2,M2 (t)]T . as well as q2(t) = X I i=1 X M2 m2=1 X M1 m1=1 qTsc (t − ti,m2,m1 ). (6.22) Transmitted signal Making use of the sounding window vector q1(t), we have the explicit transmitted signal u(t) by concatenating the inputs of the M1 elements of Array 1 u(t) = q1(t)u(t). (6.23) Received signal The signal at the output of Switch 2 can be written as Y (t) = X L ℓ=1 q T 2 (t)s(t; θℓ) + r No 2 q2(t)W(t), (6.24) with s(t; θℓ) = exp(j2πνℓt)q T 2 (t)C2(Ω2,ℓ)AℓC1(Ω1,ℓ) T q1(t − τℓ)u(t − τℓ). (6.25) Implementing (6.18), we can rewrite s(t; θℓ) = exp(j2πνℓt)q T 2 (t)s(t; θℓ) = exp(j2πνℓt) · X 2 p2=1 X 2 p1=1 αℓ,p2,p1 q T 2 (t)c2,p2 (Ω2,ℓ)c T 1,p1 (Ω1,ℓ)q1(t) · u(t − τℓ). (6.26) Expression (6.26) is the extension of the first equation in (7) in Fleury et al. (2002a) to incorporate polarization. Following the same approach as in this paper, we define the M2 × M1 sounding matrices U(t; τℓ) = q2(t)q1(t) Tu(t − τℓ). (6.27) With this definition, (6.26) can be further written as s(t; θℓ) = exp(j2πνℓt) X 2 p2=1 X 2 p1=1 αℓ,p2,p1 c T 2,p2 (Ω2,ℓ)U(t; τℓ)c1,p1 (Ω1,ℓ). (6.28) We can also express s(t; θℓ) as s(t; θℓ) = X 2 p2=1 X 2 p1=1 sp2,p1 (t; θℓ), (6.29) where sp2,p1 (t; θℓ) .= αℓ,p2,p1 exp(j2πνℓt)c T 2,p2 (Ω2,ℓ)U(t; τℓ)c1,p1 (Ω1,ℓ), (6.30) an expression similar to (7) in Fleury et al. (2002a)
158 Deterministic channel parameter estimation 6.3.3 SAGE Algorithm Log-likelihood function of the complete/hidden data First,the complete/hidden data is defined as (6.31) It can be shown that the log-likelihood function of given the observation X(t)=(t)reads A(0)2Rs(:0)z(t)dt-Is(t:0)Pdt (6.32) After certain manipulations,it can be shown that A(:)x 2R(of f(B))-IPTeD(Sa..Sh)ot (6.33) The calculation ofG and Gwill be elaborated in the following Computation of G From(6.29)we can write -/会2am (6.34 Inserting(6.30)into(6.34)yields -三三/ano2m0 w.WYx.U.加0 X(r:小 (6.35) with X()denoting the M2x M dimensional matrix with entries X.ma.m(Te,ve)=exp(-j2vdtm(t:T)ze(t)dt =exp(-j2xvet)qz.ma(t)qn.m(t-rt)u"(t-r)ze(t)dt We may notice from the timing structure of the sounding and sensing system that,in each cycle,i.e.i=1,2...,only s principle yie 92m()q1,mt-T)u'(t-T)= u化-T）；tem2m4mam:+T 0 ；otherwise
158 Deterministic channel parameter estimation 6.3.3 SAGE Algorithm Log-likelihood function of the complete/hidden data First, the complete/hidden data is defined as Xℓ(t) = s(t; θℓ) + p βℓ r No 2 q2(t)W(t). (6.31) It can be shown that the log-likelihood function of θℓ given the observation Xℓ(t) = xℓ(t) reads Λ(θℓ; xℓ) ∝ 2R Z s(t; θℓ) ∗xℓ(t)dt | {z } G1 − Z |s(t; θℓ)| 2dt | {z } G2 . (6.32) After certain manipulations, it can be shown that Λ(θℓ; xℓ) ∝ 2R{α H ℓ f(θ¯ ℓ)} − IP Tsc · α H ℓ D˜ (Ω2,ℓ, Ω1,ℓ)αℓ. (6.33) The calculation of G1 and G2 will be elaborated in the following. Computation of G1 From (6.29) we can write G1 = Z X 2 p2=1 X 2 p1=1 sp2,p1 (t; θℓ) ∗xℓ(t)dt = X 2 p2=1 X 2 p1=1 Z sp2,p1 (t; θℓ) ∗xℓ(t)dt. (6.34) Inserting (6.30) into (6.34) yields G1 = X 2 p2=1 X 2 p1=1 Z sp2,p1 (t; θℓ) ∗xℓ(t)dt = X 2 p2=1 X 2 p1=1 Z α ∗ ℓ,p2,p1 exp(−j2πνℓt)c H 2,p2 (Ω2,ℓ)U ∗ (t; τℓ)c ∗ 1,p1 (Ω1,ℓ)xℓ(t)dt = X 2 p2=1 X 2 p1=1 α ∗ ℓ,p2,p1 c H 2,p2 (Ω2,ℓ) Z exp(−j2πνℓt)U ∗ (t; τℓ)xℓ(t)dt | {z } .=Xℓ(τℓ,νℓ) c ∗ 1,p1 (Ω1,ℓ), (6.35) with Xℓ(τℓ, νℓ) denoting the M2 × M1 dimensional matrix with entries Xℓ,m2,m1 (τℓ, νℓ) = Z exp(−j2πνℓt)U ∗ m2,m1 (t; τℓ)xℓ(t)dt = Z exp(−j2πνℓt)q2,m2 (t)q1,m1 (t − τℓ)u ∗ (t − τℓ)xℓ(t)dt. We may notice from the timing structure of the sounding and sensing system that, in each cycle, i.e. i = 1, 2, ..., I, only when the mth 2 transmitter antenna and the mth 1 receiver antenna are active, the product between q2,m2 (t) and q1,m1 (t) gives non-zero result. Applying this principle yields q2,m2 (t)q1,m1 (t − τℓ)u ∗ (t − τℓ) = ( u ∗ (t − τℓ) ; t ∈ [ti,m2,m1 , ti,m2,m1 + Tsc] 0 ; otherwise
Deterministic channel parameter estimation 159 Therefore (p-jant)e (n(t. (6.36) =1,m2m1 Using parameter substitutiontin(6.36)yields Xm2.m(e,以)= exp(-j2m2.m+t)u（传m2m+-Tr6.m2m+td -ep(-j2rvdimm) 1 exp(-j2mvet)u'(ti.mzm+t-Te)re(ti.ma.m+)dt'. (6.37) Here u(t)is a periodic function with period T and ta is integer times of T,therefore u(m+)- u(t).Inserting this expression in(6.37)and substituting twitht,we obtain Xi.mm()= ep(-2mvtm2,m）u'化-Tep(-2mve)zt+tma,mdt Equation(6.35)can be reformulated as =月 Cai2nfn(a） (6.38) beg@=c.x.emd成=nut aVec(Af)=[a,1,a12,at,21,a,2.2 (6.39 [c(2dX(,c1,(,小 fa)= (6.40) c52n2)x(,w)c12n10 we obtain for G GI=ap f(0t) (6.41)
Deterministic channel parameter estimation 159 Therefore Xℓ,m2,m1 (τℓ, νℓ) = X I i=1 Z ti,m2,m1 +Tsc ti,m2,m1 exp(−j2πνℓt)u ∗ (t − τℓ)xℓ(t)dt. (6.36) Using parameter substitution t ′ = t − ti,m2,m1 in (6.36) yields Xℓ,m2,m1 (τℓ, νℓ) = X I i=1 Z Tsc 0 exp(−j2πνℓ(ti,m2,m1 + t ′ ))u ∗ (ti,m2,m1 + t ′ − τℓ)xℓ(ti,m2,m1 + t ′ )dt′ = X I i=1 exp(−j2πνℓti,m2,m1 ) · Z Tsc 0 exp(−j2πνℓt ′ )u ∗ (ti,m2,m1 + t ′ − τℓ)xℓ(ti,m2,m1 + t ′ )dt′ . (6.37) Here u(t) is a periodic function with period Tsc, and ti,m2,m1 is integer times of Tsc, therefore u(ti,m2,m1 + t ′ − τℓ) = u(t ′ − τℓ). Inserting this expression in (6.37) and substituting t ′ with t, we obtain Xℓ,m2,m1 (τℓ, νℓ) = X I i=1 exp(−j2πνℓti,m2,m1 ) Z Tsc 0 u ∗ (t − τℓ) exp(−j2πνℓt)xℓ(t + ti,m2,m1 )dt. Equation (6.35) can be reformulated as G1 = X 2 p2=1 X 2 p1=1 α ∗ ℓ,p2,p1 c H 2,p2 (Ω2,ℓ)Xℓ(τℓ, νℓ)c1,p1 (Ω1,ℓ) ∗ , = X 2 p2=1 X 2 p1=1 α ∗ ℓ,p2,p1 fp2,p1 (θ¯ ℓ) (6.38) where fp2,p1 (θ¯ ℓ) .= c H 2,p2 (Ω2,ℓ)Xℓ(τℓ, νℓ)c1,p1 (Ω1,ℓ) ∗ and θ¯ ℓ = [Ω1,ℓ, Ω2,ℓ, τℓ, νℓ]. Defining αℓ .= Vec(AT ℓ ) = [αℓ,1,1, αℓ,1,2, αℓ,2,1, αℓ,2,2] T , (6.39) f(θ¯ ℓ) = c H 2,1 (Ω2,ℓ)Xℓ(τℓ, νℓ)c1,1(Ω1,ℓ) ∗ c H 2,1 (Ω2,ℓ)Xℓ(τℓ, νℓ)c1,2(Ω1,ℓ) ∗ c H 2,2 (Ω2,ℓ)Xℓ(τℓ, νℓ)c1,1(Ω1,ℓ) ∗ c H 2,2 (Ω2,ℓ)Xℓ(τℓ, νℓ)c1,2(Ω1,ℓ) ∗ , (6.40) we obtain for G1 G1 = α H ℓ f(θ¯ ℓ). (6.41)
160 Deterministic channel parameter estimation Computation of G Inserting(6.30)into (6.32)yields for G2 G2=Is(t:0c)Pdt -/店2 -/小2三onm293ef Vi.p2.(t;Da.eini.ere) =/儿店三exp =1 D52p1a4n） Dcan be further simplified Dmm(7)=Vrni(i.S)"Vimm(S.d =(SU(t:T)eL ('ein (U(tT)em(udt =cr(.UH(Goei.(2小ep(2U6,em,小dt, (6.42) where U(t:)=q2(t)q(t)u(t-).Replacing U(t:e)by this expression in n of(6.42)yields n=c (S)a(t)i(t)e2 (S2)cp:(S2.)qz(t)gn(t)"cim(S)lu(t-re)dt. (6.43) a scalar a scalar Since the product of each pair of terms in(6.43)is a scalar,we can rearrange(6.43)as ncm(z.)q2(t)q(t)c(z)cin(Se)a(t)a(t)cm()u(t-)Pdt, (6.44) where q2()q5因=【2（0,92.2(0，…2,6间j][2,1国，92.20,92，Ma( =diag2.1(t),,g2.u(t)2 =diag21(),92.M(, (6.45) with diagrepresenting a diagonal matrix with entries in the argument and(t(t).Similarly,we have q1d)q1()F=[g1,()91,2（t)91,M,（tj[g1,1()q1,2()91,,(（t) =diagg1.1(),,91.M(tl (6.46)
160 Deterministic channel parameter estimation Computation of G2 Inserting (6.30) into (6.32) yields for G2 G2 = Z |s(t; θℓ)| 2 dt = Z X 2 p2=1 X 2 p1=1 sp2,p1 (t; θℓ) 2 dt = Z X 2 p2=1 X 2 p1=1 αℓ,p2,p1 exp(j2πνℓt) c T 2,p2 (Ω2,ℓ)U(t; τℓ)c1,p1 (Ω1,ℓ) | {z } .=Vℓ,p2,p1 (t; Ω2,ℓ,Ω1,ℓ,τℓ) 2 dt = Z X 2 p2=1 X 2 p1=1 αℓ,p2,p1 Vℓ,p2,p1 (t; Ω2,ℓ, Ω1,ℓ, τℓ) 2 · | exp(j2πνℓt)| 2 | {z } =1 dt = X 2 p ′ 2=1 X 2 p ′ 1=1 X 2 p2=1 X 2 p1=1 α ∗ ℓ,p′ 2 ,p′ 1 αℓ,p2,p1 Z Vℓ,p′ 2 ,p′ 1 (t; Ω2,ℓ, Ω1,ℓ, τℓ) ∗Vℓ,p2,p1 (t; Ω2,ℓ, Ω1,ℓ, τℓ)dt | {z } .=Dp′ 2 ,p′ 1 ,p2,p1 (Ω2,ℓ,Ω1,ℓ,τℓ) . Dp ′ 2 ,p′ 1 ,p2,p1 (Ω2,ℓ, Ω1,ℓ, τℓ) can be further simplified Dp ′ 2 ,p′ 1 ,p2,p1 (Ω2,ℓ, Ω1,ℓ, τℓ) = Z Vℓ,p′ 2 ,p′ 1 (t; Ω2,ℓ, Ω1,ℓ, τℓ) ∗Vℓ,p2,p1 (t; Ω2,ℓ, Ω1,ℓ, τℓ)dt = Z c T 2,p′ 2 (Ω2,ℓ)U(t; τℓ)c1,p′ 1 (Ω1,ℓ) ∗ c T 2,p2 (Ω2,ℓ)U(t; τℓ)c1,p1 (Ω1,ℓ)dt = Z c H 1,p′ 1 (Ω1,ℓ)U H (t; τℓ)c ∗ 2,p′ 2 (Ω2,ℓ)c T 2,p2 (Ω2,ℓ)U(t; τℓ)c1,p1 (Ω1,ℓ)dt | {z } .=η , (6.42) where U(t; τℓ) = q2(t)q1(t) Tu(t − τℓ). Replacing U(t; τℓ) by this expression in η of (6.42) yields η = Z c H 1,p′ 1 (Ω1,ℓ)q1(t) | {z } a scalar q T 2 (t)c2,p′ 2 (Ω2,ℓ) | {z } a scalar c T 2,p2 (Ω2,ℓ)q2(t) | {z } a scalar q1(t) T c1,p1 (Ω1,ℓ) | {z } a scalar |u(t − τℓ)| 2 dt. (6.43) Since the product of each pair of terms in (6.43) is a scalar, we can rearrange (6.43) as η = Z c T 2,p2 (Ω2,ℓ)q2(t)q T 2 (t)c2,p′ 2 (Ω2,ℓ)c H 1,p′ 1 (Ω1,ℓ)q1(t)q1(t) T c1,p1 (Ω1,ℓ)|u(t − τℓ)| 2dt, (6.44) where q2(t)q T 2 (t) = q2,1(t), q2,2(t), ..., q2,M2 (t) T q2,1(t), q2,2(t), ..., q2,M2 (t) = diag[|q2,1(t)| 2 , ..., |q2,M2 (t)| 2 ] = diag[q2,1(t), ..., q2,M2 (t)], (6.45) with diag[...] representing a diagonal matrix with entries in the argument and |q2,m2 (t)| 2 = q2,m2 (t). Similarly, we have q1(t)q1(t) T = q1,1(t), q1,2(t), ..., q1,M1 (t) T q1,1(t), q1,2(t), ..., q1,M1 (t) = diag[q1,1(t), ..., q1,M1 (t)]. (6.46)
