Open loop cyclostationarity based timing recovery for accelerated timing acquisition in frequency selective channels

[0001] The present application claims the benefit of Provisional Application Serial No. 60/336,071 filed on Oct. 25, 2001.

[0002] The present invention relates to timing recovery in digital receivers.

[0003] Digital receivers typically sample the received signal, and process these samples in order to perform a variety of functions. It is important in such receivers to sample the received signal at the instants when a data element such as a symbol is present. Accordingly, receivers recover the timing at which data elements are received and use this timing to control sampling of the received signal.

[0004] Timing recovery deals with finding the ideal sample times to recover the received symbols from a received signal. Mismatch between the transmitter and receiver clocks complicates the recovery process because slight differences between these clocks causes a signal to be sampled at the wrong times, contributing to intersymbol interference.

[0005] One popular technique for recovering timing in a receiver is known as the Gardner technique. FIG. 1 illustrates a waveform representing data that might be received at a receiver. In the limited case of FIG. 1, the data is 1 1 1 1. This waveform assumes that the transmitted data is confined to a binary constellation containing only values of 1 and 1, although it can be shown that the Gardner technique works as well for larger constellations. As shown by FIG. 1, the received signal is sampled at two times the data rate.

[0006] The Gardner technique essentially forms a timing error according to b(ac). As long as the sampling is exactly synchronized with respect to the received data, the quantity b(ac) is zero. For example, in the case of a transition from 1 to 1 as shown in FIG. 1, b is zero so that the quantity b(ac) is zero. In the case of no transition so that two adjacent bits have the same value such as 1 and 1 as shown in FIG. 2, b is no longer zero. However, a and c have the same value so that the quantity b(ac) is again zero.

[0007] If the sampling timing is not exactly synchronized to the received data, the quantity b(ac) is not zero but, instead, takes on a value that is a measure of how far off the sampling is from the optimal sampling instants. This value is the timing error and may be used to re-synchronize the sampler to the received data.

[0008] Although the Gardner technique works reasonably well in the case of white Gaussian noise such as electronic thermal noise or atmospheric noise, it does not work as well in the case of multi-path reception of the transmitted signal.

[0009] The present invention relates to the accurate recovery of timing even in the presence of multi-path reception of a transmitted signal.

[0010] In accordance with one aspect of the present invention, a method is provided to determine an error between the rate f_sym at which data are received and the rate f_s at which the data are sampled in a receiver. The method comprises the following: processing a received signal by a nonlinear operator; performing a DFT on the processed signal to produce a plurality of DFT bins each characterized by a respective frequency; determining a dominant spectral component k₀ from at least two of the DFT bins whose frequencies are substantially close to the frequency of the dominant spectral component k₀; and, determining the data rate f_sym from the dominant spectral component k₀.

[0011] In accordance with another aspect of the present invention, a method is provided to correct the sampling rate f_s of a receiver so that the sampling rate f_s matches the rate f_sym at which data are received by the receiver. The method comprises the following: processing a received signal by a nonlinear operator; performing a DFT on the processed signal to produce a plurality of DFT bins each characterized by a respective frequency; determining a dominant spectral component k₀ from two of the DFT bins whose frequencies are substantially close to the frequency of the dominant spectral component k₀; determining an error based on the dominant spectral component k₀, wherein the error indicates an offset between f_s and f_sym; and, adjusting the sampling rate f_s to reduce the error.

[0012] In accordance with still another aspect of the present invention, a method is provided to determine a relationship between a sampling rate f_s at which a receiver samples received data and a data receive rate f_sym at which the data are received by the receiver. The method comprises the following: processing a received signal by a nonlinear operator; sampling the processed received signal to produce signal samples; performing a DFT on the signal samples to produce a plurality of DFT bins each characterized by a corresponding frequency; determining a dominant spectral component k₀ from two of the DFT bins whose corresponding frequencies are equal or closest to f_sym; and, determining the relationship between the sampling rate f_s and the data receive rate f_sym based on the dominant spectral component k₀.

