605

IEEE TRANSACTIONS ON COMMUNICATIONS, VOL 36, NO. 5 , MAY 1988

Digital Filter and Square Timing Recovery MARTIN OERDER,

STUDENT MEMBER, IEEE, AND

Abstruet-Digitalrealizationsof timing recovery circuits for digital data transmission are of growing interest.In this paper, we present a new digital algorithm which can be implemented very efficiently also at high data rates. The resulting timing jitter has been computed and verified by simulations. In contrast to other known algorithms, the one presented here allows free running sampling oscillators and a new planar filtering method which prevents synchronization hangups.

HEINRICH MEYR,

FELLOW, IEEE

11. TIMING ESTIMATION Here we considerthetimingrecoveryfordigitaldata transmission by linearmodulationschemes(PAM,QAM, PSK). The rcccived signal (PAM) or the equivalent low-pass signal (QAM, PSK) can bc written as m

r(t)=

angT(t-nT-c(t)T)+n(t) n=-w

I . INTRODUCTION

=U ( t ) +

n(t).

D.

(1)

IGITALrealizationsofreccivcrsforsynchronousdata Where a, are the complex valucd transmitted symbols with slgnals-baseband as well as QPSK or QAM signals-are of growing interest as the capabilitiesof signal processors (for mean power 1 (e.g., + 1, f i with QPSK), g T ( t ) is the low data rates) and application specific intcgratcd circuits (for transmission signal pulse, T i s the symbol duration, n(r) is the high data rates)increase.Thesereceivers have to include channel noise which is assumed to be white and Gaussian with algorithmsfortimingrecovery.Several such discrete-time power density NO.and t ( t ) is an unknown,slowlyvarying time delay. algorithms have been proposed during the last few years [ 11Now timing recovery means the estimation of thc delay E ( t ) [ 3 ] .The majority of these solutions, howcvcr. include only the integration of one part of the timing synchronization, namely, to enable the optimal dctcction of the data. Because E varies in a digital realization, we canprocessthe the generation of some kind of timing error signal, into the veryslowly, digital part of the receiver. This error signal is then typically rcccived signal section by section. For each section A m , we used to control an analog VCO which generates the sampling can assume E to be constant and obtain an estimate ?m. This estimate must then be combined with the previous estimates strobes. Due to the advantages of an integrated realization, however, (i.e., it must be filtered) such that the optimal estimate Cm is obtained. The latter canbe used to control an analog or digital as much of the receiver as possiblc should bc digital. This means that the input signal should be sampled at fixed a rate by sampler for the detection. Below we consider a special type of timing estimator which a free running oscillator and all further processing should then is particularly suited for digital realization. It is similar to the be done digitally using these samples. For symbol detection, in the input this mcans that thc optimum decision metrics must be continuous-timc filter and square synchronizer that generated from the given samples by some sortof interpolation signal is squared and the resulting spectral component at the which is controlled by an estimate of the current timing offset symbol rate is extracted by a filtering operation. In Fig. 1 the algorithm is shown. After a receiving filter [impulse response [4]. Thcreforc, wc need an algorithm which determines this ) sampled at rate N / T absolute timing offset(not only a timing error signal) from the gR(t)1 the signal f ( r ) = r ( t ) * g ~ ( 1 is ("*" denotes a convolution). We thus have samples given samples of the signal. Such an algorithm is proposed in this paper. It is the digital fk=f(kT/N). counterpart of thewell-knowncontinuous-timefilterand squaretimingrecovery [SI, [6],hut it extractsthetiming The sequcncc informationfromthesquaredsignal in a new way.The analysis of the timing jitter presented in this paper leads to results that are similar to the continuous-time case, although the method of analysis is different. Another main contribution of the paper is a new method of hangup-freefiltering of thetimingsignal. With allother known timing recovery methods a major problem is that the synchronization loop can get stuck at an unstable equilibrium point. In this paper, we show how this can be avoided through represents the samples of the filtered and squared input signal planar filtering of two-dimensional timing estimates. a spectralcomponentat 1 / 7 . Thisspectral andcontains The final section of the paper presents a digital realization of component, which in a conventional synchronizer is extracted the timing detector whichis suitable for VLSI integration also by a PLL or a narrow-band filter is here determined for every at high data rates. section o f length L 7 (i.e., from L N samples) by computing the complex Fourier coefficient at the symbol ratc Paper approved by the Editor lor Signal Design, Modulation, and Detection of the IEEE Comrnunmtions Society. Manuscript rcceived April 18. 1987; revised September 23. 1987 This paper was presented at GLOBECOM '87, Tokyo, Japan. November 1987. The authors are with Lehrstuhl fiir Elektrische Regelungstechnlk, Aachen University of Technology, West Germany IEEE Log Number 8820399.

k-mLN

As is shown in the next section, the normalized phase tm = - 1 /27r arg (X,,,) of this coefficient is an unbiased estimate for E.

0090-6778/88/0500-0605$01 .OO @ 1988 IEEE

606

IEEk 'TRANSACTIONS COMMUNICATIONS, ON

NO.

VOL. 36,

5 . MAY 19

Using the identity (A6) from Appendix A , we obtain for tl expectation of X LN

Fig. 1 .

~

E [ X ]=

Dlacrete-time filter and square estimator.

I

E [ x ~ ] ~ - J ~ " ~ ' ~ (1 k=O

The sampling rate must be such that the spectral component at 1ITcan still be represented. i.e., N / T > 2 / T . We use N = 4 forpracticalreasons. In thecase of bandwidthefficient modulation with a single-sided bandwidthof less than 1/T, the receiving filter gR(f) also has a single-sided bandwidth of less has a single-sided than 1/T and thusthesquaredsignal bandwidth of lessthan 2 / T . Therefore, with N = 4, the sequence xk completely describes the underlying continuoustime signal.

111. STATISTICS OF THE ESTIMATE In this scction, we compute the statistics of the estimate Fm asa function of thepulseform g ( t ) and the noisepower density NOof the addjtivc noise. We assume m = 0 and omit the index m for the sake of simpler notation. A . Meun The mean of the estimate is

For small variance of the estimates we can linearize the arg- and thus operation.

-1

E [ ; ] = - arg ( E [ X I ) 2a -1

arg

=-

28

(

LN- I

=e--

k=O

1

2*

arg P o ( l / T ) .

(1

Therefore, under the assumption The linearization is valid, of course, only for larg(X)I 4 a. subsequent arg operations, filter the are which toHowever, due discussed in Section rv, this is the only case of interest. We have Io compute expectation Of Ihe signal

E[x,]=E

!ni[

u,,g(kT/N-nT-cT)+A(kT/N) m

I'I

.

[1 nzm

a,,g(kT/N- n T - e T )

(7)

I

2 2 m

n=-a

g R ( t )= gT( - t + a T) (generalized matched filter).

(1

*

(1

P(t)=C a,go(t-nT-ET-arT)+nl(t).

(.

g(t)=g,(t-aT)

with g o ( t ) = g T ( t ) g f ( - t )

Since g o ( t )is symmetrical, the optimal sampling instant is go(? = 0 ) , i.e., for the symbol a,, at

(1

t,p,,n=nT+tT+aT, i.e., the required sampling offset is then yields

+E[Inl(kT/N)l21 =

(1

We then have

Theexpectationmustbetakenwithrespecttothejoint distribution of thc symbols a, and the noise n ( f ) .Noise and symbols are independent of each other. Therefore, and with E[n^(t)]= 0, the cross term of the binomial in (7) vanishes. The remaining terms are

E [ x k ]= E

Po(l/T)=O,

g is an estimate of the clock phase E . But even if is not valid, the meanoft^exactly equals the required sampli offset as we show below. We assume

(E

+ a ) T . Evaluating (

m

m=

E[a,a~]g(kT/N-nT-eT) e4

. g*(kT/N-mT-ET)+E[Jnl(kT/N)12]. (8) With noise power u2 and independently distributed symbolsof mean power 1, we have m

Ig(kT/N-nT- €T)12+u2.

E[Xk]= "=-a

(9)

Since go(t) is symmetrical. thc Fouricr transform rcal . Therefore, we have

in (21)

and thus

E[i]=~+a which is exactly what is required for symbol detection.

c

MEYR: OERDER AND

TIME DIGITAL SQUARE FII.TER AND

607

RECOVERY

B. Variance Here we determine the variance of the random variable i , i.e., the mean square error of the estimation. We assume

The three terms [30a)-c)]representthepartsofthetiming jitter that are generated by (signal X signal), (signal x noise), and (noise x noise) interaction.

C. Conditions for Asymptotically Jitter-Free Timing Recovery In this section, we study the conditions to be fulfilled by the arg P o ( l / T ) = O (24) transmit and receive filters necessaryfor the timing estimateto to simplify thenotation. (It can easily be shown that the results have zero variance in the noiseless case (No = 0). i.e., the are valid for arbitrary E and arg Po.) We then have conditions for the s x s-portion of the variance E=O

var

121 = E [ ? ’ ] 1 l (27r)’ L

u:xs = --

E [(Irn X ) 2 ] (2.)’ ( E [Re X ] ) ’ ’

z -

1

(Im P m ( 1 / m 2 m

(Po(l/T))’

(34)

to be zero. We have Pm(t)=g(t)g*(t-mT)

The latter approximationis valid since the imaginary partof X has zero mean and the variances of both imaginary and real part are small compared to the squared real mean. From (14) and (24) follows

LN E [Re XI =EIX1=--Po(l/T). T

P,(f) = G ( f ) =

lm

(35)

* (G*(- f ) e - j z T f m r )

G ( f - v)G*( - v)e-’2nvmT dv

(36)

-m

(26) and use of the abbreviation

The variancc of the imaginary part is

H(f)=G(l/T-f)G*(-f)

(37)

yields

E [Urn m 2 1

jm H(v)e-jZnVm7dv.

P,(l/T)=

(38)

-m

For real valued g ( t ) , i.e., symmetrical joint transfer function of the transmit and the receive filter, we have

Lh’-1 L N - I

E [ x k x k r ]sin (27rk/N) sin (27rk’/N)

= X=O

(27)

k=O

with

xk=

1

l2

a,g(kT/N- nT)+A(kT/N) n= - m

.

G* ( - f )

=G ( f )

(39)

H(f)=G(l/T-f)G(f).

(40)

That means that H(f) is symmetrical around 1/2T and

(28)

By using some approximations which are valid for large L , the expectation can be computed (Appendix B). If the results are used in ( 2 5 ) , we obtain

Irn P,,,(l/T)= Jm Irn H ( v ) cos (27rurnT) du.

(41)

-m

Therefore, a sufficient condition for zero jitter is

Im H ( f )= 0

(42)

which can be obtained, for example. with with

g ( t ) = g(

-I)

(symmetrical pulse shape)

(43)

and also of course with all linear-phase pulses

g(7+t)=g(7-t)

rn

(44)

as they act like the corresponding symmetrical pulse g(t - T) withanadditionaltimingdelay to = 7. The conditions (43) and ( M ) ,however, show that the optimal receive filter in the bynchronization path is a matched filter

gR(t)=gT(-t).

(45)

These results are valid, of course, only with the approxima111-B, in particular,onlyforlarge tionsmadeinSection estimationintervals L T . Inthecasc of shortintervals,the estimate exhibits jitter, but the spectrum of the jitter hasa zero at the origin and can thus be suppressedby low-pass filtering. True absenceof jitter can be obtaincd in general only by using nonoverlapping pulses. This is in contrast to the conventional continuous-time filter and square timing recovery. In the continuous-time case, the timing is determined by detecting thezeros of the timing wave.

608

IEEE TIRANSACTIONS ONCOMMUNICATIOUS.

-

VOI. 36. NO i M A Y 198

Therefore, true jitter-free timing [recovery is possible if the timing wave exhibits only amplitude jitter,but no phasc jittcr. The latter can be achicved. for example, with locally symmetric pulscs [ 7 ] . In o u r case, however. the estimation is based on samples which have an arbitrary offset fromthe zeros and thus exhibitrandomamplitudefluctuations.Therefore, only asymptotically jitter-free recovery can be obtained.

I

D . Sirnululion Results Fig. 2 shows the variance of the estimatcs i (29) for scvcral estimation intervals L where both transmit and receive filter are fourth-order Butterworth filterswith corner frequency 0.7/ Ixn 5 x” ? 1E-BbT and themodulationformat is 8PSK (solid lines).The markers show the results of Monte Carlo simulations (5000 lE-87 I I I I I I I estimates for cach point). In addition, for L = 64, the three E 5 18 15 28 25 30 35 48 parts of (29) are shown by the dotted lines. The simulations are E/No very close to the theoretical results. Only for L i4 are there errors due to the approximations in thc computation of the Fig. 2. Computed and simulated variance of the timing estimate for sever8 estimation intervals L. Transmit and recelve filter are fourth-orde variance which are valid for large L only. For L = 1 and E / Butterworth, .f, = 0.71T. Dotted lines: s X s, s X n, and n x n-part fo NO= 0 dB the simulation result is smaller than the theoretical L = 64. result. This is due to the finite range of E . The variancc tcnds to 1/12 when iis uniformly distributcd in the estimation range. linear phase Fig. 3 showsthe correspondingcurvesfor filters *ith a transfer function amplitude similar to thc above 1E-82 Butterworth filter. For L z 16 and with moderate E/,%‘”the E simulation results match the theory very well. In particular, the predicted missing of an s X s-portion of the variance can 1E-83 m bc seen. Because the absolute variances are much smallerthan t in Fig. 2 , the effects of the finite observation intcrvalsL T are : li-84 much more visiblc here, especially Tor large E / h i , . c E. Frequency Offset Since we usc a free running sampling oscillator, a frequency offset between transmit and receivetiming may be present resulting in a continuously rising or falling t. In contrast to carrierrecovery,however, the frequencyoffset in timing recovery is very small . . . lo-’ of thesymbolrate). We can thcrefore always find an observation length L = L’ for which we can consider E to be approximately constant. Thcn all considerations of the previous sections apply. For L > L ‘ inspection of the estimation algorithm revealsthat the estimate x is nothing but the avcrageoverestlmatesfromshorter intervals, Therefore, also the m e a n of the estimatc i is just the average of the timing delay t over the observation interval, as long as the variationof e is smallerthan 7 1 2 . Thelatter condition, of course, limits the possible observation length L . Similarly, for small frequency offsets, thc variance or t h e estimates can be expectedtobe nearly independentofthe frequency offset. For larger frequcncy offsets, one would have to examine whether the algorithm behaves like its continuoustime counterpart that exhibits a significant incrcase in timing jittcr in the presence of frequency offset.

Iv. PLANAR FILTERING OF THE ESTrhlATES Due to frequency offset and random variations of the delay E , the observation length L is limited. The estimation, however, can be significantly improved if the knowledge of the statistical properties of E is used to postfilter the estimates. For example, a simple “random walk” model forthetime delay E leads to a first-order Kalman filter. The varianceof the filtered estimates can be computed from the variance of the [X]. unfilteredestimatesandtherandomwalkparamctcrs Since the rangeoftheestimates tm isfinite,thefilter innovation must be reduced to the range [ - 0.5, 0.51 as shown in Fig. 4. With this kind of filtcring. however, the following situation can occur. If the true value e is at a value 0.5 distant from the liltered value t,the estimates i also vary aroundthis value and thus the innovation is at &0.5 and vanishes in thc mean. Then

X

1E-85

?

1E-Eb

\

E

I8

5

15

28

25

38

35

40

€/NO

Fig 3

As Fig. 2 . but Imear-pha.;e filters

-

I

I

I

I

I

i

Fig. 4 . F~ltel-ingor the estimates. the filter is in an unstable equilihrium and it can remain ther some time in spite of the large error. This is a problem which arises also for other models an’ filters and with almost all knownsynchronization mcthods even if they-like the one presented above-could be c a l k “open loop” at first sight.Because of the periodic behaviorc the filrcred delay “hangups” can occur in the filter loop. Relow we present a realization which avoidsthese prok lems. The central idea is to filter a complex phasor instead c the corresponding (periodic) angle. The term X,,, from Fig. is such a phasor. Instead of first determining thc angle of thi phasor and filtering the angle, we can apply a Kalman filter t the phasor itself and use the angle of the filtered phasor t control thc sampling (Fig. 5 ) . In Figs. 6 and 7, a situation for filteringo f e and filtering c X is shown. In the first case a hangupcan occur if the error i Cis approximately 0.5 (corresponding to an angle of T ) . Wit planar filtering, however. the filtered value movcs correct1 (with astep width which depend5 o n t h e filter coefficient present)towardsitsplace.Thus,hangupproblems cannc occur any morc.

x

609

OFRDER ANT) M E Y R DIGITAL FILTER AND SQUARE TIME RECOVER1

I

Fig 6

.

Euarnple for the tra.jectory of estimate and filtered value wlth filtering of the delay values

Fig. 8. Fast digltdl realization of the detector

1 '

0

0x2

0

x3

VI. CONCLUSIONS The proposed timing recovery enables aVLSI realization of digital rcccivcrs which can operate on a sampled input signal without any feedback to the sampling dcvice. The latter can operateata fixed rate with a freerunningoscillator.The planar filtering algorithm results in very fast and hangup-free timing recovery.

x,

Fig. 7 . TIP. 6, but planar filtcring

Thc filtered estimate C,, of Fig. 5 has of course the finite range L - 0.5, 0.51 again and with only small variation ip jumps can occur bctwccn i 0 . 5 In i,,, . The interpolation unit. however, which is controlled by Cm can easily discriminate these "wrap-around" jumps from true variations of the time delay and thereforedeterminetheunderlyinginfiniterange estimate and correctly compute the decision metric [4]. As a final remark. let us note that with the digital filter and square timing estimation, the planar filtering is nothing but a (wcightcd)summationofsucccssivc valucs X,,z:and this is merely anextension of theestimationinterval L T inthe X , (Fig.1) with an additional algorithmforcomputing wcighting of thc tcrms.

z,,!

v. REALIZATION OF THE DETECTOR Fig. 8 shows a possible realizationof the computationofX,,, which allows high dataratesthrough the use of parallel processing and pipelining. the quadruplcs of samples With a doublesetoflatches. belongingtoanestimation interval of length L = 1 are on only values 0 collected. Since the sin and cos functions take and f 1, no multiplications arc nccessary. The samples can then be processed at the symbol rate 1/T rather than at the4 / T sampling rate. Squaring and additioncan be divided by latches into furthcr pipeline stages. Thus, the fact that the estimation algorithm needs 4 samples per symbol (instead of one or two as other algorithms do) is relevant virtually only to the A / D converter and therefore the estimator can be used even at high data rates. We are currently incorporating the detector into a CMOS standard cell chip which will run at about 10 Mbits/s. In caseswhere a lownumber of samplesper symbol is important (e.g., when adaptivc ccho cancellers are used), the actual sampling rate can be reduced to 2 samples per symbol to generate the missing by using asimpleall-passfilter samples [2].

APPENDIX A Equivalcncc of discrete-time and continuous-time computation of the Fouriercoefficients of periodic band-limited functions. Assuming x ( t ) to be a Tperiodic and N/2T-bandlimited signal, we show that

To do so, we start with the integral form. Due to the low-pass limitation, we can rewrite the signal by using sinc-interpolation (sinc x = sin x/x)

610

IEEE TRANSACTIONS ON COMMUNICATIONS, VOL. 36, NO 5 , MAY S I

Inparticular,with ko+L,V-l

x(t) =

_m

y ( t - n T ) , we havc

LN x ( k T / N ) e - J 2 r k " NY =(-l / T ) .

(A6)

r

k-kg

E[anIa;Jg(kT/N-nl T) "2

"I

. g*(kT/N-n2T)E[n(k'T/N)n*(k'T/N)] =No

We first compute the expectation

"I

[(

x

(Bll

E 3E=[ aC, , ~ ~ , l g ( k T / N - n l T )

E[XkXk,1 =E

2 p o (ikTT)/pN(-O )

,=- m

APPENDIX B

"2

. g*(k' T/N-n2 T ) E [ n ( k ' T / N ) n * ( k T / N ) ] a,,g(kT/N-nT)+A(kT/N)

m

g ( k=TN / No- i T ) g * ( k ' T / N - i T )

'1 I

i=

a,,g(kT/N- nT)+ A ( k T / N )

-m

. p((k-k')T/N)

(B1

E,=E[n(kT/N)n*(kT/N)n(k'T/N)n*(k'T/N)]

a,,g(k'T/N-nT)+A(k'T/N)

=Wop2(0)+%p2((k-k')T/N). a,,g(k'T/N-nT)+ri(k'T/N)

1*I

.

xx

LN- I LN- I

S I I=

k=O

With

k'=O

5

m

Era;]= O

=

E[uja;] = 6,,

[

LN-I k=O

y 0

for i = j # k = l

for i = t # j = k for i = j = k = l

sin (27rk/N) sin ( 2 a k ' / N )

E I r n PO(l/T) T LN-I

5

m

SI*=

E[a,a,]= 0

po(kT/N-iT)

/=-a,

I _ - _

. po(k'T/N-jT)

E[a,A(t)]= 0

1 1 E[U,~;Q~U =, ? ]

(B1

Thecorrespondingterms in (27) are now computed. TI approximations are valid for large estimation intervals LT.

k'=O

I=--

(B1: p,(kT/N-iT)

/=-m

. p , ( k ' T / N - i T ) sin (27rk/N) sin (27rk'/N) 036)

2 2 m

ZL

otherwise

m

m

pj(kT/N)pj(k'T/N)

jm

~ l n ( t ) n * ( t + 7 ) 1 = ~ O ( D ( r ) = ~ og R ( t ) g ; ( t + T ) dt -m

(Bl,

037) 10 of thc 16 terms which result from (BI) vanish.The following terms remain:

E[X,X~,]=EI+~E,+~E,+E,

k=O

(B8)

k'=O

. sin ( 2 a k / N ) sin ( 2 a k ' / N )

with

m

E I = ~ 711

"2

E[anla:2an3ai41

. g ( k T / N - n, T)g*(kT/N-n2 T ) . g(k'T/N-n3T)g*(k'T/N-n4T)

m

po(kT/N-iT) k=-m

"4

"3

m

z(y-2)L k'=-m

I=--

. p o ( k ' T / N - i T ) sin ( 2 a k / N ) sin (27rk'IN) = L ( N / T ) 2 ( y - 2 )(Im P o ( l / T ) ) ' = 0 .

039)

Therefore, =

2

po(kT/N-iT)po(k'T/N-jT)

El 1

,=-m,=-m

+

m

m

p , ( k T / N - i T ) p ,E( kI 2 'T/N-iT) i=--,=

x m

+(y-2)

and

m

,= - m

p o ( k T / N - i T ) p o ( k ' T / N - i TE )l

k=O

k'=O

I - - _

. sin (27rk/N) sin ( 2 a k ' / N ) = 0

(B1

61 1

OERDER AND MEYR. DIGITAL FILTER AND SQUARE TIME RECOVERY

with k ’ - 0 z--w

k=O

. rp((k- k ’ )T / N ) sin (27rk/N) sin (27rk’/N) m

m

=LNo

g(kT/N)g*(k‘T/N) k = - m k =-m

. p((k - k ’ ) -L(N/T)’N,

T / Nsin ) ( 2 a k / N )sin (27rk’/N)

[Ie

g(t)g*(f’)p(t-f’) -m

. sin (27rUT) sin (27rt’/T) dt dt’ : = L(N/T)’NOI3

(B23) k=O

h’=O

. sin ( 2 a k / N )sin ( 2 r k ’ / N ) L,V- I

L’v- I

p’((k-k‘)T= / NN) ’ , k=O

k’=O

. sin ( 2 ~ k / Nsin) (21rk’/N) ZLh-2

zNi m=o