Block decision feedback equalization method and apparatus

[0002] Cross-reference is made to a related U.S. patent application by Joseph T. Graf and Thomas L. Tapp, filed on even date herewith, and entitled “Block Decision Directed Equalization Method and Apparatus.”

[0003] The present invention generally relates to an improved data communication system and technique, more particularly relates to an improved equalization scheme for a digital communication system, and even more particularly relates to the block decision feedback equalization of signals transmitted via a dynamic dispersive media.

[0004] Equalization is a technique often used in wireless communication applications. For example, when digital data is transmitted over a dispersive communication channel, an equalizer can be used to improve the estimate of the originally transmitted symbols. More particularly, equalizers are becoming increasingly important in high frequency (HF) radio communication applications. Such applications are requiring increasingly higher data rates in a fixed bandwidth; i.e. increased bandwidth efficiency. For example, many current applications require performance sufficient to meet the specifications for a 64-QAM (quadrature-amplitude modulation) signal.

[0005] As bandwidth efficiency is increased, however, information reliability is typically reduced. The reduction in reliability results, at least in part, from various channel effects. Equalization techniques can be used to mitigate many of the various channel effects. Use of an equalization technique can thereby make the communication system more reliable.

[0006] Several different channel effects can cause a reduction in reliability. A description of some of the various channel effects follows. Constrained channel bandwidth, for example, can cause time dispersion and intersymbol interference (ISI). Thermal noise and interference can lower the signal-to-noise ratio and thereby impact performance. Multipath interference can cause fading and time dispersion. Fading can also result from the physical characteristics of the communication channel being used. Delay and Doppler shift, which also result from the physical characteristics of the communication channel being used, can cause, respectively, time and frequency shifts.

[0007] Transmitter and receiver characteristics, which also affect the signal, must be accounted for as well. For example, the transmitter can include transmit filters and automatic level control (ALC). Likewise, the receiver can include receive filters and automatic gain control (AGC). These system influences should be addressed since they are also included in the channel response of any real system.

[0008] The equalizer's role is to compensate for such effects. FIG. 1 depicts a block diagram of a basic communication system. Data from a source 100 is fed to a modulator 102 that in turn passes the modulated data to a wireless communication component 104. The communication signal is then transmitted by a transmitter antenna 106 to a receiver antenna 108. The link between the two antennas 106, 108 is known as the propagation channel 110. The receiver antenna 108 passes the received signal to the receiver's wireless communication component 112 that in turn sends the signal to a receiver-based demodulator 114. The demodulated signal is then sent for further processing, handling or storage as represented by sink 116.

[0009] The channel between the two wireless communication components 104, 112 is known as the radio channel 118. The channel between the modulator 102 and the demodulator 114 is known as the modulation channel 120. The channel between the source 100 and the sink 116 is known as the digital channel 122. The equalizer typically compensates for the influences found across the modulation channel 120. Thus, when used for such a purpose, it can be used to compensate for the channel and system effects described above.

[0010] FIG. 2 depicts a use of an equalizer as a filter compensating for the effects influencing the modulation channel. The effects of the transmitter 200 can be represented as H_Tx(f). The effects of the propagation channel 202 can be represented as C(f). The effects of the receiver 204 can be represented as H_Rx(f). For example, the ideal response of a linear equalizer 206 can be represented as a filter that mitigates the effects of the modulation channel by having a response such as [H_Tx(f)C(f)H_Rx(f)]⁻¹.

[0011] Typical equalizer implementations require knowledge of the characteristics of the communication channel. Several different techniques have been used to estimate or model these characteristics. For example, the initial channel characteristics can be estimated using a training sequence. The estimate can then be updated using a least-mean squares (LMS) algorithm or a recursive-least squares (RLS) algorithm. Another technique employs periodically inserted correlator probes, which are used to guide the channel updates. Further, use of a carefully designed probe can enable an efficient estimation of the channel characteristics. One example of such a probe is provided in the MIL-STD-188-110B Appendix C standard, wherein a simple correlator can be used to obtain the channel estimate.

[0012] An example of a typical signaling sequence is depicted in FIG. 3. An initial preamble 300 of “x” symbols is used for acquisition and synchronization. In this example, the data packets 302, 304 each include thirty-two unknown symbols. The probes 306, 308 each include a sequence of sixteen known symbols. Alternatively, different numbers of symbols can be communicated in the unknown data segments 302, 304 and probes 306, 308. For example, in the more recently developed MIL-STD-188-110B Appendix C standard, each known probe segment includes thirty-one known symbols and each unknown data segment 302, 304 includes 256 unknown symbols. Each “unknown” data segment 302, 304 is thus encapsulated between two “known” probe packets 300, 306, 308. Various modulation techniques, such as phase-shift keying (PSK), quadrature amplitude modulation (QAM) and quadrature phase-shift keying (QPSK), for example, can be used to transmit the various segments of the signaling sequence.