[0013] These and other features and advantages will become more apparent from a detailed consideration of the invention when taken in conjunction with the drawings in which:

[0014] FIGS. 1 and 2 illustrate waveforms useful in explaining a prior art approach to the correction of the sampling rate in a receiver;

[0015] FIG. 3 illustrates an 8 VSB spectrum at baseband;

[0016] FIG. 4 illustrates a timing error detector according to an embodiment of the present invention;

[0017] FIG. 5 illustrates a DFT (Discrete Fourier Transform) of a received signal;

[0018] FIG. 6 illustrates a DFT of a squared received signal;

[0019] FIG. 7 illustrates a DFT of a received signal raised to the fourth power;

[0020] FIG. 8 illustrates an arrangement for correcting the sampling frequency of a receiver based on the timing error detected by the timing error detector shown in FIG. 4; and,

[0021] FIG. 9 illustrates a loop filter that can be used on the arrangement of FIG. 8.

[0022] While the timing estimator disclosed herein is equally applicable to other linear modulation techniques, the present invention is disclosed herein with particular reference to 8 VSB modulation, such as that currently employed in digital television. The 8 VSB signal is a linearly modulated 8-ary PAM signal with real-valued symbols bk and a complex pulse shape p_v(t). The transmitted signal may be expressed by the following equation: $\begin{array}{cc}{x}_v\left(t\right)=\underset{k}{}b\left[k\right]p_v\left(t-{\mathrm{kT}}_{\mathrm{sym}}\right).& \left(1\right)\end{array}$

[0023] The symbols sent every T_sym seconds are assumed to be independent and identically distributed. These symbols are taken from the alphabet having the following amplitudes:

[0024] where 2d is the distance between adjacent symbols. The 8 VSB pulse shape p_v(t) is the complex root raised cosine pulse shape whose spectrum is shown in FIG. 3.

[0025] The 8 VSB signal x_v(t) is passed through a channel denoted c(t) and is received at the receiver to yield a noisy signal according to the following equation:

[0026] where * denotes convolution and (t) is the zero-mean complex additive white Gaussian noise with independent and identically distributed real and imaginary components, each with a variance _n². In all-digital receivers, the received signal is, prior to matched filtering, oversampled by some multiple (such as two) of the nominal symbol rate f_sym(1/T_sym) to give the digital signal r_vn.

[0027] The cyclostationarity technique disclosed herein recovers information about the symbol rate by first passing the received signal through a non-linearity to detect a discrete spectral component corresponding to the symbol rate and then filtering or using some other technique to determine the frequency of the spectral component.

[0028] In one embodiment of the present invention, this non-linearity is squaring. Accordingly, the received signal represented by equation (2) is squared to produce the following equation: $\begin{array}{cc}\begin{array}{c}{r}_v^2\left(t\right)={\left(x_v\left(t\right)*c\left(t\right)+\left(t\right)\right)}^2\\ ={\left(x_v\left(t\right)*c\left(t\right)\right)}^2+2\left(t\right)\left(x_v\left(t\right)*c\left(t\right)\right)+v^2\left(t\right)\end{array}& \left(3\right)\end{array}$

[0029] An interesting result of simply squaring the signal rather than squaring the magnitude of the signal is that the complex-valued noise averages to zero. This result means that the expected value of the noise, v²(t), is zero. Moreover, because the noise is assumed zero-mean and independent of the data, the expected value of the cross terms in equation (3) is also zero. Therefore, taking the expected value removes all terms involving the noise v(t) according to the following equation:

[0030] Combining equations (1) and (4) produces the following equation: $\begin{array}{cc}\left\{r_v^2\left(t\right)\right\}={}_b^2\underset{k}{}{p}_c^2\left(t-{\mathrm{kT}}_{\mathrm{sym}}\right)& \left(5\right)\end{array}$

