Blind Channel Estimation Using Fractional Sampling - CiteSeerX
1
Blind Channel Estimation Using Fractional Sampling Sarod Yatawatta, Student Member, IEEE, Athina P. Petropulu, Senior Member, IEEE, and Riddhi Dattani, Student Member, IEEE,
Abstract— We consider the problem of blind estimation of a communication channel by oversampling the channel output. We propose a non-parametric approach which, based on the cyclic spectrum, finds the channel phase response without the need of phase unwrapping nor channel length information. For bandlimited channels, the cyclic spectrum is not supported for the full bandwidth. For this case, we propose an approximation for the discretized phase of the cyclic spectrum, which, under certain conditions, results in a simpler channel estimation method. The proposed approach is applied to simulated data and real recordings, and is compared to existing methods. Index Terms— Blind channel estimation, cyclostationary inputs, cyclic spectrum.
I. I NTRODUCTION
T
HE communications channel sets a limit to the transmission rate of the communication system. Signal pulses are broadened in time as they travel through the channel (multipath effect), leading to intersymbol interference (ISI) at the receiver. In order to remove the effect of ISI, prior knowledge of the channel is required. Generally this information is obtained by transmitting training symbols. However, this approach implies waste of bandwidth, which can be quite significant in cases where the channel is time varying. An alternative approach is blind channel estimation and equalization. Within short time intervals, estimation of the communication channel based on the received data can be viewed as a blind linear time invariant (LTI) system estimation problem. In general, for the estimation of an LTI system excited by an unknown stationary non-Gaussian input, based on the system output, one has to use higher-order statistics (HOS) of the output sampled at baud rate. When the input is cyclostationary, which is the case of communications signals, system estimation can be carried out using second-order cyclic statistics of the oversampled output [1], [2]. Based on this premise, several methods have been proposed for blind channel estimation (for a review of existing results the reader is referred to [8]). This work was supported by ONR under grant ONR-N00014-20-1-0137 The authors are with The Department of Electrical & Computer Engineering, Drexel University, Philadelphia, PA 19104
They can be categorized as time-domain and frequencydomain, with reference to the transform domain in which the problem is formulated. A serious deficiency in time domain-methods ([2], [7]) is their sensitivity to channel order estimation errors. Most often in practice the channel (the combined effect of transmit-receive filters and multipaths) is a continuous function of time and its magnitude will not suddenly drop to zero at one point. The notion of channel length depends on the application. Some methods have been proposed to overcome the need for channel length information [11], [14], [12], yet this could still pose a problem when the channel is time varying. On the other hand, frequency-domain methods do not require channel length, but cannot be applied to band-limited channels. In that case, the phase of the full spectrum cannot be recovered because the cyclic spectrum has limited support [5]. In practice, most channels are bandlimited, both by the physical nature of the channel and by the pulse shaping employed. In this paper we propose a novel frequency-domain blind channel estimation approach. The channel phase frequency response is obtained in closed form in terms of the appropriately discretized cyclic phase spectrum. For the case of bandlimited channels with a small rolloff, we propose an approximate expression for the channel phase response, which involves low computational complexity. The resulting channel estimate will always contain an error, however the error will be negligible if the channel has negligible energy within the stopband. Despite this bias, the proposed approach can still be of great practical use since it involves low cost. As a matter of practice, choosing the correct algorithm for a particular purpose does not depend on the accuracy of the result only, but also on economy in hardware implementation. With the introduction of wireless communication capabilities to most electronic devices in the offing, the proposed approach could provide a viable alternative in reducing cost. The paper is organized as follows. Section II provides the problem formulation and summarizes relevant literature. Section III details the proposed approach. In section III.A, the general method is described and in section III.B, the simplification applicable to bandlimited channels is given.
2
Section IV provides results on the performance of the proposed methods based simulations and also based on real data recordings. Finally, concluding remarks are given in Section V. II. P ROBLEM F ORMULATION Let us consider the baseband representation of a digital data communication system. Assuming that the channel is LTI and BIBO (Bounded Input producing Bounded Output) stable, the received signal is given by: x(t) =
∞ X
sν h(t − νT ) + w(t), ν ∈ Z
(1)
ν=−∞
where sν is the transmitted signal, which is assumed to be zero-mean i.i.d. with unit variance; h(t) is the channel impulse response; T is the symbol interval, and w(t) is stationary, zero-mean, white noise with variance No and is independent of sν . The signal x(t) is cyclostationary with cyclic period T [1]. Sampling x(t) with sampling interval ∆ = Tp , p ∈ Z, yields: 4
x(n) = x(n∆) =
∞ X
sν h(n − νp) + w(n), n ∈ Z (2)
ν=−∞
where h(n) = h(n∆), w(n) = w(n∆). The discrete time signal x(n) is a cyclostationary process with cyclic period p. Cyclostationarity implies that the autocorrelation function is periodic with period p. The Fourier series coefficients of the autocorrelation are referred to as cyclic autocorrelation, Rxkα (ν), and the Fourier Transform of those as the cyclic spectrum, Sxkα (ejω ) [1], i.e., Rxkα (ν) =
p−1 X
Rx (n + ν, n)e−j2πkαn , α =
n=0
1 , p
(3)
where ψ(ω), φ(ω) are the phases of Sxkα (ω) and H(ejω ), respectively. It has been shown in [4], that the above equation cannot recover the phase response of an arbitrary channel, since when evaluated for different ω’s it always leads to an underdetermined system of equations. Viewing ψ(ω), φ(ω) as Fourier Transforms of sequences Ψ(τ ) and Φ(τ ), respectively, it holds [5] that: Φ(τ ) =
jΨ(τ ) , τ 6= νT, ν ∈ Z. 2sin(πτ /T )
(7)
Since Φ(τ ) is infinite at τ = νT , the phase φ(ω) cannot be obtained from (7). However, if Φ(τ ) is bounded at τ = νT , the missing information can be recovered via interpolation, assuming continuity of Φ(τ ). The case where Φ(τ ) is not bounded, i.e., contains impulses at τ = νT occurs when φ(ω) contains sinusoidal components with frequency νT . This observation led to the identifiability condition [3] stating that discrete channels with zeros uniformly 2π/p spaced on a circle are not identifiable because they result in oscillatory phase. However, even if the channel is identifiable, (7) will be inaccurate for τ ≈ νT . Parametric approaches that obtain the channel parameters from cepstral information obtained for τ away from the points νT were proposed in [5]. III. P ROPOSED M ETHOD We first consider an FIR channel that has a full passband. Later, in section B, we consider the situation where the channel is bandlimited. A. FIR Channels
First we make the following assumption: (A1) The channel satisfies identifiability conditions given in n, ν, k ∈ Z = [3] and [9], i.e., there are no uniformly 2π/p spaced ν=−∞ zeros on a circle and there are no spectral nulls in the where Rx (ν, n) denotes autocorrelation of x(n). non-rolloff passband. It holds that: We next show that an appropriately discretized version kα jω Sx (e ) = H(ejω )H ∗ (ej(ω−2πkα) ) + pNo δ(k), k ∈ Z. of (6) can lead to a solution of the channel phase. (4) Let ω = 2π N l, l = 0, ..., N − 1. Then, This is a key relationship that can lead to the estimation of ψ(l) = φ(l + m) − φ(l − m) (8) the channel frequency response H(ejω ). The magnitude response of the channel can be easily or, recovered from the power-spectrum of x(n). To obtain ψ(l + m) = φ(l + 2m) − φ(l) (9) phase information, we can rewrite (4) as: Sxkα (ejω )
4
Sxkα (ω) = =
∞ X
Rxkα (ν)e−jνω ,
Sxkα (ej(ω+πkα) ), (5) H(ej(ω+πkα) )H ∗ (ej(ω−πkα) ), k 6= 0
from which we get: ψ(ω) = φ(ω + πkα) − φ(ω − πkα), k = ±1, ±2, .... (6)
where m is the closest integer to N kα/2 = N k/2p. Let us fix k to the value k = 1 for maximum support [1]. Based on (9) evaluated for l = 1, ..., N − 1, and taking φ(0) = 0, we can formulate the system of equations: AN,m Φ = Ψ
(10)
3
where Φ = [φ(1) ... φ(N −1)]T , Ψ = [ψ(m+1) ... ψ(m+ (A2) The transition band of the channel is negligible compared to the passband, i.e., the rolloff is very small. N − 1)]T , and ψ(l) = ψ(l + N ). AN,m is a (N − 1) × (N − 1) sparse matrix. It can be (A3) The energy within the stopband of the channel is negligible. viewed as consisting of the columns 2 to N and rows 2 to N of an N × N circulant matrix, B, whose first row (A4) Let Lc = bωc N/2πc (cutoff frequency index), where ωc is the cutoff frequency, then equals: b = [−1, 0T2m−1 , 1, 0TN −2m−1 ]. (11) Lc < N/4. (14) Here 0Tk denotes a row vector containing k zeros. The other (A5) The number of zeros of the channel on the unit circle Nz satisfies, rows of B are circularly shifted versions of b. N − 2Lc A matrix similar to AN,m was also encountered in [13] Nz > L . (15) N in a different context, where it was shown that for N and where L is the channel length. 2m co-prime numbers, it is full rank with determinant equal Assumption (A4) implies that ωc < π/2 rad/cycle. For to one. Since N is a number that we select, we can give it a any system that transmits at or below Nyquist rate [15], the value so that it is co-prime with 2m = 2bN/2pc, thus channel sampled at symbol rate (p = 1) has full bandwidth. guaranteeing full rank for A. The phase vector Φ can then Hence to satisfy this condition we need p ≥ 3. In general, be obtained by solving (10). Ideally, we would like the we have no prior knowledge of Nz to check the validity of phase expression to involve the modulo 2π phases only, so assumption (A5). However, by careful design of the shaping that no phase unwrapping is necessary. It turns out that function, we can ensure that the minimum number of Nz Φ can actually be based on modulo 2π phases, and this is is greater than the one required by (15). Multipaths will due to the special properties of matrix AN,m [13]. It can only add more zeros thus the condition will hold under multipaths if the shaping function itself satisfies it, unless be shown that, the estimate: the channel becomes unidentifiable according to [9]. 4 ˜ b = Let us discretize ψ(ω) in (6) on a grid with spacing (12) Φ A−1 N,m Ψ 2π/N and select N and p so that whenever φ(ω) falls in ˜ is the modulo 2π phase corresponding to Ψ, the passband, φ(ω + 2πkα) falls in the stopband. We take where Ψ satisfies: k = 1 again for maximum support. This can be achieved 2π T b = Φ+φ(0)1(N −1)×1 + ν[1, ..., N −1] +2πK (13) as long as, Φ Lc < m < N/2 − Lc (16) N where K is a vector of size (N − 1) × 1 with integer elements, and ν is an unknown integer. Combining the phase estimate of (12) with the correct magnitude spectrum, results in the channel impulse response with an unknown scalar multiplicative constant, and a circular time shift equal to ν. We should note here that in comparison with the methods of [5],[6], the proposed method does not require phase unwrapping. B. FIR Channels with pulse shaping In this section, we consider the case where the channel is bandlimited by pulse shaping. The objective of pulse shaping is to restrict the spectrum of the signal within a fixed frequency interval [15]. Hence, by its design, the magnitude within the stopband is low compared to the passband. Moreover, (4) will have limited support, thus (9) will result in an underdetermined system of equations. In [5], it has been suggested to use parametric methods to overcome this difficulty. Here we propose a non parametric approach that yields an approximate channel estimate. In addition to assumption (A1) we make the following assumptions:
where m is the closest integer to αN/2 or N/2p. A value for m greater than 1 is guaranteed by N > 4Lc . Proposition 1: Let h(n) be an FIR impulse response of d length L, which satisfies assumptions (A1)-(A5). Let h(n) be the channel impulse response constructed based on the correct magnitude response of h(n), i.e., |H(l)|, and phase d constructed using the estimated phase of the response φ(l), \ cyclic spectrum, ψ(l + m), as: \ −ψ(l + m); 4 d = φ(l) 0 ≤ l ≤ Lc − 1, N − Lc ≤ l ≤ N − 1 0; Lc ≤ l ≤ N − Lc − 1 (17) with N and p satisfying (16). Then, it holds: d h(n) (18) d
=
4
IDF T {|H(l)|ej φ(l) }
=
ah(n + b) ⊗
2π
−2N e−j N ko n 2π
Λ(1 − ej N n ) ! Λ X (2i + 1)N δ n− × + E(n) 2Λ i=0
4
where E(n) =
a N
N −L Xc −1
2π
H(l)ej N ln
− H(l)e
j{ 2π N (b+n)l+π
P(N −1)/Λ−1 i=Nz
a = e−jNz (N/2−m) , b = Λ=b
(19)
l=Lc
N − 2Lc c, Nz − 1
u(l+ko −iΛ)}
!
,
Nz , ko = 2m − Lc , 2
NX z −1 ω π − ωc φ(ω) = Nz (π− )+φt (ω)−π ) u(ω−ωc −i 2 2 Nz − 1 i=0 (23) where u(ω) is the step function. Using (23), we can rewrite φ(ω + 2πα) in discrete form (ω = 2π N l) as:
and ⊗ denotes N point circular convolution. Proof: The channel Frequency Response can be given as: H(ω) =
NY total
each zero on the unit circle. Note that these jumps will not occur within the passband of the channel. Hence, if Nz is large enough, within the stopband, φz (ω) will dominate the overall phase response because all its derivatives become large whenever ω ≈ θi for any i ∈ [0, Nz − 1]. If the phase due to zeros away from the unit circle is φt (ω), the phase of the channel can be given as
φ(l + 2m) (1 − zi z −1 )|z=ejω
(24)
Nz π l + φt (l + 2m) = Nz (N/2 − m) − N NX z −1 −π u(l + 2m − Lc − iΛ)
(20)
i=0
Here the channel has Ntotal zeros, out of which Nz lie on the unit circle.
i=0
Zi r
where l is the discrete frequency index. Under the assumptions made earlier, φ(l) will lie on the passband while φ(l + 2m) will lie on the stopband. Hence the jumps of π will dominate (24) thus we can ignore the contribution from φt (l + 2m) to simplify our analysis. Hence we can write (9) as
φi di ωc − ω r≈1
−ψ(l + m)
(25)
Nz π ≈ φ(l) + l − Nz (N/2 − m) N NX z −1 +π u(l + 2m − Lc − iΛ).
Fig. 1. Calculation of Phase Due to Channel Zeros Corresponding to the Stopband. As ω increases, the phase contribution from these zeros change as if the location of the zeros rotate clockwise.
i=0
Let us refer to Fig. 1 to calculate the phase of the terms of (20) that contain the zeros on the unit circle. As ω increases, the phase contributed by this set of zeros changes as if the zeros themselves are rotated by ω clockwise. Let φi be the phase of (1 − zi e−jω ) where |zi | = 1. It holds: φi (ω) =
π ωc − ω π − ωc + +i , 0 ≤ i < Nz . 2 2 Nz − 1
(21)
The sum of the phases of all terms (1−zi e−jω ) for which |zi | = 1, i.e. φz (ω) will then be equal to Nz (π − ω2 ). For ω ≈ θi , it holds: r 1 dφi (ω) = ((1−r)2 − (1+r)(θi −ω)2 ), ω ≈ θi dω (1 − r)3 2 (22) i (ω) We see that as r → 1, dφdω becomes infinite, thus indicating discontinuity of φi at ω = θi . We see from Fig. 1 that at this discontinuity, the overall change in φi (ω) is π. Hence, the overall phase response is linear with jumps of π at frequencies corresponding to
Then we can re write (17) as: φ(l) + NNz π l − Nz (N/2 − m) PNz −1 + π i=0 u(l + 2m − Lc − iΛ); d ≈ φ(l) 0 ≤ l ≤ L c − 1, N − Lc ≤ l ≤ N − 1 0; Lc ≤ l ≤ N − Lc − 1. (26) d d = IDF T {|H(l)|ej φ(l) Then h(n) } yields:
d = h(n) a N
+
LX c −1
(27)
2π
H(l)ej{ N (b+n)l+π
l=0 N −L Xc −1
PNz −1
2π
H(l)ej{ N (b+n)l+π
i=0
u(l+ko −iΛ)}
P(N −1)/Λ−1 i=Nz
u(l+ko −iΛ)}
l=Lc
+
N −1 X
l=N −Lc
+ E(n)
2π
H(l)ej{ N (b+n)l+π
PNz −1 i=0
u(l−N +1+ko −iΛ)}
5
=
N −1 2π a X H(l)ej{ N (b+n)l} G(l) + E(n) N l=0
with G(l) = 2
M X
rect(
i=0
l + ko − i2Λ ) − 1, M ∈ Z Λ
(28)
where E(n) can be deduced from (27) and rect( Λτ ) is a rectangular pulse of width Λ. We see that the IDFT of G(l) is 2π
g(n) =
2π
e−j N ko n (1 − e−j N Λn ) (1 − ej2πn ) (1 + e−j
2π N Λn
2π
(1 − ej N n )
)
(29)
N N where g(n) 6= 0 for n = 2Λ , 3 2Λ , . . .. By simplification of (29), (18) follows. (E.O.P) d will contain delayed Equation (18) suggests that h(n) copies of the original channel separated be N Λ . If condition (15) holds, these copies will be separated well enough to allow one to extract h(n). Note also that the error E(n) will be small (by assumption (A3)) if the energy within the stopband is negligible, i.e., N −L −1 2 Xc |E(n)| ≤ |H(l)|. (30) N l=Lc
The magnitude of each replica appearing in (27) is determined by the magnitudes of the impulses in (18), which in turn is determined by f (n) =
1 2π
|1 − ej N n |
.
(31)
This expression has its maxima at n = 0 and n = N , is monotonically decreasing as n goes from 0 to N/2 and, is monotonically increasing as n goes from N/2 to N . Hence the first and last replicas will have the highest magnitudes. As a consequence, to extract one copy of the channel d we could do the following. Let the estimate out of h(n) anticipated channel length be Le . We find the segment with maximum energy out of the first Le samples and the last Le samples. After extracting the maximum energy sequence, we obtain a channel estimate with a scaling ambiguity and an unknown delay. IV. S IMULATION R ESULTS Simulations were done to verify the performance of the proposed methods and to compare with other methods. In section A, we give the results obtained for FIR channels with full bandwidth support, and in section B we give results for bandlimited channels. Finally, in section C, we give results of the application of the method in section III.B to real measurements.
The basic procedure is P as follows. First, we generated ∞ a QAM signal x(n) = ν=−∞ sν h(n − νp) + w(n). The input symbols {sν } were from an i.i.d. QAM source. The noise w(n) was stationary, white and Gaussian. The Spectral Correlation Density (SCD) or cyclic spectrum was estimated as follows. First we estimated the unbiased autocorrelation as: N −1 X bx (n + ν, n) = 1 R x(n + ν + pl)x∗ (n + pl), n, ν ∈ Z. N l=0 (32) bxkα (ν) Next, we estimated the cyclic autocorrelation R using equation (4). If some estimate of the channel length, ˆ is available, then we can window R bxkα (ν) by a i.e., L ˆ rectangular window of length 2L. Finally, we estimated the SCD by taking an N point FFT, i.e.,
bxkα (ν)), α = 1 , k = 1. Sbxkα (l) = FFT(R p
(33)
b The phase of the channel was estimated using ψ(l), the modulo 2π phase of Sbxα (l), according to methods proposed in sections III.A and III.B. In order to make b ψ(l) symmetric, it was shifted circularly by b−N/2pc. (0) The power spectral density, Sbx (l) was used for amplitude estimation. b b = |Sbx(0) (l)|1/2 exp(j φ(l)) H(l)
(34)
Finally, the impulse response was computed as: b b h(n) = IFFT(H(l)).
(35)
The performance was measured using and the Overall Normalized Mean Square Error (ONMSE) and also Bit Error Rate (BER). We tested a total of Mh different channels. For each channel we performed Mc Monte Carlo runs corresponding to independent input realizations. The ONMSE was defined as: Mh Mc P L j 2 bj 1 X 1 X n=0 |hi (n) − h (n)| ONMSE = PL j 2 Mh j=1 Mc i=1 n=0 |h (n)| (36) ˆ j is the estimate of the i-th Monte Carlo run of the where h i j-th channel. Based on each channel estimate we implemented a zero forcing equalizer and obtained the input estimate. The BER rate was computed for each channel and then averaged over the different channels tested. A. FIR Channels Each channel was taken to be FIR of length 8, and its coefficients were generated randomly. For each channel we performed Mc = 25 Monte Carlo runs, using the estimation method of Section III.A. Finally the ONMSE and BER
6
B. FIR Channels with pulse shaping We generated channels of the form h(t) = a1 c(t + 1.5T, 0.11) + a2 c(t − T1 , 0.11) (37) + a3 c(t − 1.5T, 0.11) W6T (t)
where a1 , a2 , a3 are normally distributed random numbers with mean zero and variance one and T1 is uniformly distributed in [−T, T ]. The shaping function c(t − τ, β) is a
−0.3
10
with windowing without windowing
−0.4
ONMSE
10
−0.5
10
Proposed Cepstrum based
−0.6
10
−0.7
10
0
5
10
15
20
25
30
35
SNR/(dB)
(a) 0
10
with windowing without windowing Proposed Cepstrum based
BER
were obtained as the average over Mh = 100 different channels. For comparison purposes, we also implemented the nonparametric method (cepstrum based) of [6], where the phase of the channel is estimated using the Fourier Series expansion of the unwrapped phase of the SCD. bxα (ν) was obtained with a rectangular The estimate of R windowing of length 2L as well as without windowing. The obtained estimate, usually contains more than one replicas of the channel response. Although the precise reason for such behavior needs to be further investigated, simulations suggest that this is a result of phase errors and the way they propagate to the final solution. For large N , i.e., N > 3L, there is enough room for the replicas to be separated. One could then extract the channel by identifying the segment with maximum energy, with length equal to the anticipated channel length, or length equal to the window size used in the equalizer. There will be at least one such disjoint segment in the estimate where the energy will be distinctly greater than the rest of the estimate. In the proposed method we chose N = 131 and p = 4. In the simulations we extracted the segment with maximum energy of length L, and we have given results in Fig. 2. ONMSE results are shown in Fig.2(a). From Fig. 2(a), bxα (ν), was not we see that the windowing used to estimate R that critical in the final estimate. Since the window requires some knowledge of channel length information, the above observation confirms that the methods tested here are not very sensitive to the channel length information. Corresponding BER results are shown in Fig. 2(b). Here we used channel length knowledge in order to compute the Zero Forcing equalizer. The results were obtained for 1000 symbols. The above results suggest that the proposed approach is only slightly better than that of [6] in terms of ONMSE and has same performance in BER. However, in terms of complexity, the proposed approach is less intensive than that of [6], since it requires no phase unwrapping. The matrix A is independent of the data, hence we only need its inverse and moreover, the inverse only has ones and zeros. Hence solving equation (12) does not involve multiplications. In contrast, the cepstrum based method [6] requires phase unwrapping and two additional FFTs.
−1
10
−2
10
0
5
10
15
20
25
30
35
SNR/(dB)
(b) Fig. 2. (a) ONMSE and (b) BER for cepstrum based estimates and proposed estimates given in section III.A. We generated 100 random FIR channels of length 8 and for each channel, obtained the estimates for 25 Monte Carlo runs. 1000 symbols were used in channel estimation and perfect channel length knowledge was assumed. In (a), we have given the ONMSE results where both a rectangular window of length 2L was used, and a window was not used.
raised cosine [15] with delay τ and rolloff β, and W6T (t) is a rectangular window of width 6 symbol intervals T . We here considered Mc = 25, Mh = 100, N = 128 and p = 4. We compared the method proposed in section III.B with the high complexity subspace based method of [7]. In Fig. 3 we have given ONMSE and BER results obtained with perfect channel length knowledge, and in Fig. 4 we have given results where the channel length was overestimated by 1 symbol (4 samples). When perfect channel length knowledge is available, the method of [7] performs exceptionally well for high SNR, while the proposed method performs slightly better for low SNR. On the other hand, the proposed approximate method outperforms the former when there is a slight mismatch in channel length estimation.
7
0
0
10
10
Proposed Subspace
Proposed Subspace
−1
ONMSE
ONMSE
10
400 Symbols 1000 Symblos
−2
10
400 Symbols 1000 Symblos
−3
10
0
−1
10
−2
5
10
15
20
25
30
10
35
0
5
10
SNR/(dB)
15
20
25
30
35
SNR/(dB)
(a)
(a)
0
0
10
10
Proposed Subspace
Proposed Subspace −1
BER
BER
10
−2
10
−1
10
400 Symbols 1000 Symblos
400 Symbols 1000 Symblos
−3
10
−4
10
0
−2
5
10
15
20
25
30
35
SNR/(dB)
10
0
5
10
15
20
25
30
35
SNR/(dB)
(b)
(b)
Fig. 3. (a) ONMSE and (b) BER for the method proposed in section III.B and for the subspace method. We generated 100 channels with pulse shaping and for each channel, 25 Monte Carlo runs were taken at each SNR. The modulation scheme used was 4 QAM and perfect channel length knowledge was used for equalization.
Fig. 4. (a) ONMSE and (b) BER for the method proposed in section III.B and for the subspace method. We generated 100 channels with pulse shaping and for each channel, 25 Monte Carlo runs were taken at each SNR. The modulation scheme used was 4 QAM and the channel length was overestimated by one symbol interval.
C. Performance with experimental data To check whether the channels given in (38) satisfy (15) we can do the following. If we oversample by 4 (p = 4, T = 1), we see that L = 24. Moreover, Lc ≈ 20, and for N = 128 (15) gives Nz > 16.5. By evaluating the zeros of the channel we find Nz = 18, indeed satisfying (15). Finally, we investigated the effect of rolloff to the error in approximation. We generated 30 channels of the form (38) with variable rolloff. We varied the rolloff from 0 to 0.4, and performed 10 Monte Carlo simulations for each value of the rolloff applying the proposed method. From Fig. 5 we see that the ONMSE is at acceptable levels for rolloffs below 0.3.
We also used experimentally obtained data to verify the performance of the proposed method. Since the real channels are bandlimited we applied the method proposed in Section III.B. The experimental set-up included an RF transmitter/receiver pair, in the indoor environment of the Communications and Signal Processing Laboratory of Drexel University. The equipment consisted of an Agilent ESG 4431B Vector Signal Generator (250 kHz to 6.0 GHz), an Agilent VSA 89640 Vector Signal Analyzer (dc to 2.7GHz), a VSA 89640 Analyzer software on a laptop, and two omnidirectional antennae. The test data consisted of a known sequence of 1088
8
0.25
Fig. 6. We also equalized the received signal using both channels, as seen in Fig. 7. We see that the results have close agreement.
0.2
Scatter plot
ONMSE
0.15 1
0.1 Quadrature
0.5
0.05
0
−0.5
0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
−1
0.4
roloff −1
−0.5
0 In−Phase
Fig. 5. Variation of the ONMSE with rolloff for the channels given in equation (38).
0.5
1
(a) Scatter plot
Scatter plot
2 1
1.5
1 0.5
Quadrature
0.5 Quadrature
4-QAM modulated data symbols. The Vector Signal Generator uses a Fujitsu MB86060 D/A converter chip and implements square root raised cosine pulse shaping with rolloff 0.08.
0
0
−0.5 −0.5 −1
channel obtained using training estimated channel
−1.5
−1
0.8 −2 −2
Real Imaginary
h(n)
0.6
−1.5
−1
−0.5
0 In−Phase
0.5
1
1.5
2
−1
(b)
0.4
−0.5
0 In−Phase
0.5
1
(c)
Fig. 7. Equalized data: (a) Received constellation (b) Equalized constellation using channel obtained via training and (c) Equalized constellation using the estimated channel. 400 symbols were used in the estimation and 2000 symbols were used in the equalization.
0.2
0
−0.2 2
4
6
8
10
12
14
n
Fig. 6. Channel obtained using experimental data. Lines with circles gives the channel obtained by training. Lines with stars gives the mean of the estimated channel using 400 symbols. Solid lines give the real part while broken lines give the imaginary part of the channel.
The data was transmitted at a frequency of 2.4 GHz and at a data-rate of 12 Msps. This sequence was transmitted repeatedly in a continuous mode. The receiver local oscillator was synchronized with the transmitter local oscillator by connecting the Reference clock-out of the transmitter to the reference clock-in of the receiver. To evaluate the obtained channel estimate we compared against the channel estimated based on a training approach. In the training approach we estimated the channel by cross correlating the known input sequence with the received signal. When implementing the proposed methods we used N = 196 and p = 4. The obtained results are give in
V. C ONCLUSIONS We presented a method for blindly identifying a linear time-invariant system excited by cyclostationary input, within a scalar ambiguity and an unknown delay. The basic information used is the appropriately discretized modulo2π phase of the cyclic spectrum of the system output. No system length information is needed, and no phase unwrapping is required. We proposed a general method when the channel has full bandwidth and, considering the inherent inability of frequency domain methods to work when the channel is bandlimited, proposed an approximate method for bandlimited channels. Simulation results show that they offer competitive cost-performance tradeoffs compared to more complex alternatives. VI. ACKNOWLEDGEMENT The authors would like to thank the editor and the anonymous reviewers for careful review and valuable comments.
9
R EFERENCES [1] W. A. Gardner, “Exploitation of spectral redundancy in cyclostationary signals,” IEEE Signal Processing Mag., vol. 8, pp. 14-36, Apr. 1991. [2] L. Tong, G. Xu, and T. Kailath, “Blind identification and equalization based on second order statistics: A time domain approach,” IEEE Trans. Inform. Theory, vol. 40, no. 2, pp. 340-349, Mar. 1994. [3] L. Tong, G. Xu, B. Hassibi, and T. Kailath, “Blind channel identification based on second order statistics: A frequency domain approach,” IEEE Trans. Inform. Theory, vol. 41, no. 1, pp. 329-334, Jan. 1995. [4] Y. Chen and C. L. Nikias, “Blind identification of a band-limited noninimum phase system from its output autocorrelation,” in Proc. IEEE Int. Conf. Acoust. Speech and Signal Processing, 1993, vol. 04, pp. 444-447. [5] Z. Ding, “Blind channel identification and equalization using spectral correlation measurements, part I: Frequency-domain analysis,” in Cyclostationarity in Communications and Signal Processing (W.A. Gardner, Ed.). Piscataway, NJ: IEEE, 1994, pp. 417-436. [6] Y. Li and Z. Ding, “ARMA system identification based on secondorder cyclostationarity,” IEEE Trans. Signal Processing, vol. 42, pp. 3483-3494, Dec. 1994. [7] E. Moulines, P. Duhamel, J. Cardoso and S. Mayrargue, “Subspace methods for the blind identification of multichannel FIR filters,” IEEE Trans. Signal Processing, vol. 43, no. 2, pp. 516-525, Feb. 1995. [8] H. Liu, G. Xu, L. Tong, and T. Kailath, “Recent developments in blind channel equalization: from cyclostationarity to subspaces,” Signal Processing, vol. 50, pp. 83-99, Apr 1996. [9] Z. Ding, “Characteristics of bandlimited channels unidentifiable from second-order statistics,” IEEE Signal Processing Let., vol. 3, no. 5, pp. 150-152, May 1996. [10] I. Kang, M. P. Fitz, and S.B. Gelfand, “Blind estimation of multipath channel parameters: A model analysis approach,” IEEE Trans. Commun., vol. 47, no. 8, pp. 1140-1150, Aug. 1999. [11] A.P. Liavas, P.A. Regalia, J.P. Delmas, “Blind channel approximation: effective channel order determination,” IEEE Trans. Signal Processing, vol. 47, no. 12, pp. 3336-3344, Dec. 1999 [12] L. Perros-Meilhac, E. Moulines, P. Chevalier, and P. Duhamel, “A parametric subspace based blind estimation of a SIMO-FIR with unknown channel order,” in Proc. 2’nd IEEE SPAWC, pp. 227-230, 1999 [13] B. Chen and A. P. Petropulu, “Frequency domain blind MIMO system identification based on second- and higher order statistics,” IEEE Trans. Signal Processing, vol. 49, pp. 1677-1688, Aug. 2001. [14] H. Gazzah and K. Abed-Meraim, “Blind equalization with controlled delay robust to order over estimation,” in Proc. Sixth Int. Symposium on Signal Proc. and its Applications, vol. 1, pp. 268-271, 2001. [15] J.G. Proakis, “Digital Communications,” 4th ed., McGraw Hill Inc. New York, 2000.
Sarod Yatawatta Sarod Yatawatta obtained his B.Sc. in Electrical and Electronic Engineering from the University of Peradeniya, Sri Lanka in August 2000. Since September, 2001, he is with Drexel University, Philadelphia, pursuing his Ph.D.
Athina P. Petropulu Athina P. Petropulu received the Diploma in Electrical Engineering from the National Technical University of Athens, Greece in 1986, the M.Sc. degree in PLACE Electrical and Computer Engineering in 1988 PHOTO and the Ph.D. degree in Electrical and Computer HERE Engineering in 1991, both from Northeastern University, Boston, MA. In 1992, she joined the Department of Electrical and Computer Engineering at Drexel University where she is now a Professor. During the academic year 1999-2000 she was an Associate Professor at Universit Paris Sud, cole Suprieure d’Electrcit in France. Dr. Petropulu’s research interests span the area of statistical signal processing, communications, higher-order statistics, fractional-order statistics and ultrasound imaging. She is the co-author of the textbook entitled, Higher-Order Spectra Analysis: A Nonlinear Signal Processing Framework, (Englewood Cliffs, NJ: Prentice-Hall, Inc., 1993). She is the recipient of the 1995 Presidential Faculty Fellow Award. She has served as an associate editor for the IEEE Transactions on Signal Processing, IEEE Signal Processing Letters and the EURASIP Journal on Wireless Communications and Networking. She is a member-at-Large of the IEEE Signal Processing Board of Governors, a member of the IEEE Signal Processing Society Conference Board and the IEEE Technical Committee on Signal Processing Theory and Methods. She is the General Chair of the 2005 International Conference on Acoustics Speech and Signal Processing (ICASSP-05), to be help in Philadelphia PA.
Riddhi Dattani Riddhi Dattani received her M.S. degree in Electrical/Telecomm Engg. from Drexel University, Philadelphia, Pennsylvania, USA in 2003 and Bachelor of Engg.(B.E.) degree in Electronics Engg. from Mumbai University, Mumbai, India in 2001. Her research interests include applications of signal processing to wireless communications and RF circuit designing.
Blind Channel Estimation Using Fractional Sampling Sarod Yatawatta, Student Member, IEEE, Athina P. Petropulu, Senior Member, IEEE, and Riddhi Dattani, Student Member, IEEE,
Abstract— We consider the problem of blind estimation of a communication channel by oversampling the channel output. We propose a non-parametric approach which, based on the cyclic spectrum, finds the channel phase response without the need of phase unwrapping nor channel length information. For bandlimited channels, the cyclic spectrum is not supported for the full bandwidth. For this case, we propose an approximation for the discretized phase of the cyclic spectrum, which, under certain conditions, results in a simpler channel estimation method. The proposed approach is applied to simulated data and real recordings, and is compared to existing methods. Index Terms— Blind channel estimation, cyclostationary inputs, cyclic spectrum.
I. I NTRODUCTION
T
HE communications channel sets a limit to the transmission rate of the communication system. Signal pulses are broadened in time as they travel through the channel (multipath effect), leading to intersymbol interference (ISI) at the receiver. In order to remove the effect of ISI, prior knowledge of the channel is required. Generally this information is obtained by transmitting training symbols. However, this approach implies waste of bandwidth, which can be quite significant in cases where the channel is time varying. An alternative approach is blind channel estimation and equalization. Within short time intervals, estimation of the communication channel based on the received data can be viewed as a blind linear time invariant (LTI) system estimation problem. In general, for the estimation of an LTI system excited by an unknown stationary non-Gaussian input, based on the system output, one has to use higher-order statistics (HOS) of the output sampled at baud rate. When the input is cyclostationary, which is the case of communications signals, system estimation can be carried out using second-order cyclic statistics of the oversampled output [1], [2]. Based on this premise, several methods have been proposed for blind channel estimation (for a review of existing results the reader is referred to [8]). This work was supported by ONR under grant ONR-N00014-20-1-0137 The authors are with The Department of Electrical & Computer Engineering, Drexel University, Philadelphia, PA 19104
They can be categorized as time-domain and frequencydomain, with reference to the transform domain in which the problem is formulated. A serious deficiency in time domain-methods ([2], [7]) is their sensitivity to channel order estimation errors. Most often in practice the channel (the combined effect of transmit-receive filters and multipaths) is a continuous function of time and its magnitude will not suddenly drop to zero at one point. The notion of channel length depends on the application. Some methods have been proposed to overcome the need for channel length information [11], [14], [12], yet this could still pose a problem when the channel is time varying. On the other hand, frequency-domain methods do not require channel length, but cannot be applied to band-limited channels. In that case, the phase of the full spectrum cannot be recovered because the cyclic spectrum has limited support [5]. In practice, most channels are bandlimited, both by the physical nature of the channel and by the pulse shaping employed. In this paper we propose a novel frequency-domain blind channel estimation approach. The channel phase frequency response is obtained in closed form in terms of the appropriately discretized cyclic phase spectrum. For the case of bandlimited channels with a small rolloff, we propose an approximate expression for the channel phase response, which involves low computational complexity. The resulting channel estimate will always contain an error, however the error will be negligible if the channel has negligible energy within the stopband. Despite this bias, the proposed approach can still be of great practical use since it involves low cost. As a matter of practice, choosing the correct algorithm for a particular purpose does not depend on the accuracy of the result only, but also on economy in hardware implementation. With the introduction of wireless communication capabilities to most electronic devices in the offing, the proposed approach could provide a viable alternative in reducing cost. The paper is organized as follows. Section II provides the problem formulation and summarizes relevant literature. Section III details the proposed approach. In section III.A, the general method is described and in section III.B, the simplification applicable to bandlimited channels is given.
2
Section IV provides results on the performance of the proposed methods based simulations and also based on real data recordings. Finally, concluding remarks are given in Section V. II. P ROBLEM F ORMULATION Let us consider the baseband representation of a digital data communication system. Assuming that the channel is LTI and BIBO (Bounded Input producing Bounded Output) stable, the received signal is given by: x(t) =
∞ X
sν h(t − νT ) + w(t), ν ∈ Z
(1)
ν=−∞
where sν is the transmitted signal, which is assumed to be zero-mean i.i.d. with unit variance; h(t) is the channel impulse response; T is the symbol interval, and w(t) is stationary, zero-mean, white noise with variance No and is independent of sν . The signal x(t) is cyclostationary with cyclic period T [1]. Sampling x(t) with sampling interval ∆ = Tp , p ∈ Z, yields: 4
x(n) = x(n∆) =
∞ X
sν h(n − νp) + w(n), n ∈ Z (2)
ν=−∞
where h(n) = h(n∆), w(n) = w(n∆). The discrete time signal x(n) is a cyclostationary process with cyclic period p. Cyclostationarity implies that the autocorrelation function is periodic with period p. The Fourier series coefficients of the autocorrelation are referred to as cyclic autocorrelation, Rxkα (ν), and the Fourier Transform of those as the cyclic spectrum, Sxkα (ejω ) [1], i.e., Rxkα (ν) =
p−1 X
Rx (n + ν, n)e−j2πkαn , α =
n=0
1 , p
(3)
where ψ(ω), φ(ω) are the phases of Sxkα (ω) and H(ejω ), respectively. It has been shown in [4], that the above equation cannot recover the phase response of an arbitrary channel, since when evaluated for different ω’s it always leads to an underdetermined system of equations. Viewing ψ(ω), φ(ω) as Fourier Transforms of sequences Ψ(τ ) and Φ(τ ), respectively, it holds [5] that: Φ(τ ) =
jΨ(τ ) , τ 6= νT, ν ∈ Z. 2sin(πτ /T )
(7)
Since Φ(τ ) is infinite at τ = νT , the phase φ(ω) cannot be obtained from (7). However, if Φ(τ ) is bounded at τ = νT , the missing information can be recovered via interpolation, assuming continuity of Φ(τ ). The case where Φ(τ ) is not bounded, i.e., contains impulses at τ = νT occurs when φ(ω) contains sinusoidal components with frequency νT . This observation led to the identifiability condition [3] stating that discrete channels with zeros uniformly 2π/p spaced on a circle are not identifiable because they result in oscillatory phase. However, even if the channel is identifiable, (7) will be inaccurate for τ ≈ νT . Parametric approaches that obtain the channel parameters from cepstral information obtained for τ away from the points νT were proposed in [5]. III. P ROPOSED M ETHOD We first consider an FIR channel that has a full passband. Later, in section B, we consider the situation where the channel is bandlimited. A. FIR Channels
First we make the following assumption: (A1) The channel satisfies identifiability conditions given in n, ν, k ∈ Z = [3] and [9], i.e., there are no uniformly 2π/p spaced ν=−∞ zeros on a circle and there are no spectral nulls in the where Rx (ν, n) denotes autocorrelation of x(n). non-rolloff passband. It holds that: We next show that an appropriately discretized version kα jω Sx (e ) = H(ejω )H ∗ (ej(ω−2πkα) ) + pNo δ(k), k ∈ Z. of (6) can lead to a solution of the channel phase. (4) Let ω = 2π N l, l = 0, ..., N − 1. Then, This is a key relationship that can lead to the estimation of ψ(l) = φ(l + m) − φ(l − m) (8) the channel frequency response H(ejω ). The magnitude response of the channel can be easily or, recovered from the power-spectrum of x(n). To obtain ψ(l + m) = φ(l + 2m) − φ(l) (9) phase information, we can rewrite (4) as: Sxkα (ejω )
4
Sxkα (ω) = =
∞ X
Rxkα (ν)e−jνω ,
Sxkα (ej(ω+πkα) ), (5) H(ej(ω+πkα) )H ∗ (ej(ω−πkα) ), k 6= 0
from which we get: ψ(ω) = φ(ω + πkα) − φ(ω − πkα), k = ±1, ±2, .... (6)
where m is the closest integer to N kα/2 = N k/2p. Let us fix k to the value k = 1 for maximum support [1]. Based on (9) evaluated for l = 1, ..., N − 1, and taking φ(0) = 0, we can formulate the system of equations: AN,m Φ = Ψ
(10)
3
where Φ = [φ(1) ... φ(N −1)]T , Ψ = [ψ(m+1) ... ψ(m+ (A2) The transition band of the channel is negligible compared to the passband, i.e., the rolloff is very small. N − 1)]T , and ψ(l) = ψ(l + N ). AN,m is a (N − 1) × (N − 1) sparse matrix. It can be (A3) The energy within the stopband of the channel is negligible. viewed as consisting of the columns 2 to N and rows 2 to N of an N × N circulant matrix, B, whose first row (A4) Let Lc = bωc N/2πc (cutoff frequency index), where ωc is the cutoff frequency, then equals: b = [−1, 0T2m−1 , 1, 0TN −2m−1 ]. (11) Lc < N/4. (14) Here 0Tk denotes a row vector containing k zeros. The other (A5) The number of zeros of the channel on the unit circle Nz satisfies, rows of B are circularly shifted versions of b. N − 2Lc A matrix similar to AN,m was also encountered in [13] Nz > L . (15) N in a different context, where it was shown that for N and where L is the channel length. 2m co-prime numbers, it is full rank with determinant equal Assumption (A4) implies that ωc < π/2 rad/cycle. For to one. Since N is a number that we select, we can give it a any system that transmits at or below Nyquist rate [15], the value so that it is co-prime with 2m = 2bN/2pc, thus channel sampled at symbol rate (p = 1) has full bandwidth. guaranteeing full rank for A. The phase vector Φ can then Hence to satisfy this condition we need p ≥ 3. In general, be obtained by solving (10). Ideally, we would like the we have no prior knowledge of Nz to check the validity of phase expression to involve the modulo 2π phases only, so assumption (A5). However, by careful design of the shaping that no phase unwrapping is necessary. It turns out that function, we can ensure that the minimum number of Nz Φ can actually be based on modulo 2π phases, and this is is greater than the one required by (15). Multipaths will due to the special properties of matrix AN,m [13]. It can only add more zeros thus the condition will hold under multipaths if the shaping function itself satisfies it, unless be shown that, the estimate: the channel becomes unidentifiable according to [9]. 4 ˜ b = Let us discretize ψ(ω) in (6) on a grid with spacing (12) Φ A−1 N,m Ψ 2π/N and select N and p so that whenever φ(ω) falls in ˜ is the modulo 2π phase corresponding to Ψ, the passband, φ(ω + 2πkα) falls in the stopband. We take where Ψ satisfies: k = 1 again for maximum support. This can be achieved 2π T b = Φ+φ(0)1(N −1)×1 + ν[1, ..., N −1] +2πK (13) as long as, Φ Lc < m < N/2 − Lc (16) N where K is a vector of size (N − 1) × 1 with integer elements, and ν is an unknown integer. Combining the phase estimate of (12) with the correct magnitude spectrum, results in the channel impulse response with an unknown scalar multiplicative constant, and a circular time shift equal to ν. We should note here that in comparison with the methods of [5],[6], the proposed method does not require phase unwrapping. B. FIR Channels with pulse shaping In this section, we consider the case where the channel is bandlimited by pulse shaping. The objective of pulse shaping is to restrict the spectrum of the signal within a fixed frequency interval [15]. Hence, by its design, the magnitude within the stopband is low compared to the passband. Moreover, (4) will have limited support, thus (9) will result in an underdetermined system of equations. In [5], it has been suggested to use parametric methods to overcome this difficulty. Here we propose a non parametric approach that yields an approximate channel estimate. In addition to assumption (A1) we make the following assumptions:
where m is the closest integer to αN/2 or N/2p. A value for m greater than 1 is guaranteed by N > 4Lc . Proposition 1: Let h(n) be an FIR impulse response of d length L, which satisfies assumptions (A1)-(A5). Let h(n) be the channel impulse response constructed based on the correct magnitude response of h(n), i.e., |H(l)|, and phase d constructed using the estimated phase of the response φ(l), \ cyclic spectrum, ψ(l + m), as: \ −ψ(l + m); 4 d = φ(l) 0 ≤ l ≤ Lc − 1, N − Lc ≤ l ≤ N − 1 0; Lc ≤ l ≤ N − Lc − 1 (17) with N and p satisfying (16). Then, it holds: d h(n) (18) d
=
4
IDF T {|H(l)|ej φ(l) }
=
ah(n + b) ⊗
2π
−2N e−j N ko n 2π
Λ(1 − ej N n ) ! Λ X (2i + 1)N δ n− × + E(n) 2Λ i=0
4
where E(n) =
a N
N −L Xc −1
2π
H(l)ej N ln
− H(l)e
j{ 2π N (b+n)l+π
P(N −1)/Λ−1 i=Nz
a = e−jNz (N/2−m) , b = Λ=b
(19)
l=Lc
N − 2Lc c, Nz − 1
u(l+ko −iΛ)}
!
,
Nz , ko = 2m − Lc , 2
NX z −1 ω π − ωc φ(ω) = Nz (π− )+φt (ω)−π ) u(ω−ωc −i 2 2 Nz − 1 i=0 (23) where u(ω) is the step function. Using (23), we can rewrite φ(ω + 2πα) in discrete form (ω = 2π N l) as:
and ⊗ denotes N point circular convolution. Proof: The channel Frequency Response can be given as: H(ω) =
NY total
each zero on the unit circle. Note that these jumps will not occur within the passband of the channel. Hence, if Nz is large enough, within the stopband, φz (ω) will dominate the overall phase response because all its derivatives become large whenever ω ≈ θi for any i ∈ [0, Nz − 1]. If the phase due to zeros away from the unit circle is φt (ω), the phase of the channel can be given as
φ(l + 2m) (1 − zi z −1 )|z=ejω
(24)
Nz π l + φt (l + 2m) = Nz (N/2 − m) − N NX z −1 −π u(l + 2m − Lc − iΛ)
(20)
i=0
Here the channel has Ntotal zeros, out of which Nz lie on the unit circle.
i=0
Zi r
where l is the discrete frequency index. Under the assumptions made earlier, φ(l) will lie on the passband while φ(l + 2m) will lie on the stopband. Hence the jumps of π will dominate (24) thus we can ignore the contribution from φt (l + 2m) to simplify our analysis. Hence we can write (9) as
φi di ωc − ω r≈1
−ψ(l + m)
(25)
Nz π ≈ φ(l) + l − Nz (N/2 − m) N NX z −1 +π u(l + 2m − Lc − iΛ).
Fig. 1. Calculation of Phase Due to Channel Zeros Corresponding to the Stopband. As ω increases, the phase contribution from these zeros change as if the location of the zeros rotate clockwise.
i=0
Let us refer to Fig. 1 to calculate the phase of the terms of (20) that contain the zeros on the unit circle. As ω increases, the phase contributed by this set of zeros changes as if the zeros themselves are rotated by ω clockwise. Let φi be the phase of (1 − zi e−jω ) where |zi | = 1. It holds: φi (ω) =
π ωc − ω π − ωc + +i , 0 ≤ i < Nz . 2 2 Nz − 1
(21)
The sum of the phases of all terms (1−zi e−jω ) for which |zi | = 1, i.e. φz (ω) will then be equal to Nz (π − ω2 ). For ω ≈ θi , it holds: r 1 dφi (ω) = ((1−r)2 − (1+r)(θi −ω)2 ), ω ≈ θi dω (1 − r)3 2 (22) i (ω) We see that as r → 1, dφdω becomes infinite, thus indicating discontinuity of φi at ω = θi . We see from Fig. 1 that at this discontinuity, the overall change in φi (ω) is π. Hence, the overall phase response is linear with jumps of π at frequencies corresponding to
Then we can re write (17) as: φ(l) + NNz π l − Nz (N/2 − m) PNz −1 + π i=0 u(l + 2m − Lc − iΛ); d ≈ φ(l) 0 ≤ l ≤ L c − 1, N − Lc ≤ l ≤ N − 1 0; Lc ≤ l ≤ N − Lc − 1. (26) d d = IDF T {|H(l)|ej φ(l) Then h(n) } yields:
d = h(n) a N
+
LX c −1
(27)
2π
H(l)ej{ N (b+n)l+π
l=0 N −L Xc −1
PNz −1
2π
H(l)ej{ N (b+n)l+π
i=0
u(l+ko −iΛ)}
P(N −1)/Λ−1 i=Nz
u(l+ko −iΛ)}
l=Lc
+
N −1 X
l=N −Lc
+ E(n)
2π
H(l)ej{ N (b+n)l+π
PNz −1 i=0
u(l−N +1+ko −iΛ)}
5
=
N −1 2π a X H(l)ej{ N (b+n)l} G(l) + E(n) N l=0
with G(l) = 2
M X
rect(
i=0
l + ko − i2Λ ) − 1, M ∈ Z Λ
(28)
where E(n) can be deduced from (27) and rect( Λτ ) is a rectangular pulse of width Λ. We see that the IDFT of G(l) is 2π
g(n) =
2π
e−j N ko n (1 − e−j N Λn ) (1 − ej2πn ) (1 + e−j
2π N Λn
2π
(1 − ej N n )
)
(29)
N N where g(n) 6= 0 for n = 2Λ , 3 2Λ , . . .. By simplification of (29), (18) follows. (E.O.P) d will contain delayed Equation (18) suggests that h(n) copies of the original channel separated be N Λ . If condition (15) holds, these copies will be separated well enough to allow one to extract h(n). Note also that the error E(n) will be small (by assumption (A3)) if the energy within the stopband is negligible, i.e., N −L −1 2 Xc |E(n)| ≤ |H(l)|. (30) N l=Lc
The magnitude of each replica appearing in (27) is determined by the magnitudes of the impulses in (18), which in turn is determined by f (n) =
1 2π
|1 − ej N n |
.
(31)
This expression has its maxima at n = 0 and n = N , is monotonically decreasing as n goes from 0 to N/2 and, is monotonically increasing as n goes from N/2 to N . Hence the first and last replicas will have the highest magnitudes. As a consequence, to extract one copy of the channel d we could do the following. Let the estimate out of h(n) anticipated channel length be Le . We find the segment with maximum energy out of the first Le samples and the last Le samples. After extracting the maximum energy sequence, we obtain a channel estimate with a scaling ambiguity and an unknown delay. IV. S IMULATION R ESULTS Simulations were done to verify the performance of the proposed methods and to compare with other methods. In section A, we give the results obtained for FIR channels with full bandwidth support, and in section B we give results for bandlimited channels. Finally, in section C, we give results of the application of the method in section III.B to real measurements.
The basic procedure is P as follows. First, we generated ∞ a QAM signal x(n) = ν=−∞ sν h(n − νp) + w(n). The input symbols {sν } were from an i.i.d. QAM source. The noise w(n) was stationary, white and Gaussian. The Spectral Correlation Density (SCD) or cyclic spectrum was estimated as follows. First we estimated the unbiased autocorrelation as: N −1 X bx (n + ν, n) = 1 R x(n + ν + pl)x∗ (n + pl), n, ν ∈ Z. N l=0 (32) bxkα (ν) Next, we estimated the cyclic autocorrelation R using equation (4). If some estimate of the channel length, ˆ is available, then we can window R bxkα (ν) by a i.e., L ˆ rectangular window of length 2L. Finally, we estimated the SCD by taking an N point FFT, i.e.,
bxkα (ν)), α = 1 , k = 1. Sbxkα (l) = FFT(R p
(33)
b The phase of the channel was estimated using ψ(l), the modulo 2π phase of Sbxα (l), according to methods proposed in sections III.A and III.B. In order to make b ψ(l) symmetric, it was shifted circularly by b−N/2pc. (0) The power spectral density, Sbx (l) was used for amplitude estimation. b b = |Sbx(0) (l)|1/2 exp(j φ(l)) H(l)
(34)
Finally, the impulse response was computed as: b b h(n) = IFFT(H(l)).
(35)
The performance was measured using and the Overall Normalized Mean Square Error (ONMSE) and also Bit Error Rate (BER). We tested a total of Mh different channels. For each channel we performed Mc Monte Carlo runs corresponding to independent input realizations. The ONMSE was defined as: Mh Mc P L j 2 bj 1 X 1 X n=0 |hi (n) − h (n)| ONMSE = PL j 2 Mh j=1 Mc i=1 n=0 |h (n)| (36) ˆ j is the estimate of the i-th Monte Carlo run of the where h i j-th channel. Based on each channel estimate we implemented a zero forcing equalizer and obtained the input estimate. The BER rate was computed for each channel and then averaged over the different channels tested. A. FIR Channels Each channel was taken to be FIR of length 8, and its coefficients were generated randomly. For each channel we performed Mc = 25 Monte Carlo runs, using the estimation method of Section III.A. Finally the ONMSE and BER
6
B. FIR Channels with pulse shaping We generated channels of the form h(t) = a1 c(t + 1.5T, 0.11) + a2 c(t − T1 , 0.11) (37) + a3 c(t − 1.5T, 0.11) W6T (t)
where a1 , a2 , a3 are normally distributed random numbers with mean zero and variance one and T1 is uniformly distributed in [−T, T ]. The shaping function c(t − τ, β) is a
−0.3
10
with windowing without windowing
−0.4
ONMSE
10
−0.5
10
Proposed Cepstrum based
−0.6
10
−0.7
10
0
5
10
15
20
25
30
35
SNR/(dB)
(a) 0
10
with windowing without windowing Proposed Cepstrum based
BER
were obtained as the average over Mh = 100 different channels. For comparison purposes, we also implemented the nonparametric method (cepstrum based) of [6], where the phase of the channel is estimated using the Fourier Series expansion of the unwrapped phase of the SCD. bxα (ν) was obtained with a rectangular The estimate of R windowing of length 2L as well as without windowing. The obtained estimate, usually contains more than one replicas of the channel response. Although the precise reason for such behavior needs to be further investigated, simulations suggest that this is a result of phase errors and the way they propagate to the final solution. For large N , i.e., N > 3L, there is enough room for the replicas to be separated. One could then extract the channel by identifying the segment with maximum energy, with length equal to the anticipated channel length, or length equal to the window size used in the equalizer. There will be at least one such disjoint segment in the estimate where the energy will be distinctly greater than the rest of the estimate. In the proposed method we chose N = 131 and p = 4. In the simulations we extracted the segment with maximum energy of length L, and we have given results in Fig. 2. ONMSE results are shown in Fig.2(a). From Fig. 2(a), bxα (ν), was not we see that the windowing used to estimate R that critical in the final estimate. Since the window requires some knowledge of channel length information, the above observation confirms that the methods tested here are not very sensitive to the channel length information. Corresponding BER results are shown in Fig. 2(b). Here we used channel length knowledge in order to compute the Zero Forcing equalizer. The results were obtained for 1000 symbols. The above results suggest that the proposed approach is only slightly better than that of [6] in terms of ONMSE and has same performance in BER. However, in terms of complexity, the proposed approach is less intensive than that of [6], since it requires no phase unwrapping. The matrix A is independent of the data, hence we only need its inverse and moreover, the inverse only has ones and zeros. Hence solving equation (12) does not involve multiplications. In contrast, the cepstrum based method [6] requires phase unwrapping and two additional FFTs.
−1
10
−2
10
0
5
10
15
20
25
30
35
SNR/(dB)
(b) Fig. 2. (a) ONMSE and (b) BER for cepstrum based estimates and proposed estimates given in section III.A. We generated 100 random FIR channels of length 8 and for each channel, obtained the estimates for 25 Monte Carlo runs. 1000 symbols were used in channel estimation and perfect channel length knowledge was assumed. In (a), we have given the ONMSE results where both a rectangular window of length 2L was used, and a window was not used.
raised cosine [15] with delay τ and rolloff β, and W6T (t) is a rectangular window of width 6 symbol intervals T . We here considered Mc = 25, Mh = 100, N = 128 and p = 4. We compared the method proposed in section III.B with the high complexity subspace based method of [7]. In Fig. 3 we have given ONMSE and BER results obtained with perfect channel length knowledge, and in Fig. 4 we have given results where the channel length was overestimated by 1 symbol (4 samples). When perfect channel length knowledge is available, the method of [7] performs exceptionally well for high SNR, while the proposed method performs slightly better for low SNR. On the other hand, the proposed approximate method outperforms the former when there is a slight mismatch in channel length estimation.
7
0
0
10
10
Proposed Subspace
Proposed Subspace
−1
ONMSE
ONMSE
10
400 Symbols 1000 Symblos
−2
10
400 Symbols 1000 Symblos
−3
10
0
−1
10
−2
5
10
15
20
25
30
10
35
0
5
10
SNR/(dB)
15
20
25
30
35
SNR/(dB)
(a)
(a)
0
0
10
10
Proposed Subspace
Proposed Subspace −1
BER
BER
10
−2
10
−1
10
400 Symbols 1000 Symblos
400 Symbols 1000 Symblos
−3
10
−4
10
0
−2
5
10
15
20
25
30
35
SNR/(dB)
10
0
5
10
15
20
25
30
35
SNR/(dB)
(b)
(b)
Fig. 3. (a) ONMSE and (b) BER for the method proposed in section III.B and for the subspace method. We generated 100 channels with pulse shaping and for each channel, 25 Monte Carlo runs were taken at each SNR. The modulation scheme used was 4 QAM and perfect channel length knowledge was used for equalization.
Fig. 4. (a) ONMSE and (b) BER for the method proposed in section III.B and for the subspace method. We generated 100 channels with pulse shaping and for each channel, 25 Monte Carlo runs were taken at each SNR. The modulation scheme used was 4 QAM and the channel length was overestimated by one symbol interval.
C. Performance with experimental data To check whether the channels given in (38) satisfy (15) we can do the following. If we oversample by 4 (p = 4, T = 1), we see that L = 24. Moreover, Lc ≈ 20, and for N = 128 (15) gives Nz > 16.5. By evaluating the zeros of the channel we find Nz = 18, indeed satisfying (15). Finally, we investigated the effect of rolloff to the error in approximation. We generated 30 channels of the form (38) with variable rolloff. We varied the rolloff from 0 to 0.4, and performed 10 Monte Carlo simulations for each value of the rolloff applying the proposed method. From Fig. 5 we see that the ONMSE is at acceptable levels for rolloffs below 0.3.
We also used experimentally obtained data to verify the performance of the proposed method. Since the real channels are bandlimited we applied the method proposed in Section III.B. The experimental set-up included an RF transmitter/receiver pair, in the indoor environment of the Communications and Signal Processing Laboratory of Drexel University. The equipment consisted of an Agilent ESG 4431B Vector Signal Generator (250 kHz to 6.0 GHz), an Agilent VSA 89640 Vector Signal Analyzer (dc to 2.7GHz), a VSA 89640 Analyzer software on a laptop, and two omnidirectional antennae. The test data consisted of a known sequence of 1088
8
0.25
Fig. 6. We also equalized the received signal using both channels, as seen in Fig. 7. We see that the results have close agreement.
0.2
Scatter plot
ONMSE
0.15 1
0.1 Quadrature
0.5
0.05
0
−0.5
0 0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
−1
0.4
roloff −1
−0.5
0 In−Phase
Fig. 5. Variation of the ONMSE with rolloff for the channels given in equation (38).
0.5
1
(a) Scatter plot
Scatter plot
2 1
1.5
1 0.5
Quadrature
0.5 Quadrature
4-QAM modulated data symbols. The Vector Signal Generator uses a Fujitsu MB86060 D/A converter chip and implements square root raised cosine pulse shaping with rolloff 0.08.
0
0
−0.5 −0.5 −1
channel obtained using training estimated channel
−1.5
−1
0.8 −2 −2
Real Imaginary
h(n)
0.6
−1.5
−1
−0.5
0 In−Phase
0.5
1
1.5
2
−1
(b)
0.4
−0.5
0 In−Phase
0.5
1
(c)
Fig. 7. Equalized data: (a) Received constellation (b) Equalized constellation using channel obtained via training and (c) Equalized constellation using the estimated channel. 400 symbols were used in the estimation and 2000 symbols were used in the equalization.
0.2
0
−0.2 2
4
6
8
10
12
14
n
Fig. 6. Channel obtained using experimental data. Lines with circles gives the channel obtained by training. Lines with stars gives the mean of the estimated channel using 400 symbols. Solid lines give the real part while broken lines give the imaginary part of the channel.
The data was transmitted at a frequency of 2.4 GHz and at a data-rate of 12 Msps. This sequence was transmitted repeatedly in a continuous mode. The receiver local oscillator was synchronized with the transmitter local oscillator by connecting the Reference clock-out of the transmitter to the reference clock-in of the receiver. To evaluate the obtained channel estimate we compared against the channel estimated based on a training approach. In the training approach we estimated the channel by cross correlating the known input sequence with the received signal. When implementing the proposed methods we used N = 196 and p = 4. The obtained results are give in
V. C ONCLUSIONS We presented a method for blindly identifying a linear time-invariant system excited by cyclostationary input, within a scalar ambiguity and an unknown delay. The basic information used is the appropriately discretized modulo2π phase of the cyclic spectrum of the system output. No system length information is needed, and no phase unwrapping is required. We proposed a general method when the channel has full bandwidth and, considering the inherent inability of frequency domain methods to work when the channel is bandlimited, proposed an approximate method for bandlimited channels. Simulation results show that they offer competitive cost-performance tradeoffs compared to more complex alternatives. VI. ACKNOWLEDGEMENT The authors would like to thank the editor and the anonymous reviewers for careful review and valuable comments.
9
R EFERENCES [1] W. A. Gardner, “Exploitation of spectral redundancy in cyclostationary signals,” IEEE Signal Processing Mag., vol. 8, pp. 14-36, Apr. 1991. [2] L. Tong, G. Xu, and T. Kailath, “Blind identification and equalization based on second order statistics: A time domain approach,” IEEE Trans. Inform. Theory, vol. 40, no. 2, pp. 340-349, Mar. 1994. [3] L. Tong, G. Xu, B. Hassibi, and T. Kailath, “Blind channel identification based on second order statistics: A frequency domain approach,” IEEE Trans. Inform. Theory, vol. 41, no. 1, pp. 329-334, Jan. 1995. [4] Y. Chen and C. L. Nikias, “Blind identification of a band-limited noninimum phase system from its output autocorrelation,” in Proc. IEEE Int. Conf. Acoust. Speech and Signal Processing, 1993, vol. 04, pp. 444-447. [5] Z. Ding, “Blind channel identification and equalization using spectral correlation measurements, part I: Frequency-domain analysis,” in Cyclostationarity in Communications and Signal Processing (W.A. Gardner, Ed.). Piscataway, NJ: IEEE, 1994, pp. 417-436. [6] Y. Li and Z. Ding, “ARMA system identification based on secondorder cyclostationarity,” IEEE Trans. Signal Processing, vol. 42, pp. 3483-3494, Dec. 1994. [7] E. Moulines, P. Duhamel, J. Cardoso and S. Mayrargue, “Subspace methods for the blind identification of multichannel FIR filters,” IEEE Trans. Signal Processing, vol. 43, no. 2, pp. 516-525, Feb. 1995. [8] H. Liu, G. Xu, L. Tong, and T. Kailath, “Recent developments in blind channel equalization: from cyclostationarity to subspaces,” Signal Processing, vol. 50, pp. 83-99, Apr 1996. [9] Z. Ding, “Characteristics of bandlimited channels unidentifiable from second-order statistics,” IEEE Signal Processing Let., vol. 3, no. 5, pp. 150-152, May 1996. [10] I. Kang, M. P. Fitz, and S.B. Gelfand, “Blind estimation of multipath channel parameters: A model analysis approach,” IEEE Trans. Commun., vol. 47, no. 8, pp. 1140-1150, Aug. 1999. [11] A.P. Liavas, P.A. Regalia, J.P. Delmas, “Blind channel approximation: effective channel order determination,” IEEE Trans. Signal Processing, vol. 47, no. 12, pp. 3336-3344, Dec. 1999 [12] L. Perros-Meilhac, E. Moulines, P. Chevalier, and P. Duhamel, “A parametric subspace based blind estimation of a SIMO-FIR with unknown channel order,” in Proc. 2’nd IEEE SPAWC, pp. 227-230, 1999 [13] B. Chen and A. P. Petropulu, “Frequency domain blind MIMO system identification based on second- and higher order statistics,” IEEE Trans. Signal Processing, vol. 49, pp. 1677-1688, Aug. 2001. [14] H. Gazzah and K. Abed-Meraim, “Blind equalization with controlled delay robust to order over estimation,” in Proc. Sixth Int. Symposium on Signal Proc. and its Applications, vol. 1, pp. 268-271, 2001. [15] J.G. Proakis, “Digital Communications,” 4th ed., McGraw Hill Inc. New York, 2000.
Sarod Yatawatta Sarod Yatawatta obtained his B.Sc. in Electrical and Electronic Engineering from the University of Peradeniya, Sri Lanka in August 2000. Since September, 2001, he is with Drexel University, Philadelphia, pursuing his Ph.D.
Athina P. Petropulu Athina P. Petropulu received the Diploma in Electrical Engineering from the National Technical University of Athens, Greece in 1986, the M.Sc. degree in PLACE Electrical and Computer Engineering in 1988 PHOTO and the Ph.D. degree in Electrical and Computer HERE Engineering in 1991, both from Northeastern University, Boston, MA. In 1992, she joined the Department of Electrical and Computer Engineering at Drexel University where she is now a Professor. During the academic year 1999-2000 she was an Associate Professor at Universit Paris Sud, cole Suprieure d’Electrcit in France. Dr. Petropulu’s research interests span the area of statistical signal processing, communications, higher-order statistics, fractional-order statistics and ultrasound imaging. She is the co-author of the textbook entitled, Higher-Order Spectra Analysis: A Nonlinear Signal Processing Framework, (Englewood Cliffs, NJ: Prentice-Hall, Inc., 1993). She is the recipient of the 1995 Presidential Faculty Fellow Award. She has served as an associate editor for the IEEE Transactions on Signal Processing, IEEE Signal Processing Letters and the EURASIP Journal on Wireless Communications and Networking. She is a member-at-Large of the IEEE Signal Processing Board of Governors, a member of the IEEE Signal Processing Society Conference Board and the IEEE Technical Committee on Signal Processing Theory and Methods. She is the General Chair of the 2005 International Conference on Acoustics Speech and Signal Processing (ICASSP-05), to be help in Philadelphia PA.
Riddhi Dattani Riddhi Dattani received her M.S. degree in Electrical/Telecomm Engg. from Drexel University, Philadelphia, Pennsylvania, USA in 2003 and Bachelor of Engg.(B.E.) degree in Electronics Engg. from Mumbai University, Mumbai, India in 2001. Her research interests include applications of signal processing to wireless communications and RF circuit designing.