[0013] An example of a communication system using an encapsulated unknown data/known probe approach, such as that depicted in FIG. 3, can be found in U.S. Pat. No. 4,365,338 (Technique For High Rate Digital Transmission Over A Dynamic Dispersive Channel) to McRae et al. The technique disclosed in U.S. Pat. No. 4,365,338, however, uses a normal equations approach to solve the over-determined least-squares problem. New methods and approaches are needed for determining the symbols of the unknown data segment in environments wherein the channel characteristics can vary over the duration of the unknown data segment.

[0014] Consequently, there exists a need for equalization approaches capable of taking changing characteristics into consideration when determining symbols transmitted in an unknown data segment. These needs and others are fulfilled by the invention disclosed in the following detailed description.

[0015] It is an object of the present invention to provide an equalizer for a communication system transmitting data via a dispersive media.

[0016] It is one feature of the present invention to utilize QR factorization to solve the over-determined least-squares problem.

[0017] It is an advantage of the present invention to provide an accurate and reliable determination of the symbols of the unknown data segment even in situations wherein the channel characteristics vary over the duration of the unknown data segment.

[0018] Thus, the present invention involves an improved block equalization apparatus and method that relies on a QR factorization approach to the least-squares detection problem, rather than the normal equations.

[0019] The invention may be more fully understood by reading the following description of the invention, in conjunction with the appended drawings wherein:

[0020] FIG. 1 is a block diagram of the basic structure of a wireless communication system.

[0021] FIG. 2 is a block diagram depicting the some of the basic components of a wireless communication system and the filtering effects associated with those components.

[0022] FIG. 3 depicts a waveform having unknown data segments separated by known data segments.

[0023] FIG. 4 depicts a more detailed view of a decision block structure having unknown data segments separated by known data segments.

[0024] FIG. 5 depicts a convolution matrix related to an embodiment of the present invention.

[0025] FIG. 6 is a flowchart showing steps involved in an embodiment of the present invention.

[0026] FIG. 7 depicts an embodiment using a T/2 spaced fractional sampling scheme.

[0027] FIG. 8 depicts a channel impulse response estimate with a time resolution of T/2.

[0028] FIG. 9 depicts a representation of a matrix for use with a T/2 spaced fractional sampling scheme.

[0029] FIG. 10 depicts another representation of a matrix for use with a T/2 spaced fractional sampling scheme.

[0030] FIG. 4 depicts an embodiment of a decision block structure of a type suitable for use with the present invention. Various known communication components can be used to transmit the data and other information described throughout this specification. Systems such as those outlined in FIGS. 1 and 2, for example, as well as many other known communication systems, can be used in conjunction with the present invention. Such systems typically include a transceiver component and a processing component. Also included is a memory system storing, among other things, received data, data to be transmitted and programs controlling data transmission, data reception and network operations.

[0031] The communicated sequence includes alternating blocks of known and unknown symbols such that each unknown data block or segment is both preceded and followed by a segment of symbols that is known to the receiver prior to, or at the time of, reception. For example, in FIG. 4, each known data segment 400, 402 is labeled “probe” and each unknown data segment 404 is labeled simply as “data.” A known data segment 400, 402 can be a concatenation of a PROBE A 406 portion (from 0 to Y−2) with a PROBE B portion 408 (from Y−1 to N−1). The variable “Y” is related to the length of the channel impulse response such that Y+1 equals the length of the channel response. In addition, the length of the known data segment is represented as “N” and the length of the unknown data segment is “M” in FIG. 4.

[0032] The known and unknown data segments can be further subdivided. Each probe or known data segment can be subdivided into a first portion, “PROBE A” 406, and a second portion, “PROBE B” 408. Further, each PROBE A 406 portion can be subject to ISI 410 from symbols, either preamble symbols or unknown data segment symbols, being transmitted just prior to PROBE A 406. Likewise, each PROBE B 408 portion can be affected by ISI 412 from the PROBE A 406 portion. This effectively isolates the PROBE B portion of the known data segment from the channel effects of the preceding decision block (assuming the channel response is no longer than the length of PROBE A).

[0033] Each unknown data segment 404 can, for discussion purposes, also be subdivided. In FIG. 4, one portion of each unknown data segment is labeled “post-data” 414 and a second portion is labeled “pre-data” 416. Each “post-data” portion can be subject to ISI 418 from PROBE B and from other data symbols of the unknown data segment 404. Each “pre-data” portion can be subject to ISI 420 from “post-data” symbols and from other data symbols within the unknown data segment 404.