[0031] where p_c(t)p_v(t)*c(t)

[0032] Equation (5) describes a signal that is periodic with a period T_sym. Therefore, equation (5) can be expressed as a Fourier series according to the following equation: $\begin{array}{cc}\left\{r_v^2\left(t\right)\right\}=\underset{k}{}{c}_k{}^{\mathrm{j2}{\mathrm{kf}}_{\mathrm{sym}}}& \left(6\right)\end{array}$

[0033] where c_k are the Fourier coefficients.

[0034] For a band limited signal, squaring the signal doubles the bandwidth of the signal. For the baseband 8 VSB signal, the baseband spectrum extends from 0.31 MHz to 5.69 MHz. Thus, the squared signal has energy in the frequency interval from 0.62 MHz to 11.38 MHz. Because f_sym10.72 MHz for the 8 VSB signal, all Fourier series coefficients are zero except c₀ and c₁.

[0035] While the relationship between the receiver's sampling rate f_s and the true symbol rate f_sym is not known precisely at the receiver, it is reasonable to assume that the receiver has a good initial estimate of the true symbol rate. Use of a crystal oscillator having an accuracy of 100 parts per million (ppm) means that it can reasonably be assumed that the receiver can determine the symbol rate f_sym of the received signal relative to its own clock within 100 ppm. Based on these assumptions, if a data signal has a symbol rate of 10 MHZ, the symbol rate is known within 1000 Hz at the receiver. However, even this small offset can cause severe error performance over time unless this error is corrected.

[0036] In one embodiment of the present invention, the received signal is oversampled to produce samples r_vn. As shown in FIG. 4, the oversampled signal r_vn is supplied to a timing error detector 10 that includes a bandpass filter 12 operating at half the nominal symbol rate (f_sym/2). The bandpass filter 12 bandpass filters the oversampled signal r_vn. The impulse response of the bandpass filter 12, for example, may be given by the following equation: $\begin{array}{cc}{h}_{\mathrm{bpf}}\left[n\right]=\cos \left(\frac{n}{2K}\right)\cos \left(\frac{n}{U_{\mathrm{os}}}\right)& \left(7\right)\end{array}$

[0037] where K is the filter length and U_OS is the oversampling factor.

[0038] The resulting signal r_bpfn is passed through a non-linear operator 14 that imposes a non-linear operation, such as squaring, on the signal r_bpfn in order to boost the spectral line at the digital frequency corresponding to the symbol rate f_sym relative to the other spectral lines. For example, FIG. 5 shown a DFT of a received signal without the non-linear operation. The symbol rate f_sym is buried among the other frequency components of the received signal and cannot be readily distinguished. FIG. 6 shows the effects on the DFT of squaring the received signal samples. The symbol rate f_sym rises above the other frequency components and, thereby, is more readily distinguished. As also can be seen by comparing FIGS. 5 and 6, the bandwidth has doubled when the samples of the received signal are squared.

[0039] The non-linearity is preferably chosen with the nominal sampling rate in mind so that the signal's spectrum is not aliased near the symbol rate f_sym. For an 8 VSB signal that is squared, sampling at twice the nominal symbol rate is enough to prevent aliasing of the squared signal.

[0040] The signal xn resulting from the non-linear operator 14 is modeled over a small frequency interval around the true symbol rate as a single complex exponential plus white noise according to the following equation:

[0041] where $\begin{array}{cc}z\left[n\right]={}^{j\frac{2k_{0}n}{N}}& \left(8\right)\end{array}$

[0042] where is a complex scalar, and where vn is the noise component. It is then possible to compute two Discrete Fourier Transform (DFT) values {tilde over (Z)}(k) and {tilde over (Z)}(k1) near the nominal value of k₀ that corresponds to the rate f_sym at which data is received, to determine k₀ from these values, and to determine the relationship between f_sym and f_s from k₀. Subspace based averaging may be used to determine these values such that the effect of noise is minimized.

