Pilot-Aided Multicarrier Wireless Channel
Estimation via MMSE Polynomial Interpolation
Kun-Chien Hung and David W. Lin
Department of Electronics Engineering and Center for Telecommunications Research National Chiao Tung University, Hsinchu, Taiwan 30010, ROC
E-mails: hkc.ee90g@nctu.edu.tw; dwlin@mail.nctu.edu.tw
Abstract—Orthogonal frequency division multiplexing (OFDM) and multiple access (OFDMA) signal receivers often employ pilot-aided channel estimation. For it, an often considered technique is frequency-domain polynomial interpolation, due to its simplicity. However, the performance of polynomial interpolators suffers in channels with large delay spreads due to modeling error. The problem can be remedied by including a linear phase to the interpolator. In this paper, we derive a method to estimate the optimal phase shift that minimizes the mean-square channel estimation error. We further consider adaptive selection of the interpolation order for best performance. As a practical application, we adapt the proposed channel estimation technique to Mobile WiMAX downlink transmission and examine the resulting performance.
I. INTRODUCTION
Polynomial interpolation is a simple and well-established approach to channel estimation in pilot-aided multicarrier communication [1]. In this approach, one may first estimate the channel responses at the pilot subcarriers by some means and then interpolate between them to obtain the responses at the intervening data subcarriers. Extrapolation is needed if some data subcarriers lie beyond the frequency range covered by the two outermost pilot subcarriers. But we shall refer to the overall approach as interpolation for convenience.
Of all polynomial interpolation schemes, linear (or first-order) interpolation is one of the simplest. However, linear interpolation may suffer great performance loss in poor symbol synchronization [2] or under a large channel delay spread [3]. As a result, some have proposed to alleviate the problem by resorting to high-order interpolation [1]. The way how inter-polation order influences the channel estimation performance can be appreciated via its effect in the time domain. Specifi-cally, polynomial interpolation is effectively a linear filtering operation [4]. Since convolution in the frequency domain corresponds to multiplication in the time domain, different orders of interpolation in the frequency domain correspond to different kinds of windowing on the channel impulse response. Fig. 1 is a conceptional illustration of the situation. Linear interpolation (dashed line) has a more tapered window than quadratic interpolation (dash-dot line) and hence results in a greater distortion of the channel response. However, if one can shift the window corresponding to linear interpolation by some amount to better capture the time range of the significant This work was supported by the National Science Council of R.O.C. under grant no. NSC 96-2219-E-009-003.
channel response samples (solid line), then the performance could improve. Higher-order interpolation is similar.
Three questions arise. First, how to effect the window shift in pilot-aided channel estimation where, instead of the channel impulse response, one only has some frequency response samples at the pilot subcarriers to work with? To answer, note that time shift of a waveform amounts to imposing a linear phase shift (or rotation) on its frequency spectrum. But to realize the phase rotation would cost one complex multiplication per subcarrier. However, if the required time shift is an integer multiple of samples, then the effect can be achieved equivalently and simply by shifting the received signal circularly in the opposite direction by that amount before feeding it into the discrete Fourier transform (DFT).
The second question is how to obtain the optimal window shift for given interpolation order. Hsieh and Wei [1] adopt the single frequency estimators proposed in [5]. However, for channel estimation this is not optimal in the mean-square error (MSE) sense, as later derivation will show. We will derive the optimal window shift that achieves minimum MSE (MMSE) in channel estimation for arbitrary polynomial order.
Thirdly, subsequent analysis will show that the MMSE in channel estimation depends on the amount of channel noise. In high channel noise, higher-order interpolation may be affected so more adversely than lower-order interpolation as to yield worse channel estimates. Therefore, it should be desirable to make the interpolation order adaptive according to the channel condition to attain the best possible channel estimation performance. We will propose a scheme for this.
In what follows, we first consider how to find the MMSE window or phase shift under given polynomial order in Section II. Then, as a practical application, in Section III we adapt the proposed technique to suit Mobile WiMAX downlink spec-ifications and examine the resulting performance. Adaptive selection of interpolation order is investigated in this context.
II. MMSE WINDOW/PHASESHIFTING
Consider a generic discrete-time, equivalent lowpass multi-path Rayleigh wireless channel with its impulse response and (normalized) frequency response given by
h(t) = L l=1 αlδ(t − τl), H(f) = L l=1 αlWτlf, (1)
−800 −60 −40 −20 0 20 40 60 80 0.2 0.4 0.6 0.8 1 1.2 1.4 Sample Time Amplitude Response Linear Interpolation Quadratic Interpolation Shifted Linear Interp. Channel Response
Fig. 1. Comparison of different ways of interpolation in terms of equivalent time-windowing effects.
where L is the number of paths, αl is the complex Gaussian gain of the lth path, τl is the associated path delay, and W = exp(−j2π/Ns) with Ns being the DFT size used in multicarrier transformation. Both t and τl are integers whose unit is the sampling period. And H(f ) has period Ns.
Consider comb-type pilot allocation and let the mth pilot subcarrier be located at frequency p0 + mF where F is the spacing of pilot subcarriers. Let ˆH(p0+ mF ) be some estimated channel responses at the pilot locations; how they are estimated is of no concern at this point. Then conventional polynomial interpolation between two pilot locations p and p + F can be written as ˆ H(p + k) = N n=0 Cn,kH(p + xnˆ ) (2) where 0 < k < F , N is the interpolation order, xn defines the nth pilot location used in interpolation, and Cn,k is the corresponding interpolation coefficient. The interpolation coefficients are real and are well-known [6] to be given in the Vandermonde form by Ck [C0,k C1,k · · · CN,k] =1 k · · · kN V (k) ⎡ ⎢ ⎢ ⎢ ⎣ 1 x0 · · · xN0 1 x1 · · · xN1 .. . ... . .. ... 1 xN · · · xNN ⎤ ⎥ ⎥ ⎥ ⎦ X −1 (3)
or in the Lagrange form by
Cn,k = N m=0,m=n k − xm xn− xm. (4)
Now consider phase-shifted interpolation corresponding to τc samples of window shift. The interpolator coefficients are
given byWτc(k−xn)Cn,k, yielding ˆ H(p + k) = N n=0 Wτc(k−xn)Cn,k H(p + xnˆ ) = Wτc(p+k) N n=0 Cn,k W−τc(p+xn)H(p + xnˆ ). (5)
Interpreted, the right-hand side (RHS) of (5) outlines a pro-cedure for channel estimation under a given τc as follows: 1) obtain channel estimates at the pilot subcarriers by some means, 2) phase-rotate the above channel estimates by an amount corresponding to time advance byτc by carrying out the computation indicated by the last bracket in (5), 3) perform conventional polynomial interpolation as indicated by the last summation in (5), and 4) phase-derotate the channel estimates to undo the earlier time advance by premultiplication with Wτc(p+k). As indicated previously, when τc is an integer, a
simple time-domain equivalent can be used to replace steps 2 and 4. But the theory is unaffected. Further, though we have limited the treatment to comb-type pilot allocation, the results can be readily extended to noncomb-type pilot allocation. But the equations become somewhat cumbersome.
Now we turn to the key issue of determining τc. For this we first find how the channel estimation MSE depends on it. A. Dependence of Channel Estimation MSE on Window Shift For convenience, let ¯H(f ) be the frequency response asso-ciated with theτc-early channel, i.e.,
¯ H(f ) W−τcfH(f ) = L l=1 αlWΔτlf (6) where Δτl= τl− τc, and let
ˆ¯ H(f ) = W−τcfH(f ).ˆ (7) Then ˆ¯ H(p + k) = N n=0 Cn,kH(p + xnˆ¯ ). (8) The MSE in ˆH(p + k) is the same as that in ˆ¯ H(p + k) as they are related by mere phase rotation. Ignore the effect of channel noise and consider only the modeling error for now, so that ˆH(p + xn¯ ) = ¯H(p + xn). The Nth-order Taylor series expansion of the τc-advanced channel response at frequency p + k about a pilot subcarrier p is given by
¯ H(p + k) = N n=0 ¯ H(n)(p)kn n! + RN(p + k) (9)
where ¯H(n)(p) is the nth derivative of ¯H(p) and RN(p + k) is the remainder term which has an integral form given by RN(p + k) = k 0 ¯ H(N+1)(p + u) N ! (k − u) Ndu = L l=1 αl −j2πΔτl Ns N +1 WΔτlp k 0 WΔτlu N ! (k − u) Ndu.
If k is small enough compared to the coherence bandwidth such that WΔτlu≈ 1 for 0 ≤ u ≤ k, then
k 0 WΔτlu N ! (k − u) n d u ≈ k N +1 (N + 1)! (10) and thence RN(p + k) ≈ k N +1 (N + 1)! L l=1 αl −j2πΔτl Ns N +1 WΔτlp ξN(p) . (11) Therefore, the channel estimation error atp + k is given by
¯ H(p + k) − V (k)X−1H(p)ˆ¯ = RN(p + k) − V (k)X−1RN(p) ≈ 1 (N + 1)! kN +1− V (k)X−1xξN(p) (12) where ˆ¯ H(p) = ˆ¯ H(p + x0), ˆ¯H(p + x1), · · · , ˆ¯H(p + xN) T , RN(p) = [RN(p + x0), RN(p + x1), · · · , RN(p + xN)]T, x =xN +10 , xN +11 , · · · , xN +1N T . (13)
It can be shown that
G(k) kN +1− V (k)X−1x = N n=0
(k − xn) . (14) Consequently, the desired relationship between the MSE and τc is given by σe2(k) = 1 [(N + 1)!]2G2(k) |ξN(p)|2 = 1 [(N + 1)!]2G2(k) (2π)2N+2 Ns2N+2 L l=1 |αl|2Δτ2N+2 l σ2 ξ(N,τc) (15)
where · denotes averaging over all pilot subcarriers. B. Estimation of Optimal Window/Phase Shift
From the above, the channel estimation MSE is minimized by minimizing σξ2(N, τc). For this, note that since
¯ H(N+1)(f) = L l=1 αl(−j2πΔτl/Ns)N +1WΔτlf, (16) by Parseval’s theorem we get
σξ2(N, τc) = | ¯H(N+1)(f)|2 (17) where · denotes averaging over all frequencies. Thus the minimization of σξ2(N, τc) can be replaced by that of | ¯H(N+1)(f)|2. While ¯H(N+1)(f) is usually not available, either, it can be approximated by theN +1st-order difference. For example, by using the N + 1st forward difference, it can be shown that the cost function is approximately given by
σξ2(N, τc) ≈ 1 F2N+2 [PNR0+ 2IN(θ)] (18) where PN = 2N +1(2N + 1)!!/(N + 1)!, θ = 2πF τc/Ns, Rn = ˆH(p + nF ) ˆH∗(p), and IN(θ) = N +1 m=1 (−1)mAN,me−jmθRm , (19) withm!! = m(m − 2)(m − 4) · · · x where x = 2 for even m andx = 1 for odd m and
AN,m= 2
N +1(2N + 1)!!(N + 1)!
(N + 1 − m)!(N + 1 + m)!. (20) Hence an estimate of the optimal window shift is
ˆτc= Ns· arg min
θ IN(θ)/2πF. (21)
For linear, quadratic and cubic interpolators, we have I1(θ) = {e−j2θR2− 4e−jθR1},
I2(θ) = {−e−j3θR3+ 6e−j2θR2− 15e−jθR1},
I3(θ) = {e−j4θR4− 8e−j3θR3+ 28e−j2θR2− 56e−jθR1}. One good way to find ˆτc is to search over [0, 2πF τmax/Ns] for value of θ that minimizes IN(θ).
III. APPLICATION TOWIMAX DOWNLINK
The Mobile WiMAX system [7], [8] furnishes an interesting example to test the performance of the proposed technique. In its downlink (DL), the subcarriers in an orthogonal frequency-division multiple access (OFDMA) symbol are divided into “clusters” that contain 14 consecutive subcarriers each. Alter-nating patterns of pilot subcarriers are placed in temporally successive symbols. Fig. 2 illustrates the cluster structure in successive symbols, where the dark circles indicate pilot subcarriers. For convenience, in this section letH(s, k) denote the channel response at the kth subcarrier of some cluster in symbol s. We now tailor the proposed channel estimation technique to suit the WiMAX DL signal structure and inves-tigate the resulting channel estimation performance, assuming perfect carrier frequency and symbol timing synchronization. We also propose a method to adaptively select the interpolation order for minimizing the channel estimation MSE.
A. Channel Estimator Design
There are only two pilot subcarriers per cluster and their spacing varies with symbol index. In a typical system with 10 MHz bandwidth andNs= 1024, the adjacent subcarriers are approximately 10 kHz apart. Considerations over coherence bandwidth and interpolator property lead to the following method of channel estimation. For convenience, we describe it for symbols 2 and 3 (enclosed in the dashed box in Fig. 2) together. To save computation, we estimate ˆτc only once per DL subframe based on the preamble symbol that begins the subframe [7], using the method of Section II-B. In addition, experience shows that an accuracy in ˆτcof 2 to 6 samples can be enough. Therefore, we let ˆτc be an integer.
1) Circularly shift the signal samples in symbols 1 to 4 by −ˆτc before taking their DFT. This is equivalent to having a channel with phase-rotated frequency response
Fig. 2. Cluster structure in Mobile WiMAX downlink and corresponding channel estimation method.
¯
H(s, k) = W−ˆτcf (k)H(s, k) where f (k) denotes the normalized frequency of subcarrierk. (We have used an adapted version of the notations in (6) here.)
2) Do LS channel estimation for pilot subcarriers in sym-bols 1 to 4 (which simply takes the ratio of the received signal value at each pilot subcarrier to the known BPSK pilot value there). This yields, for each pilot with time-frequency index (s, k), a phase-rotated channel estimate
ˆ¯
H(s, k) = W−ˆτcf (k)H(s, k). (We have used an adaptedˆ version of the notations in (7) here.)
3) Linearly interpolate along the time axis to obtain ˆ¯
H(2, 4),H(2, 8),ˆ¯ H(3, 0) andˆ¯ H(3, 12).ˆ¯
4) Perform conventional interpolation along the frequency axis to obtain channel estimates for the remaining sub-carriers in symbols 2 and 3. This corresponds to the operation described in (8), or equivalently, that in the last summation in (5). To minimize the modeling error, in N th-order interpolation we use the N + 1 nearest pilots of each data subcarrier for it.
Note that we do not have to phase-derotate the result as indicated by the premultiplication withWτc(p+k) in the RHS
of (5), for the circular shifting in step 1 makes the phase-rotated channel the target of estimation. In other words, the phase-rotated channel is all that is needed for signal detection. Note also that the temporal linear interpolation in step 3 results in four reference data points per cluster. Hence in step 4 we may employ an interpolation order up to three.
B. Performance Analysis
Four factors contribute to the channel estimation error. They are, in order of their appearance in the channel estimation steps: 1) suboptimality in the estimated window shift ˆτc, 2) the channel noise (introduced in step 2), 3) modeling error due to time-domain interpolation (introduced in step 3), and 4) modeling error due to frequency-domain interpolation (introduced in step 4). Item 1 is difficult to analyze but, fortunately, constitutes a minor contribution in the total MSE. Thus only the other three factors need to be analyzed. We assume that these three kinds of error are uncorrelated.
First, consider the channel noise, assumed additive white Gaussian (i.e., AWGN). It enters the channel estimator com-putation through the LS channel estimation conducted at the pilot subcarriers. Since the pilots are BPSK-modulated, let σn2 denote the estimation noise variance at any pilot. This
estimation noise propagates to other subcarriers in subse-quent temporal and frequency interpolations. Via the tem-poral interpolation, it contaminates ˆH(s, k) where (s, k) ∈¯ {(2, 4), (2, 8), (3, 0), (3, 12)}. Since ˆ¯H(s, k) = [H(s − 1, k) +ˆ¯
ˆ¯
H(s + 1, k)]/2 where H(s − 1, k) andˆ¯ H(s + 1, k) containˆ¯ independent noise, the noise variance in ˆH(s, k) is given by¯
σw2(s, k) = 1 2 2 σ2n+ 1 2 2 σn2= 1 2 σ2n. (22) Via the frequency interpolation, it results in an estimation noise variance σ2w(s, k) = N n=0 Cn,k2 σ2w(s, xn) (23) wheres ∈ {2, 3}, k ∈ {1, 2, 3, 5, 6, 7, 9, 10, 11, 13}, and xn ∈ {0, 4, 8, 12}, with Cn,k indexed similarly to (2). Hence the average noise variance over all data subcarriers is given by
σw,12 = 0.5130σn2, σ2w,2= 0.6699σ2n, σ2w,3= 0.8121σ2n, where the second subscript toσ denotes interpolation order.
Secondly, consider the modeling error in temporal interpo-lation. Its mean-square value is defined by
σD2(s, k) = E
¯H(s,k) −12 ¯H(s − 1, k) + ¯H(s + 1, k)
2
for (s, k) ∈ {(2, 4), (2, 8), (3, 0), (3, 12)}. Assume Rayleigh faded paths. Then [3], [9]
σD2(s, k) = L l=1 E|αl|23 2 −2J0(2πfl) + 1 2 J0(4πfl) ≈3π24 L l=1 fl4E|αl|2 σ2D (24) whereflis the peak Doppler shift of pathl times the OFDMA symbol period, J0(·) is the Bessel function of the first kind of order 0, and the approximation is obtained by expanding the Bessel function into a second-order Taylor series. Assume that the channel responses at different subcarriers are uncor-related. Then a relation similar to (23) exists concerning the propagation of the temporal modeling error in the frequency domain via channel estimation step 4, and we get the average MSE over all data subcarriers as
σd,12 = 0.4531σ2D, σ2d,2= 0.5577σD2, σ2d,3= 0.6525σD2, where the second subscript toσ again gives the interpolation order.
Finally, consider the modeling error in frequency interpola-tion. From (15), straightforward numerical calculation yields the following average MSE over data subcarriers:
σ2i,1= k1σξ2(1, τco), σ2i,2= k2σξ2(2, τco), σ2i,3= k3σ2ξ(3, τco), where k1 = 2.6458, k2 = 12.8125, k3 = 93.0820, and the second subscript toσ once more gives the interpolation order.
5 10 15 20 25 30 35 40 −40 −35 −30 −25 −20 −15 −10 −5 Pilot SNR (dB) NMSE (dB) Simulation: Linear Analysis: Linear Simulation: Quad Analysis: Quad. Simulation: Cubic Analysis: Cubic
Fig. 3. NMSE of channel estimation in WiMAX DL transmission over SUI-4 channel at 100 km/h mobile speed with different orders of interpolation.
Putting all together, we obtain the overall average channel estimation MSE as
σ2H,j = σ2w,j+ σd,j2 + σi,j2 (25) for interpolation orderj ∈ {1, 2, 3}.
C. Simulation Results
We simulate a system with carrier frequency = 2.5 GHz, bandwidth = 10 MHz, DFT size = 1024, and cyclic prefix length = 128. And we let a DL subframe contain 24 OFDMA symbols following the preamble. Below we present some results under the SUI-4 power-delay profile (PDP) [10] with block-static fading at a rate corresponding to 100 km/h of mobile speed. The channel power profile is [0, −4, −4] (in dB) and the delay profile is [0, 14, 36] (in sample periods).
Fig. 3 show some results. The approximate analysis matches the simulation results very closely. (For channels with larger delay spread, the approximate analysis will be less accurate, but will still follow the general behavior of the simulation results.) Compared to quadratic and cubic interpolations, linear interpolation exhibits a higher MSE floor (due to frequency interpolation) in high pilot SNR, but it performs somewhat better in low SNR where the AWGN effect is more prominent. D. Adaptive Selection of Interpolation Order
That lower-order interpolators incur smaller error due to AWGN whereas higher-order ones have smaller modeling error suggests adaptive selection of the interpolator order. This can be accomplished if we can estimateσ2n, σξ2(N, τco), and σD2 from the preamble symbol and use the results to predict the MSE via (25) and determine the interpolation order.
It is relatively easy to estimateσn2 using the null subcarriers in the preamble symbol. Further σ2ξ(N, τco) can be estimated based on (18), withR0estimated by| ˆH(f )|2−ˆσn2where ˆσ2n is the estimate of σ2n. (Note thatP1= 6, P2= 20, and P3= 70.) But it is hard to estimate σ2
D, which depends on the time-variation of the channel response, using only the preamble
10 15 20 25 30 35 40 −45 −40 −35 −30 −25 −20 −15 −10 Pilot SNR (dB) NMSE (dB) Simul.: Linear Simul.: Quad. Simul.: Cubic Simul.: Adaptive Simul.: Approx. MMSE Simul.: Exact MMSE
Fig. 4. Comparison of adaptive interpolation with fixed interpolation over the SUI-4 PDP under 100 km/h mobile speed.
symbol. Fortunately, it is the least dominating term of all and is thus disregarded.
Fig. 4 compares the MSE performance of adaptive interpo-lation with that of fixed interpointerpo-lation in the SUI-4 channel. The performance of adaptive interpolation is nearly optimal. For comparison, we also plot the performance of 4-tap MMSE channel estimation obtained with the exact channel correlation function (marked “exact MMSE”) and that with the correlation function corresponding to a uniform PDP of length equal to the cyclic prefix (marked “approx. MMSE”). The exact MMSE estimator performs better than polynomial interpolation, but is impractical. On the other hand, the adaptive interpolation scheme outperforms the approximate MMSE estimator.
REFERENCES
[1] M.-H. Hsieh and C.-H. Wei, “Channel estimation for OFDM systems based on comb-type pilot arrangement in frequency selective fading channels,” IEEE Trans. Consumer Electron., vol. 44, no. 1, pp. 217– 225, Feb. 1998.
[2] J. Park et al., “Performance analysis of channel estimation for OFDM systems with residual timing offset,” IEEE Trans. Wireless Commun., vol. 5, no. 7, pp. 1622–1625, July 2006.
[3] K.-C. Hung and D. W. Lin, “Optimal delay estimation for phase-rotated linear interpolative channel estimation in OFDM and OFDMA systems,”
IEEE Signal Processing Lett., vol. 15, pp. 349–352, 2008.
[4] R. W. Schafer and L. R. Rabiner, “A digital signal processing approach to interpolation,” Proc. IEEE, vol. 61, no. 6, pp. 692–702, June 1973. [5] S. Kay, “A fast and accurate single frequency estimator,” IEEE Trans.
Acoust. Speech Signal Processing, vol. 37, no. 12, pp. 1987–1990, Dec.
1989.
[6] E. K. Blum, Numerical Analysis and Computation: Theory and Practice. Reading, Mass.: Addison-Wesley, 1972.
[7] IEEE Stds. 802.16e-2005 and 802.16-2004/Cor1-2005, IEEE Standard
for Local and Metropolitan Area Networks — Part 16: Air Interface for Fixed and Mobile Broadband Wireless Access Systems — Amendment 2: Physical and Medium Access Control Layers for Combined Fixed and Mobile Operation in Licensed Bands and Corrigendum 1. New York:
IEEE, Feb. 28, 2006.
[8] WiMAX Forum, Mobile WiMAX — Part I: A Technical Overview and
Performance Evaluation. WiMAX Forum White Paper, Aug. 2006.
[9] T. S. Rappaport, Wireless Communications Principles and Practice, 2nd
ed. Upper Saddle River, New Jersey: Prentice Hall, 2002.
[10] V. Erceg et al., “Channel models for fixed wireless applications,” standards contribution IEEE 802.16.3c-01/29r1, Feb. 23, 2001.