[0034] For the equalizer, the goal is to estimate the received segment of unknown data symbols, which is surrounded by known probe sequences. As is well known in the art, several different channel estimation algorithms, such as the least squares, LMS, RLS, Fast RLS, and Lattice algorithms, can be used to generate a channel estimate from a known probe segment. The channel impulse response can be modeled as a transversal filter with Y+1 taps,

[0035] f=[f₀,f₁, . . . ,f_Y]^T,

[0036] where the superscript ‘T’ indicates (real) transposition. Except where noted, the sampling rate is the same as the baud rate.

[0037] Block equalization is based on the least-squares solution to the matrix equation M₁{circumflex over (B)}=C, where matrix M₁is constructed from the estimate of the channel impulse response, {circumflex over (B)} is the unknown data segment to be detected, and the vector C is a function of all of the received samples affected by the unknown data vector B. As is well known in the art, several different channel estimation algorithms, such as the least squares, LMS, RLS, Fast RLS, and Lattice algorithms, can be used to generate f. In U.S. Pat. No. 4,365,338 to McRae et al., the channel estimate is considered to be fixed over the duration of the unknown data segment.

[0038] In general, using the reversed vector,

[0039] ^T=[f_Y,f_Y−1, . . . ,f₀],

[0040] wherein Y is related to the length of the channel response such that Y+1 equals the length of the channel response, the convolution can be represented as r′=MB+n or in matrix form as: $r^{\prime }=\left[\begin{array}{ccccc}{\underset{\_}{f}}^T& & & & 0\\ & {\underset{\_}{f}}^T& & & \\ & & {\underset{\_}{f}}^T& & \\ & & & ⋰& \\ 0& & & & {\underset{\_}{f}}^T\end{array}\right]\left[\begin{array}{c}{A}_3\\ B\\ A_4\end{array}\right]+n$

[0041] where r′ consists of all received symbols affected by the unknown data segment B, and n is a vector of additive white Gaussian noise samples. The last “Y” symbols of the preceding probe sequence form the vector A₃, while the first “Y” symbols of the trailing probe sequence form the vector A₄. These are the probe symbols that cause intersymbol interference (ISI) with symbols in the unknown data segment.

[0042] In addition, the channel estimate can vary over the duration of the unknown data block. One simple approximation to the channel impulse response within the data segment can be obtained by a linear interpolation between the channel estimates obtained by the probes at either end. Other approximations may include quadratic or cubic interpolations using the channel estimates from the probes bordering the adjacent data segments.

[0043] The convolution matrix can be partitioned as shown in the equation r′=M₁B+M₂A₃+M₃A₄+n, where $M_1=\left[\begin{array}{cccc}{f}_0& & & 0\\ f_1& f_0& & \\ ⋮& ⋮& ⋰& \\ f_Y& f_{Y-1}& & f_0\\ & f_Y& & f_1\\ & & ⋰& ⋮\\ 0& & & f_Y\end{array}\right],\left(\left(M+Y\right)\times M\right)$ $M_2=\left[\begin{array}{cccc}{f}_Y& f_{Y-1}& \dots & f_1\\ & f_Y& & f_2\\ & & ⋰& ⋮\\ & & & f_Y\\ & & & \\ 0& & & \end{array}\right],\left(\left(M+Y\right)\times Y\right)$ $M_3=\left[\begin{array}{cccc}& & & 0\\ & & ⋰& \\ f_0& & & \\ f_1& f_0& & \\ ⋮& & ⋰& \\ f_{Y-1}& f_{Y-2}& \dots & f_0\end{array}\right],\left(\left(M+Y\right)\times Y\right).$

[0044] The numbers in parentheses indicate the number of rows and columns of the associated matrix. Note that M₂consists of the first “Y” columns of the original convolution matrix, while M₃consists of the last “Y” columns. The M₁matrix consists of the columns located between those columns used for the M₂and M₃matrices.

[0045] FIG. 5, for example, depicts a representation of an M matrix 500 and the equation MX=r′ wherein Y=1 (and therefore f^T=[f₁, f₀]). The superscripts on the f coefficients indicate the time of the channel estimate. For example, f^0T=[f⁰₁, f⁰₀] is the estimate of the channel for the time r₀is received. In this example, the first column 502 of the M matrix 500 becomes the M₂matrix, the second through fourth columns 504 become the M₁matrix and the fifth column 506 becomes the M₃matrix. The X matrix is composed of the last Y (in this example, 1) symbols of the expected “left-hand” probe, A₃, the unknown data segment B (which in this case includes designations for three unknown symbols b₀, b₁and b₂to be determined) and the first Y (in this example, 1) symbols of the expected “right-hand” probe segment, A₄. The A₃and A₄probe segments are replicas, generated by the receiving apparatus, of portions of the probes expected to be transmitted in the communication sequence. In some applications the preceding and trailing probes will be congruent (each having the same sequence of symbols) while in other applications they will not be congruent (each being composed of a different sequence of symbols). Further, it will be appreciated that “Y” can be any of several different values and that the actual value used will be a function of the length of the channel response.