[0043] Accordingly, a timing error estimator 16 performs a DFT on xn. A conjugate centrosymmetrized version of the kth bin of the resulting N-point DFT is given by the following equation: $\begin{array}{cc}\tilde{X}\left(k\right)=\underset{n=0}{\overset{N-1}{}}x\left[n\right]{}^{-j\frac{2k}{N}\left(n-\frac{N-1}{2}\right)}& \left(9\right)\end{array}$

[0044] where xn are the received signal samples. Equation (9) can be re-written as a product of vectors according to the following equation:

[0045] where $\begin{array}{cc}{\tilde{w}}_k^H=\left[{}^{j\frac{2k}{N}\left(\frac{N-1}{2}\right)}{}^{j\frac{2k}{N}\left(\frac{N-3}{2}\right)}{}^{-j\frac{2k}{N}\left(\frac{N-1}{2}\right)}\right]C^N& \left(11\right)\end{array}$

[0046] and where

[0047] C^N represents the N dimensional complex vector space comprising complex vectors each of length N. The rotated DFT vector {tilde over (w)}_k is termed conjugate centrosymmetric because its mth element is equal to the conjugate of its (N-m)th element for m0, . . . , N1. {tilde over (X)}(k) may be found by first finding the DFT value X(k) and by then post multiplying this value by ${}^{j\frac{2k}{N}\left(\frac{N-1}{2}\right)}.$

[0048] If zn is the single complex exponential in equation (8), then the kth bin of the conjugate centrosymmetrized DFT of zn evaluated at DFT bin k is given by the following equation: $\begin{array}{cc}\tilde{Z}\left(k\right)=\frac{\sin \left(\left(k-k_0\right)\right)}{\sin \left(\frac{}{N}\left(k-k_0\right)\right)}{}^{j\frac{2k_0}{N}\left(\frac{N-1}{2}\right)}& \left(13\right)\end{array}$

[0049] and the conjugate centrosymmetrized DFT of zn evaluated at DFT bin k1 is given by the following equation: $\begin{array}{cc}\tilde{Z}\left(k+1\right)=\frac{\sin \left(\left(k+1-k_0\right)\right)}{\sin \left(\frac{}{N}\left(k+1-k_0\right)\right)}{}^{j\frac{2k_0}{N}\left(\frac{N-1}{2}\right)}& \left(14\right)\end{array}$

[0050] Using the trigonometric identity sin(AB)sin(A)cos(B)cos(A)sin(B), the numerator of equation (14) becomes sin((kk₀)). Therefore, dividing equation (13) by equation (14) yields the following equation: $\begin{array}{cc}\frac{\tilde{Z}\left(k\right)}{\tilde{Z}\left(k+1\right)}=-\frac{\sin \left(\frac{}{N}\left(k+1-k_0\right)\right)}{\sin \left(\frac{}{N}\left(k-k_0\right)\right)}& \left(15\right)\end{array}$

[0051] Applying the trigonometric identity given above to equation (15) produces the following equation: $\begin{array}{cc}{k}_0=\frac{N}{}{\tan }^{-1}\left(\frac{\tilde{Z}\left(k\right)\sin \left(\frac{}{N}k\right)+\tilde{Z}\left(k+1\right)\sin \left(\frac{}{N}\left(k+1\right)\right)}{\tilde{Z}\left(k\right)\cos \left(\frac{}{N}k\right)+\tilde{Z}\left(k+1\right)\cos \left(\frac{}{N}\left(k+1\right)\right)}\right)& \left(16\right)\end{array}$

[0052] Based on k₀, the true symbol rate, i.e., the symbol rate of the received data, can be determined from the following equation: $\begin{array}{cc}{f}_{\mathrm{sym}}=\frac{k_0}{N}{f}_s& \left(17\right)\end{array}$