[0046] The channel estimate is used to construct the M₂and the M₃matrices, which are in turn used to subtract the ISI introduced by the probe symbols. The “precursor” and “postcursor” ISI is subtracted from the receive vector r′ to form

[0047] C=r′−M₂A₃−M₃A₄,  Equation 3.

[0048] Using Equation 3, the least squares solutions for {circumflex over (B)} can be obtained from

[0049] C=M₁{circumflex over (B)},  Equation 4.

[0050] The equations lead one to determine the {circumflex over (B)} that minimizes the mean-squared error (ignoring noise),

[0051] e=(M₁{circumflex over (B)}−C)^*T(M₁{circumflex over (B)}−C).

[0052] Prior approaches, such as U.S. Pat. No. 4,365,338 to McRae et al., have obtained the least squares (LS) solution {circumflex over (B)} by solving the normal equations

[0053] (M₁^*TM₁){circumflex over (B)}=M₁^*TC,  Equation 5,

[0054] where the asterisk denotes the complex conjugate operator and wherein M₁^*TM₁has dimension m×m (wherein m equals the number of transmitted symbols).

[0055] The present invention uses a different method to solve the over-determined least-squares problem. Rather than solving the normal equations, the present invention relies on QR factorization. An overview of QR factorization can be found, for example, in Section 5.3.3 of G. H. Golub and C. F. van Loan, Matrix Computations,3^rd. Ed., Baltimore, Johns Hopkins University Press, 1996. In the present invention, the least squares solution, {circumflex over (B)}, to M₁{circumflex over (B)}=C, is found via the QR factorization of M₁.

[0056] Any permutation of the rows of the M₁matrix with the identical permutation of the rows of the C vector should be recognized as a mere re-ordering of the underlying arithmetic equations, and identical mathematically to the original matrix equation. Similarly, any permutation of the columns of the M₁matrix with the identical permutation of the rows of the {circumflex over (B)} vector should be recognized as a mere re-ordering of the terms in each of the underlying arithmetic equations, and again, identical mathematically to the original matrix equation. There may be implementation advantages, however, to applying such permutations before performing the QR factorization algorithm.

[0057] FIG. 6 outlines the steps performed in the QR factorization method of the present invention. The elements of matrices Q and R are found such that Q is unitary (Q^*T=Q⁻¹), R is upper triangular (or alternatively, as discussed herein, lower triangular), and Q·R=M₁600. Since M₁has dimension (M+Y)×M, it follows that Q has dimension (M+Y)×(M+Y) and R has dimension (M+Y)×M.

[0058] After performing the factorization, R₁is extracted from R 602, where $R=Q^{*T}{M}_1=\left[\begin{array}{c}{R}_1\\ 0\end{array}\right]\begin{array}{c}M\mathrm{rows}\\ Y\mathrm{rows}\end{array}.$

[0059] Next, c is computed 604 using $Q^{*T}C=\left[\begin{array}{c}c\\ d\end{array}\right]\begin{array}{c}M\mathrm{coefficients}\\ Y\mathrm{coefficients}\end{array},$

[0060] wherein Q^*TC is a length (M+Y)vector. The lower case c is used to designate the first M components, which are the components of interest here. R₁is the upper square portion of R. Finally, using back substitution, the upper triangular equation, R₁{circumflex over (B)}=c, is solved for {circumflex over (B)} 606 to yield the least squares solution. Since R₁is upper triangular, the solution is simplified. If desired, when solving the triangular matrix equation R₁{circumflex over (B)}=c 606, instead of feeding back the computed values of {circumflex over (b)}₁without modification, the “decisions” based on the computed values can be fed back.

[0061] As an example, consider the 3×3 upper triangular matrix equation Rb=C, that is, $\left[\begin{array}{ccc}{r}_{1,1}& r_{1,2}& r_{1,3}\\ 0& r_{2,2}& r_{2,3}\\ 0& 0& r_{3,3}\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]=\left[\begin{array}{c}{c}_1\\ c_2\\ c_3\end{array}\right].$

[0062] The bottom row of this matrix equation is the simple equation r_3,3b₃=c₃, with the solution b₃=c₃/r_3,3. The second row gives the equation r_2,2b₂+r_2,3b₃=c₂. Since b₃is now known, we can solve for b₂. That is b₂=(c₂−r_2,3b₃)/r_2,2. Finally, the top row yields r_1,1b₁+r_1,2b₂+r_1,3b₃=c₁. From this equation, we get b₁=(c₁−r_1,2b₂−r_1,3b₃)/c₃, using the previously computed values of b₂and b₃. This procedure can be extended to arbitrary dimension matrices. It can also be extended to lower triangular matrices, simply by starting with the top row and working down, rather than starting at the bottom and working up.

[0063] The vector b represents the transmitted symbols that we are attempting to detect. After computing the algebraic solution for each b_i, a “decision” of the closest symbol in the signal constellation can be made, {circumflex over (b)}_i, and used in the solution for all subsequent b components. In the example above, after first computing b₃, the closest constellation point to this value is determined, which will be denoted {circumflex over (b)}₃. Then, {circumflex over (b)}₃, rather than b₃, is used to solve for b₂. That is, b₂=(c₂−r_2,3{circumflex over (b)}₃)/r_2,2. Similarly, a decision is made on this computed value of b₂, denoting the closest constellation point {circumflex over (b)}₂. The calculation of b₁uses both {circumflex over (b)}₂and {circumflex over (b)}₃. (i.e., b₁=(c₁−r_1,2{circumflex over (b)}₂−r_1,3{circumflex over (b)}₃)/c₃) It will be appreciated that this procedure can be extended to arbitrary dimension matrix equations and to lower triangular matrix equations.

[0064] In a related embodiment, a factorization is used wherein Q is still unitary, but R is lower triangular, such that QR=M₁. The performance of this type of factorization, which can be alternatively referred to as a QL factorization (wherein L designates the use lower triangular matrix), is equivalent to that of the upper triangular factorization embodiment. It will be appreciated that implementation of this embodiment just requires some changes to the indices used in the upper triangular factorization embodiment.

[0065] By way of further example, three methods for computing the QR (or a related) factorization of M₁will be presented herein. The approaches presented below are based on Householder transformations, Givens rotations, and the “fast” Givens method (which avoids using square roots). These are described in detail in Chapter 5 of G. H. Golub and C. F. van Loan, Matrix Computations,3^rd. Ed., Baltimore, Johns Hopkins University Press, 1996, and they can be readily modified to handle complex numbers.

[0066] It will be appreciated, however, that other methods of QR factorization also exist and can be used in the present invention. For example, other approaches to computing the QR factorization include the modified Gram-Schmidt algorithm and singular value decomposition (SVD). It is intended that embodiments using the other, undescribed methods also be considered to be within the scope of the present invention.

[0067] Complexity can be greatly reduced by exploiting the fact that the M₁matrix has a regular banded structure, containing many zeros. In addition, the Q matrix need not be explicitly computed as indicated in the first step 600. Instead, the first 600 and third 604 steps of the algorithm outlined above can be combined by computing the vector Q^*TC and solving for R at the same time.

[0068] First, the QR least squares solution using Householder transformations will be described. A Householder transformation is a square matrix H of the form $H=I-\frac{2}{v^{*T}v}v{v}^{*T}$

[0069] for some vector v, where I is the identity matrix. If v is carefully chosen, it can be arranged that, for some vector x, Hx is zero in all but the first coordinate. The v that achieves this goal has its first coordinate as $v_1=x_1\pm \frac{x_1}{\left|x_1\right|}\sqrt{x^{*T}x},$

[0070] while the rest of its components are given by

[0071] v_i=x_i, i≠1.

[0072] By choosing the x vector appropriately, a QR factorization can be computed as a series of Householder transformations. A MATLABm-file that uses Householder transformations to compute R and Q^*TC is shown below. It does not rely on any black-box MATLABfunctions (although it does use the ‘\’ operator to perform back substitution, since this will be quicker than using ‘for’ loops). The variables A and b correspond to M₁and C, respectively. The structure of the M₁matrix (specifically its known regions of zero components) is exploited to greatly reduce the computational burden. The variable R corresponds to the matrix of the same name, and Q^*TC is accumulated in the variable, b1. The variable taps is the length of the channel estimate, Y+1.

[0073] The following is an example of a MATLABm-file using Householder transformations:

[0074] function x=housels_bandsp2(A,b,taps);

[0075] %househoulder LS Solution to overdetermined Ax=b

[0076] %lower bw=taps−1,upper bw=0

[0077] %#rows=#cols+taps−1

[0078] % A corresponds to M1, b to C, b1 to Q^*TC, R to R, x to bhat=LS soln.

[0079] [m,n]=size(A);

[0080] b1=b;

[0081] for j=1:n

[0082] x=A(j:j+L,j);

[0083] x1u=x(1)/abs(x(1)); %unit vector in direction of x(1)

[0084] v1=x(1)+x1u*sqrt(x′*x);

[0085] v=[1;x(2:L+1)/v1];

[0086] beta=2/(v′*v);

[0087] vvp=v(1:L+1)*v(1:L+1)′;

[0088] A(j:j+L,j:min(j+L,n))=(eye(L+1)−beta*vvp)*A(j:j+L,j:min(j+L,n));

[0089] b1(j:j+L)=(eye(L+1)−beta*vvp)*b1(j:j+L);