[0053] where N is the number of bins in the DFT and f_s is the sampling frequency of the receiver.

[0054] The relationship between f_sym and f_s should ideally be the design relationship. For example, if the sampling rate is nominally set at twice the symbol rate used by the transmitter, then ideally f_symf_s/2. If the actual symbol rate of the received data, as determined by equations (16) and (17), and f_s do not have this ideal relationship, then the sampling rate of the receiver is adjusted until the ideal relationship is achieved.

[0055] The timing error estimator 16 can also determine the timing error by subtracting k₀ from k. This timing error may be used to appropriately adjust the sampling frequency of the receiver in order to correctly sample the received signal.

[0056] As discussed above, the non-linear operator 14 may be arranged to square the signal from the bandpass filter 12. Instead, other non-linearities can be imposed on the received signal. For example, the received signal may be raised to the fourth power, i.e., (r_vn)⁴. FIG. 7 shows the effects on the DFT of raising the received signal samples to the fourth power. The symbol rate f_sym rises even more above the other frequency components and, thereby, is even more readily distinguished than is the case of simply squaring the samples. As also can be seen by comparing FIGS. 5 and 7, the bandwidth has quadrupled when the samples of the received signal are raised to the fourth power. For an 8 VSB signal that is raised to the fourth power, sampling at four times the nominal symbol rate is enough to prevent aliasing of the signal raised to the fourth power. When the non-linear operator 14 raises the received signal samples to the fourth power, the bandpass filter 12 becomes unnecessary.

[0057] It is preferable for the process described above that values for N and k be chosen so that k<k₀<k1 and so that N<fs/(2f) where f is the tolerance of the a priori estimate. For example, in the case where the samples are raised to the fourth power so that the sampling is at about four times the symbol rate, N can be set at 4096, and k can be set at 1024 representing the nominal symbol rate f_sym. The frequency associated with bin k1 is the next higher frequency. The bin k1 can be used instead of the bin k1.

[0058] FIG. 8 shows an arrangement for correcting the sampling frequency of the receiver. The received signal is re-sampled by a re-sampler 22 that receives the samples from an upstream sampler, that effectively reconstructs an analog signal from the samples, and that re-samples the reconstructed analog signal in response to a sampling frequency supplied by a numerically controlled oscillator 24. The samples at the output of the re-sampler 22 are processed by a timing error estimator 26 to detect the error between f_s and f_sym as described above, and this error is filtered by a loop filter 28. The filtered error from the loop filter 28 then controls the numerically controlled oscillator 24 to produce the correct sampling frequency for the re-sampler 22. A matched filter may be located in the output of the re-sampler 22 upstream of the timing error estimator 26, or the matched filter may instead be located upstream of the re-sampler 22.

[0059] FIG. 9 shows an example of a loop filter 40 that can be used for the loop filter 28 of FIG. 8. The loop filter 40 is a low pass filter that smoothes the error provided by the timing error estimator 26. The smoothed error is denoted in FIG. 9 as ln. The quantities ₁ and ₂ are constants.

[0060] If the noise component vn in xn is removed, {tilde over (X)}(k,n) and {tilde over (X)}(kl,n) can be used in equation (16) in place of {tilde over (Z)}(k) and {tilde over (Z)}(k1) because multiplication by a scalar has no effect on equation (16). One way of removing the noise component in xn is to perform subspace-based, or eigenvalue-based, averaging in the frequency domain.

[0061] Let {tilde over (X)}(k,n) and {tilde over (X)}(k1,n) be conjugate centrosymmetric N point DFT values of xnN1, . . . , xn at frequency bins k and k1, respectively. These DFT values can be grouped as a vector X{tilde over (X)}(k,n){tilde over (X)}(k1,n)^T. Then, the real-valued auto-correlation matrix R_xx of the received signal is simply the following expected value: XnX^Hn. This matrix has the following form: $\mathrm{Re}\left\{{\hat{R}}_{\mathrm{xx}}\left[n\right]\right\}=\left[\begin{array}{cc}a& b\\ c& d\end{array}\right]$