[0090] end

[0091] R=triu(A); %ensure no finite precision effects prevent triangularity

[0092] x=R(1:n,1:n)\b1(1:n); %back subst.—C code can take advantage of zeros return

[0093] The above function returns the linear least squares solution. The following code fragment implements decision feedback in the back-substitution step.

[0094] %block DFE “trick” in backsubstitution

[0095] xhat=zeros(n,1);

[0096] for j=n:−1:1

[0097] newx=(b1(j)−R(j,:)*xhat)/R(j,j);

[0098] xhat(j)=slicer(newx,constellation);

[0099] end

[0100] x=xhat;

[0101] This code fragment calls the function “slicer”, which returns the nearest point in the constellation of possible transmitted symbols. The set of all possible transmitted symbols is stored in the vector “constellation.”

[0102] function [value]=slicer(in,constellation)

[0103] % arbitrary slicer [value]=slicer(in,constellation)

[0104] % returns nearest value in constellation vector

[0105] [C,index]=min(abs(constellation-in));

[0106] value=constellation(index);

[0107] Next, a QR least squares solution using Givens rotations will be described. A Givens rotation differs from a Householder transformation, in that it zeros one coordinate of a vector, rather than all but the first. Thus, the QR factorization requires more Givens rotations than Householder reflections. Since the present application involves a sparse matrix, however, Givens rotations may be more efficient.

[0108] A Givens rotation is a matrix that differs from the identity matrix in four coordinates: $G\left(i,k,\theta \right)=\left[\begin{array}{ccccccc}1& & & & & & 0\\ & ⋰& & & & & \\ & & c& \dots & s& & \\ & & ⋮& ⋰& ⋮& & \\ & & -s& \dots & c& & \\ & & & & & ⋰& \\ 0& & & & & & 1\end{array}\right]\begin{array}{c}\\ \\ i\\ \\ k\\ \\ \end{array}$ $ik$