[0062] The eigenvectors associated with this matrix have the following form: $e_i={\left[\begin{array}{cc}\frac{{}_i-d}{c}& 1\end{array}\right]}^T{C}^l$ $e_i=[\begin{array}{cc}\frac{b}{{}_1-a}& {1]}^T{C}^l\end{array}$

[0063] where _l are eigenvalues given by the following equation: ${}_i=\frac{a+d\sqrt{{\left(a-d\right)}^2+4\mathrm{bc}}}{2}$

[0064] The eigenvector associated with the largest eigenvalue of the real part of the sample correlation matrix {circumflex over (R)}_xxn is a multiple of Z{tilde over (Z)}(k) {tilde over (Z)}(k1)^T. Therefore, these values of {tilde over (Z)}k and {tilde over (Z)}k1 can be derived from this eigenvector and used in equation (16) to compute k₀ which can be used in equation (17) to determine the relationship between f_s and f_sym.

[0065] Because R_xx is real valued, only the real part of the vector product XnX^Hn need be included when forming the sample correlation matrix. Because {circumflex over (R)}_xxn is

[0066] a 22 matrix, a simple closed form expression is available for the eigenvectors.

[0067] The first DFT values {tilde over (X)}(k,N1) and {tilde over (X)}(k1,N1) may be computed gradually (without the need to store or compute all of the twiddle factors $w_N^{\mathrm{nk}}={}^{-j\frac{2kn}{N}},$

[0068] n0, . . . , N1) using the well known Goertzel algorithm. Subsequent values of {tilde over (X)}(k,n) and {tilde over (X)}(k1,n) for n>N1 may be found recursively. If X(k,n) is the kth bin of the N point DFT of the current and past N1 values of xn, then X(k,n) is given by the following equation: $\begin{array}{cc}X\left(k,n\right)=\underset{m=0}{\overset{N-1}{}}x\left[m+n-\left(N-1\right)\right]{}^{-\mathrm{j2}k\frac{m}{N}}& \left(18\right)\end{array}$

[0069] and X(k,n) may be expressed in terms of X(k,n1) according to the following equation: $\begin{array}{cc}\begin{array}{c}X\left(k,n\right)=\left(X\left(k,n-1\right)-x\left[n-N\right]\right){}^{j2\frac{k}{N}}+\\ x\left[n\right]{}^{-j2\frac{k\left(N-1\right)}{N}}\end{array}& \left(19\right)\end{array}$

[0070] for n>N1. Use of equation (19) requires only two complex multiplications and two additions, which is a significant computational savings over equation (18).

[0071] Furthermore, when k is an integer, equation (19) simplifies even more and requires only one complex multiplication and two additions according to the following equation: $X\left[k,n\right]=\left(X\left[k,n-1\right]-x\left[n-N\right]+x\left[n\right]\right){}^{j\frac{2k}{N}}$

[0072] Accordingly, when averaging is used as described above, as each new sample is received, the sample is passed through the non-linear operation. Two conjugate centrosymmetric DFT values are found using the current sample and the previous N1 samples. These values are used to update the real part of the sample correlation matrix, and the eigenvector having the largest eigenvalue is extracted to give a scaled estimate of {tilde over (Z)}(k) and {tilde over (Z)}(k1). These values are used as per equations (16) and (17) to determine the symbol frequency of the received data.

[0073] Modifications of the present invention will occur to those practicing in the art of the present invention. Accordingly, the description of the present invention is to be construed as illustrative only and is for the purpose of teaching those skilled in the art the best mode of carrying out the invention. The details may be varied substantially without departing from the spirit of the invention, and the exclusive use of all modifications which are within the scope of the appended claims is reserved.