[0109] By appropriate choices of (real) c and (complex) s values, it can be ensured that component, y_k, of y=G^*Tx is zero. Specifically, $y_j=\{\begin{array}{cc}{\mathrm{cx}}_i-{\mathrm{sx}}_k& j=i\\ {\mathrm{sx}}_i+{\mathrm{cx}}_k& j=k\\ x_j& j\ne i,k\end{array}$

[0110] is obtained y_kcan be forced to zero by setting $c=\frac{x_i}{\left(1\angle x_i\right)\sqrt{x_i^2+x_k^2}}\mathrm{and}s=\frac{-x_k}{\left(1\angle x_i\right)\sqrt{x_i^2+x_k^2}}.$

[0111] Given inputs x_i(variable x1) and x_k(variable x2), the following MATLABfunction computes c and s.

[0112] function [c,s]=givens_c(x1,x2)

[0113] % compute c and s for a complex Givens rotation

[0114] if (x2==0)

[0115] c=1;

[0116] s=0;

[0117] else

[0118] den=sqrt(x1′*x1+x2′*x2);

[0119] angle_corr=conj(x1/abs(x1)); %rotation to force c to be real

[0120] c=x1*angle_corr/den;

[0121] s=−x2*angle_corr/den;

[0122] end

[0123] return

[0124] A series of Given rotations can also be used to compute the QR factorization of a matrix. The MATLABfunction presented below, which calls the above function to compute individual c's and s's, uses such a series. As with Householder transforms, Q^*TC is computed instead of Q. Again, the banded structure of the M₁matrix is used to greatly reduce the computational burden. The variables A and b correspond to M₁and C, respectively. The variable R corresponds to the matrix of the same name, and Q^*TC is accumulated in the variable qb2.

[0125] The following MATLABfunction implements the aspects discussed above:

[0126] function x=givensls_bandsp_fact(A,b,taps);

[0127] %Givens LS Solution to overdetermined Ax=b

[0128] %lower bw=taps−1,upper bw=0

[0129] %#rows=#cols+taps−1

[0130] L=taps−1;

[0131] [m,n]=size(A);

[0132] %use Givens algorithm for qr decomp

[0133] qb2=b;

[0134] for j=1:n

[0135] for i=j+L:−1:j+1

[0136] [c,s]=givens_c(A(i−1,j),A(i,j));

[0137] A(i−1:i,j:min(j+L,n))=[c s′;−s c]′*A(i−1:i,j:min(j+L,n));

[0138] qb2(i−1:i)=[c s′;−s c]′*qb2(i−1:i); %compute q′*b as method proceeds

[0139] end

[0140] end

[0141] r=triu(A);

[0142] R1=r(1:n,1:n); %already have R1, should already be upper triangular

[0143] c=qb2(1:n); %pick off top part

[0144] %d=qb1(M+1:M+L); %don't need d

[0145] x=R1\c; %R1 is upper triangular

[0146] %C code can take advantage of zeros in R1

[0147] return

[0148] Fast Givens rotations can be used to compute a factorization of a matrix A consisting of a nonsingular matrix M, a diagonal matrix D with real positive values along its diagonal, and an upper triangular matrix T, with the properties M^*TA=T and M^*TM=D. This MTD factorization is related to the QR factorization by Q=MD^−½ and R=D^½T. Fast Givens rotations are called “fast” because they do not involve square roots. Otherwise, the complexity is on the same order as the normal Givens algorithm.

[0149] The function below computes the parameters used for the fast Givens rotations (analogous to c and s for a standard Givens rotation).

[0150] function [alpha,beta,gamma,type]=fastgivens_c(x,d)

[0151] % fast givens rotations for complex x, real d

[0152] if(x(2)˜=0)

[0153] alpha=−(x(1)/x(2))′;

[0154] beta=−alpha′*d(2)/d(1);

[0155] gamma=−alpha′*beta′;

[0156] if (gamma<=1)

[0157] type=1;

[0158] tau=d(1);

[0159] d(1)=(1+gamma)*d(2);

[0160] d(2)=(1+gamma)*tau;

[0161] else

[0162] type=2;

[0163] alpha=1/alpha;

[0164] beta=1/beta;

[0165] gamma=1/gamma;

[0166] d(1)=(1+gamma)*d(1);

[0167] d(2)=(1+gamma)*d(2);

[0168] end

[0169] else

[0170] type=2;

[0171] alpha=0;

[0172] beta=0;

[0173] gamma=0;

[0174] end

[0175] return

[0176] The function below computes the MTD factorization of the M₁matrix using fast Givens transforms. The variable A contains M₁on entry and T on exit from the routine. The variable L is the number of channel taps minus 1.

[0177] function [M,A,d]=fastgivensqr_bandsp(A,L)

[0178] %fast Givens M(D{circumflex over ( )}−1)T factorization

[0179] %for a banded matrix, upper bandwidth=0; lower bandwidth=L

[0180] %#rows=#cols+L

[0181] [m,n]=size(A);

[0182] for(k=1:m)

[0183] d(k)=1;

[0184] end

[0185] M=eye(m);

[0186] for(j=1:n)

[0187] for(k=j+L:−1:j+1)

[0188] old_d=d(k−1:k);

[0189] [alpha,beta,gamma,type]=fastgivens_c(A(k−1:k,j),old_d);

[0190] if (type==1)

[0191] A(k−1:k,j:min(n,j+L))=[beta 1;1 alpha]′*A(k−1:k,j:min(n,j+L));

[0192] M(k−1:k,1:j+L)=[beta 1;1 alpha]′*M(k−1:k,1:j+L);

[0193] d(k−1)=old_d(2)*(1+gamma);

[0194] d(k)=old _d(1)*(1+gamma);

[0195] else

[0196] A(k−1:k,j:min(n,j+L))=[1 alpha;beta 1]′*A(k−1:k,j:min(n,j+L));

[0197] M(k−1:k,1:j+L)=[1 alpha;beta 1]′*M(k−1:k,1:j+L);

[0198] d(k−1)=old_d(1)*(1+gamma);

[0199] d(k)=old_d(2)*(1+gamma);

[0200] end

[0201] end

[0202] end

[0203] M=M′;

[0204] return;

[0205] A code fragment using this function to compute the least squares solution is listed below. As with the earlier methods, a more efficient implementation would compute M^*TC, rather than M. For example, a method can be used that has the same updates as for M in the routine above, that starts with C instead of I, and that does not transpose at the end.

[0206] [Mm,T,dv]=fastgivensqr_bandsp(M1_est,L);

[0207] S1=T(1:M,1:M);

[0208] Mtb=Mm′*C;

[0209] c=Mtb(1:M);

[0210] Bhat=S1\c; %back substitution

[0211] As with the other functions presented herein, the code can be modified to be more memory efficient.

[0212] Any of the equalizers discussed above can alternatively be a T/k-fractionally-spaced equalizer taking more than one sample per transmitted symbol interval (where “k” indicates the number of samples taken per transmitted symbol, and “T” indicates the duration of one symbol interval). In this case, the channel is also estimated to a time resolution of T/k, using one of several well-known techniques.

[0213] FIG. 7 depicts the sampling of a portion of a received signal 700 performed in an embodiment using a T/2 spaced fractional sampling scheme. The first sample taken of each transmitted symbol is represented by r₀702, r₁704, etc., and the second sample taken of each transmitted symbol is represented as r_0.5706, r_1.5708, etc. (The subscript indicates the number of symbol intervals, T, from the first sample, r₀.) Analogously, the T/2 spaced channel impulse response estimate may be separated into the first sample per symbol period and the second sample per symbol period.

[0214] FIG. 8 depicts a channel impulse response estimate with a time resolution of T/2. The first coefficient of the channel impulse response in a symbol interval is represented by f₀802, f₁804, etc., and the second coefficient of the channel impulse response in a symbol interval is represented as g₀806, g₁808, etc., where the subscripts indicate the symbol interval.

[0215] FIG. 9 depicts a representation of an M matrix 900 for the T/2 spaced fractional sampling scheme where the length of the channel impulse response is two symbol intervals. This is the set of convolution equations expressed in matrix form. The first column 902 of the M matrix 900 becomes the M₂matrix, the second through fourth columns 904 become the M₁matrix and the fifth column 906 becomes the M₃matrix in the equations presented above.

[0216] The M₁matrix is thus modified by interlacing the channel estimates for the first sampling time (f) and the channel estimates for the second sampling time (g). (If more than two samples per symbol (k) are taken, then the channel estimates for these sampling times are interlaced as well.) The r′ vector is modified to include all T/2 spaced samples of the received signal. As in the above-described embodiments, the probe induced ISI is calculated using the M₂and M₃matrices and subtracted from r′ to yield C. The remainder of the approach proceeds exactly as above, only using the new M₁, M₂, and M₃matrices and the new r′ (and hence C) vectors.

[0217] Further, implementation advantages may be obtained by de-interlacing the columns of M and of r′ by grouping the rows associated with each sampling time within the symbol interval. For example, the convolution matrix equation represented in FIG. 9 would become the matrix equation of FIG. 10. FIG. 10 depicts a representation of an M matrix 1000 for the T/2 spaced fractional sampling scheme. The first column 1002 of the M matrix 1000 becomes the M₂matrix, the second through fourth columns 1004 become the M₁matrix and the fifth column 1006 is again the M₃matrix Since each row of the M matrix and the associated row of the r′ vector represent one arithmetic equation, this rearrangement is merely a reordering of the equations.

[0218] As with the T-spaced equalizer, permuting the rows of the M₁matrix and the rows of the C vector results in a new matrix equation mathematically identical to the original. Also, permuting the columns of the M₁matrix and the rows of the {circumflex over (B)} vector results in a new matrix equation mathematically identical to the original. Such permutations may simplify the implementation of the QR factorization algorithm.

[0219] Further, the present invention can be used in conjunction with an approach wherein changing channel characteristics occurring over the duration of the data block are taken into account. Such an approach is disclosed in a U.S. patent application by Joseph T. Graf and Thomas L. Tapp, filed on even date herewith, and entitled “Block Decision Directed Equalization Method and Apparatus.” The entire specification, including any and all drawing figures and appendices, of the U.S. patent application by Joseph T. Graf and Thomas L. Tapp, filed on even date herewith, and entitled “Block Decision Directed Equalization Method and Apparatus”, is hereby incorporated herein by this reference.

[0220] By way of further example, a back substitution approach can be applied to the present invention in conjunction with the teachings of the incorporated patent application. In this embodiment, both upper and lower triangular solutions can be used at the same time. For example, a lower triangular back substitution can proceed from one of the known probe segments to some mid-point of the unknown data segment as an upper triangular back substitution proceeds from the other known probe segment to the same point. Although this approach requires that two QR factorizations be performed, it can still result in overall improved performance. The upper and lower triangular approach can be beneficial, for example, in circumstances where burst noise can corrupt symbols adjacent one side of the unknown data segment and thereby corrupt the channel estimate that is derived for the solution from that side. For instance, when the right-hand decision vector has been corrupted, the symbol decisions proximate to the right-hand side can be based on the solution that is proceeding from the left.

[0221] Further, a QR or combined QR and QR factorization can be applied in an iterative approach. For example, a QR factorization can be used to solve for one or several of the transmitted symbols at one end of a data block. These detected symbols can be used to update the channel impulse response, by, for example, the LMS algorithm. These detected symbols can then be treated as part of the known data (“probe”), and the process is repeated with a now shorter unknown data block.

[0222] Similarly, a QR factorization can be used to solve for one or several of the transmitted symbols at one end of a data block, and a QL factorization can be used to solve for one or several of the transmitted symbols at the other end of a data block. These detected symbols can be used to update the channel impulse response, by, for example, the LMS algorithm. These detected symbols can then be treated as part of the known data (“probe”), and the process is repeated with a now shorter unknown data block. While these iterative approaches involve several QR or QL factorizations, they provide better estimates of the channel impulse response within the block of unknown data.

[0223] It is thought that the present invention will be understood from the description provided throughout this specification and the appended claims, and that it will be apparent that various changes may be made in the form, construct steps and arrangement of the parts and steps thereof, without departing from the spirit and scope of the invention or sacrificing all of their material advantages. The forms herein described are merely exemplary embodiments thereof.