US7302018B2 - Iterative detection in MIMO systems - Google Patents

Abstract

Iterative detection for a MIMO (multiple-input, multiple-output) wireless communications system uses multiple stages for interference cancellation refinements. Detection is determined by maximizing signal to noise plus interference ratios. Equalization matrix updating uses a matrix-inversion-free recursion. For P transmitter antennas and Q receiver antennas, the soft symbol vector estimates are z=Fr−Gs^hard with F a P×Q IMMSE detector, r the received Q vector, and G a P×P interference canceller with s^hard hard decisions of prior symbols and iterations.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

The following copending application dislcoses related subject matter and has a common assignee: application Ser. No. 10/618,970, filed Jul. 14, 2003.

BACKGROUND OF THE INVENTION

The present invention relates to communication systems, and more particularly to multiple-input multiple-output wireless systems.

Wireless communication systems typically use band-limited channels with time-varying (unknown) distortion and may have multi-users (such as multiple cellphone users within a cell). This leads to intersymbol interference and multi-user interference, and requires interference-resistant detection for systems which are interference limited. Interference-limited systems include multi-antenna systems with multi-stream or space-time coding which have spatial interference, multi-tone systems, TDMA systems having frequency selective channels with long impulse responses leading to intersymbol interference, CDMA systems with multi-user interference arising from loss of orthogonality of spreading codes, high data rate CDMA which in addition to multi-user interference also has intersymbol interference.

Interference-resistant detectors commonly invoke one of three types of equalization to combat the interference: maximum likelihood sequence estimation, (adaptive) linear filtering, and decision-feedback equalization. However, maximum likelihood sequence estimation has problems including impractically large computation complexity for systems with multiple transmit antennas and multiple receive antennas. Linear filtering equalization, such as linear zero-forcing and linear minimum square error equalization, has low computational complexity but has relatively poor performance due to excessive noise enhancement. And decision-feedback (iterative) detectors, such as iterative zero-forcing and iterative minimum mean square error, have moderate computational complexity but only moderate performance.

An iterated decision approach appears in Chan and Wornell, A Class of Block-Iterative Equalizers for Intersymbol Interference Channels: Fixed Channel Results, 49 IEEE Tran. Comm. 1966 (2001) and Chan and Wornell, A Class of Asymptotically Optimum Iterated-Decision Multiuser Detectors, Proc. IEEE ICASSP 2265 (2001).

SUMMARY OF THE INVENTION

The present invention provides a method and detector for multiple-input multiple-output (MIMO) systems with progressive refinement of soft symbol estimates by iterative detection across antennas in several stages while accounting for decision-feedback error.

This has advantages including improved performance for wireless systems.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a flow diagram.

FIGS. 2 a-2 c illustrate functional blocks of detectors, receivers, and transmitters.

FIGS. 3 a-3 b simulation results.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

1. Overview

Preferred embodiment detection methods for multi-input, multi-output (MIMO) systems progressively refine soft symbol estimates by iterative detections across multiple antennas in multiple stages while accounting for the effects of decision-feedback errors. FIG. 1 is a flow diagram for first preferred embodiment methods and shows an inner loop which iteratively detects each transmitted symbol using cancellation of previously detected symbols and an outer loop of detection stages which refine the inner loop detections. Explicitly, the preferred embodiments estimate a P-vector of transmitted symbols s from a P-vector z=Fr−Gŝ where F is a P×Q IMMSE detector applied to the Q-vector of received signals r, and G is a P×P interference canceller applied to the P-vector of prior detection hard decisions ŝ; both F and G are refined for each symbol detection and also for successive stages of detections.

The usual iterative MMSE detection method estimates each symbol once using full cancellation of previously detected symbols as if the prior symbol detection were error free. In contrast, preferred embodiment detectors and methods incorporate multistage iterative MMSE detection with the following characteristics:

• (a) N≧2 detection stages. Each detection stage includes P iterations where P is the number of transmitter antennas and symbol streams. Each stage iteratively applies MMSE (IMMSE) detection to obtain soft and hard estimates of all P symbols (s₁, s₂, . . . , s_P). After all P symbols are detected, the next stage is started using the results of the prior stage as initial conditions, and the first stage is analogous to the usual IMMSE and has 0 initial conditions. That is, if (n,i) denotes the ith detection iteration during the nth stage so that ŝ^(n,i) are the current hard decisions for the detected symbols, then ŝ^(n,0)=ŝ^(n−1,P) are initial conditions for incrementing stages and ŝ^(1,0)=0 are the initial conditions for the first stage.
• (b) Account for decision feedback (interference cancellation) error. Determine the effect of decision feedback error by the correlation between the hard decision feedback and the actual symbols. This correlation is a function of the post-detection SINR which depends upon the channel state H. This correlation may be updated every iteration.
• (c) Within each stage, order the iterations according to the post-detection SINR of the symbols. The first iteration detects (hard decision on the soft estimate) the symbol with the largest SINR, and a contribution of this symbol subtracted (at least partial cancellation) for the second iteration. The second iteration detects the symbol among the remaining (non-cancelled) symbols with the largest SINR. And so forth.
• (d) Improve the reliability of the hard decision symbol estimates in the feedback (cancellations) from stage to stage. This improves the interference subtractions. The preferred embodiment approach relates the portion of the hard decision symbol estimate subtracted to the reliability of the hard decision.

Preferred embodiment wireless communication systems incorporate preferred embodiment detection methods. FIGS. 2 a-2 b illustrate functional blocks of a preferred embodiment detector and corresponding receiver for a MIMO wireless communications system; and FIG. 2 c shows a multiple-antenna transmitter for such a system. In such systems components, base stations, and mobile users, could each include one or more application specific integrated circuits (ASICs), (programmable) digital signal processors (DSPs), and/or other programmable devices with stored programs for implementation of the preferred embodiment. The base stations and mobile users may also contain analog integrated circuits for amplification of inputs to or outputs from antennas and conversion between analog and digital; and these analog and processor circuits may be integrated on a single die. The stored programs may, for example, be in external or onboard ROM, flash EEPROM, and/or FeRAM. The antennas may be parts of RAKE detectors with multiple fingers for each user's signals. The DSP core could be a TMS320C6xxx or TMS320C5xxx from Texas Instruments.

2. Single Stage Iterative MMSE Detectors

FIG. 2 b illustrates a receiver with an interference-resistant detector as could be used in a wireless communications system with P transmit antennas (P data streams) and Q receive antennas. FIG. 2 c illustrates a corresponding transmitter with P transmit antennas. The received signal in such a system can be written as:
r=Hs+w
where r is the Q-vector of samples of the received baseband signal (complex numbers) corresponding to a transmission time n:

$r=\left[\begin{array}{c}{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">r</figure-callout>}}_1\left(n\right)\\ {\mathrm{<figure-callout\; id="2"\; label="r"\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">r</figure-callout>}}_2\left(n\right)\\ ⋮\\ r_Q\left(n\right)\end{array}\right];$
s is the P-vector of transmitted symbols (complex numbers of a symbol constellation) for time n:

$s=\left[\begin{array}{c}{s}_1\left(n\right)\\ s_2\left(n\right)\\ ⋮\\ s_P\left(n\right)\end{array}\right];$
H is the Q×P channel matrix of attenuations and phase shifts; and w is a Q-vector of samples of received (white) noise. That is, the (q,p)th element of H is the channel (including multipath combining and equalization) from the pth transmit source to the qth receive sink, and the qth element of w is the noise seen at the qth receive sink.

Note that the foregoing relation applies generally to various systems with various interference problems and in which n, r, s, P, and Q have corresponding interpretations. For example:

• (i) High data rate multi-antenna systems such as BLAST (Bell Labs layered space time) or MIMO and multi-stream space-time coding: spatial interference suppression techniques are used in detection.
• (ii) Broadband wireless systems employing OFDM (orthogonal frequency division multiplex) signaling and MIMO techniques for each tone or across tones.
• (iii) TDMA (time division multiple access) systems having frequency-selective channels with long impulse response which causes severe ISI (intersymbol interference). Use equalizers to mitigate ISI.
• (iv) CDMA (code division multiple access) systems having frequency-selective channels which cause MUI (multi-user interference) as a result of the loss of orthogonality between spreading codes. For high data rate CDMA systems such as HSDPA and 1xEV-DV, this problem is more severe due to the presence of ISI. Equalizers and/or interference cancellation may be used to mitigate these impairments.
• (v) Combinations of foregoing.

P is essentially the number of symbols that are jointly detected as they interfere with one another, and Q is simply the number of collected samples at the receiver. Because there are P independent sources, Q must be at least as large as P to separate the P symbols. A detector in a receiver as in FIGS. 2 a-2 b outputs soft estimates z of the transmitted symbols s to a demodulator and decoder.

Presume that different symbols that are transmitted via P different antennas are uncorrelated and may also may utilize different modulation schemes. This implies the P×P matrix of expected symbol correlations, Λ=E[ss^H], is diagonal with entries equal the expected symbol energies (λ_k=E[|s_k|²]):

$\Lambda =\left[\begin{array}{cccc}{\mathrm{<figure-callout\; id="1"\; label="\lambda "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\lambda </figure-callout>}}_1& 0& \cdots & 0\\ 0& {\mathrm{<figure-callout\; id="2"\; label="\lambda "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\lambda </figure-callout>}}_2& \cdots & 0\\ ⋮& ⋮& ⋰& ⋮\\ 0& 0& \cdots & {\lambda }_P\end{array}\right]$

For linear filtering equalization detectors, such as linear zero-forcing (LZF) or linear minimum mean square error (LMMSE), the soft estimates, denoted by P-vector z, derive from the received signal by linear filtering with P×Q matrix F; namely, z=F r. LMMSE detection finds the matrix F by minimizing the mean square error, E[∥z−s∥²]. With perfect estimation of the channel H, the minimizing matrix F is given by:

$\begin{array}{c}F={\left[H^{H}H+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{\Lambda }^{-1}\right]}^{-1}{H}^H\\ =\Lambda {H^H\left[H\Lambda H^H+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{I}_Q\right]}^{-1}\end{array}$
where σ² is the variance per symbol of the additive white noise w and I_Q is the Q×Q identity matrix. Note F has the form of a product of an equalization matrix with H^H which is the matrix of the matched filter for the channel.

A one-stage iterative (decision-feedback) detector for blocks of P symbols has a series of P linear detectors (P iterations) with each linear detector followed by a (hard) decision device and interference subtraction (cancellation). Each of the P linear detectors (iterations) generates both a hard and a soft estimate for one of the P symbols. The hard estimate is used to regenerate the interference arising from the already-estimated symbols which is then subtracted from the received signal, and the difference used for the next linear symbol estimation. This presumes error-free decision feedback. More explicitly, let s denote the P-vector of transmitted symbols to be estimated, ŝ⁽ⁱ⁾ denote the ith iteration output P-vector of hard symbol estimates (first i components equal to the hard estimates ŝ₁, ŝ₂, . . . , ŝ_i of the first i symbols, s₁, s₂, . . . , s_i, and the remaining P−i components each equal to 0, and z⁽ⁱ⁾ denote the ith iteration output P-vector of soft estimates of s₁−ŝ₁, s₂−ŝ₂, . . . , s_i−1−ŝ_i−1, s_i, s_i+1, . . . , s_p. That is, estimates of the transmitted symbols with the already-estimated first i−1 symbols subtracted out. The hard decision for the ith symbol, ŝ_i, arises from application of a hard decision operator on the soft estimate: ŝ_i=D{z_i ⁽ⁱ⁾}. The ith iteration detection is:
z ⁽ⁱ⁾ =F r−F H ŝ ⁽ⁱ⁻¹⁾
where the second term is the soft estimation of the regenerated (propagated by H) hard decision symbol estimates of the prior i−1 iterations. For the first iteration there are no already-estimated symbols, so ŝ⁽⁰⁾=0_P (a P-vector with each component equal to 0). Of course, the ith iteration soft estimates z₁ ⁽ⁱ⁾, z₂ ⁽ⁱ⁾, . . . , z_i−1 ⁽ⁱ⁾ are ideally just estimates of channel noise because the hard estimates would exactly cancel the transmitted symbols. Thus computational simplicity suggests omitting these computations by zeroing-out the corresponding rows (columns) of the matrices. More precisely, take:
z ⁽ⁱ⁾ =F ⁽ⁱ⁾ r−G ⁽ⁱ⁾ ŝ ⁽ⁱ⁻¹⁾
where F⁽ⁱ⁾ and G⁽ⁱ⁾ are the P×Q detection matrix and the P×P interference cancellation matrix for the ith iteration, respectively, defined as:

$F^{\left(i\right)}=\left[\begin{array}{c}{0}_{\left(i-1\right)xQ}\\ {\Phi }^{\left(i\right)}\end{array}\right]$
G ⁽ⁱ⁾ =F ⁽ⁱ⁾ └B _i−10_Q×(P−i+1)┘

where the inversion matrix for IMMSE is

$\begin{array}{c}{\Phi }^{\left(i\right)}={\left(A_i^H{A}_i+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{\Lambda }_i^{-1}\right)}^{-1}{A}_i^H\\ ={\Lambda }_i{A_i^H\left(A_i{\Lambda }_i{A}_i^H+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{I}_Q\right)}^{-1}\end{array}$
with the last P−i+1 and first i−1 symbol portions of the channel matrix H defined in terms of the P column vectors h₁, h₂, . . . , h_P of H as:

• • A_i=[h_i, h_i+1, . . . , h_P] a Q×(P−i+1) matrix, and
  • B_i=[h₁, h₂, . . . , h_i−1] a Q×(i−1) matrix.
    Also, 0_(i−1)×Q is the (i−1)×Q matrix of 0s, 0_Q×(P−i+1) is the Q×(P−i+1) matrix of 0s, and Λ_i is the lower-right (P−i+1)×(P−i+1) diagonal submatrix of Λ and thus the symbol energies of the symbols not-already estimated.

Ordered detection based on the symbol post-detection signal-to-interference-plus-noise ratio (SINR) is often used to reduce the effect of decision feedback error. In particular, let the detection order be π(1), π(2), . . . , π(P) where π( ) is a permutation of the P integers {1, 2, . . . , P}; that is, the first estimated symbol (hard estimate output of the first stage of the iteration) will be ŝ_π(1) and thus also be the corresponding nonzero element of ŝ⁽¹⁾. The maximum SINR of the components of the first soft estimate z⁽¹⁾, which estimates all P symbols, determines π(1). Similarly, the SINRs of the components of z⁽²⁾, which estimates all of the symbols except the cancelled s_π(1), determines π(2), and so forth. That is, the ith iteration estimates symbol s_π(i), and modifying the foregoing to accommodate the ordering is routine but omitted for clarity in notation.

These MMSE detectors are biased in the sense that E[z_k|s_k]−s_k≠0. However, the bias of the MMSE detectors can be removed by applying a scaling factor to the soft outputs. This scaling factor does not affect post-detection SINR, yet results in increased mean square error compared to the regular biased MMSE estimate. While this unbiasing operation does not affect the performance of LMMSE detectors, it improves the performance of IMMSE detectors because the decision device that is used to generate decision feedback assumes unbiased soft output. The unbiasing operation for IMMSE detectors rescales the soft estimates as follows:
{hacek over (z)} _p =z _p/α_p
where

$\begin{array}{c}{\alpha }_p={\lambda }_p{h_p^H\left[A_p{\Lambda }_p{A}_p^H+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{I}_Q\right]}^{-1}{h}_p\\ ={\left({\left[A_p^H{A}_p+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{\Lambda }_p^{-1}\right]}^{-1}{A}_p^H{A}_p\right)}_{1,1}\end{array}$
with {hacek over (z)}_p denoting the soft output after unbiasing and the subscript 1,1 denoting the (1,1) matrix element. For unbiased IMMSE, variance-based and mean-squared-error-based normalizations are equivalent.
3. Multistage Iterative MMSE Detector Preferred Embodiments

FIG. 1 is a flow diagram for first preferred embodiment detection methods which include refinement stages (outer loop) of iterative MMSE detections (inner loop). More explicitly, the methods include the following steps in which the superscript (n,i) indicates values generated during the ith iteration within the nth stage, (n,0) denotes the initial conditions for the first iteration in the nth stage, and thus (1,0) denotes the overall initial conditions. Explanations and definitions follow this listing of the steps:

• (1) compute the P×P equalization matrix Φ^(n,i) from the estimated Q×P transmission channel matrix H and the P×P correlation matrix R^(n,i−1) which is a diagonal matrix with the (k,k)th element equal to the correlation ρ_k ^(n,i−1) between the kth symbol and the hard decision for the kth symbol available from the iteration (n,i−1) detection.
• (2) compute the forward detection P×Q matrix F^(n,i) using H and Φ^(n,i) from (1).
• (3) compute the P×P cancellation matrix G^(n,i) using H, F^(n,i) from (2), and R^(n,i−1).
• (4) compute soft estimate P-vector z^(n,i) using received Q-vector of samples r, F^(n,i) from (2), G^(n,i) from (3), and the P-vector of hard decisions, ŝ^(n,i−1), available from the iteration (n,i−1) detection.
• (5) compute the post-detection SINRs, γ_k ^(n,i), for soft estimates z_k ^(n,i) from (4) for all P−i+1 symbols not previously detected during the iteration (n,1), (n,2), . . . , (n,i−1) detection(s) of the nth stage; this uses Φ^(n,i) from (1) and R^(n,i−1).
• (6) let π(i) denote the subscript k from (5) which labels the largest γ_k ^(n,i), and detect symbol s_π(i) using soft estimate z_π(i) ^(n,i) from (4); this detection may be part of a sequence decoding and/or a hard decision decoding; also provide hard decision ŝ_π(i) ^(n,i)=D(z_π(i) ^(n,i)) for use in interference cancellation in subsequent detections.
• (7) update the P-vector of hard decisions ŝ^(n,i−1) to ŝ^(n,i) by changing the π(i) component from ŝ_π(i) ^(n,i−1) to ŝ_π(i) ^(n,i) from (6) and leaving all other components unchanged.
• (8) compute the correlation ρ_π(i) ^(n,i) and then update correlation matrix R^(n,i−1) to R^(n,i) by changing the π(i) diagonal element from ρ_π(i) ^(n,i−1) to ρ_π(i) ^(n,i) and leaving all other elements unchaged. The computation of ρ_π(i) ^(n,i) may use γ_π(i) ^(n,i) from (5)-(6) and a constellation plus noise model or may be assessed by ensemble averaging.

Steps (1)-(8) are performed P times to complete one IMMSE detection (inner loop). The initial conditions for the inner loop are taken from the results of the prior stage: R^(n,0)=R^(n−1,P) and ŝ^(n,0)=ŝ^(n−1,P). The first stage (overall) initial conditions are simply R^(1,0)=0 and ŝ^(1,0)=0.

FIG. 2 a illustrates the flow of data (solid arrows) and control (broken arrows) in the foregoing steps which, in more detail, proceed as:
z ^(n,i) =F ^(n,i) r−G ^(n,i) ŝ ^(n,i−1) for 1≦n≦N, 1≦i≦P
with ŝ ^(1,0)=0_P
ŝ ^(n,0) =ŝ ^(n−1,P) for 2≦n≦N
Find the F^(n,i) and G^(n,i) by maximizing the post-detection SINRs of the symbols as follows. First, substituting r=H s+w, yields

$\begin{array}{c}{z}^{\left(n,i\right)}=F^{\left(n,i\right)}r-G^{\left(n,i\right)}{\hat{s}}^{\left(n,i-1\right)}\\ =F^{\left(n,i\right)}\left(Hs+w\right)-G^{\left(n,i\right)}{\hat{s}}^{\left(n,i-1\right)}\end{array}$
so impose a unity gain constraint: require the diagonal elements of F^(n,i) H to equal 1. (This is analogous to the unbiasing operation for IMMSE.) Hence, for each (n,i) the soft symbol estimate may be written:
z _p ^(n,i) =s _p+η_p ^(n,i) for 1≦p≦P
where η_p ^(n,i) is the residual interference plus noise term.

Simplify the optimization to find F^(n,i) and G^(n,i) by presuming the correlation of the transmitted symbols and their hard estimates is expressed by a P×P correlation matrix:

${\Lambda }^{-1}E\left[s{\hat{s}}^{\left(n,i\right)H}\right]=R^{\left(n,i\right)}=\left[\begin{array}{cccc}{\rho }_1^{\left(n,i\right)}& 0& \cdots & 0\\ 0& {\rho }_2^{\left(n,i\right)}& \cdots & 0\\ ⋮& ⋮& ⋰& ⋮\\ 0& 0& \cdots & {\rho }_P^{\left(n,i\right)}\end{array}\right]$
Then the post-detection SINR for symbol s_k, denoted γ_k ^(n,i), is expressed as:
γ_k ^(n,i)=λ_k /E[|η _k ^(n,i)|²]
Now determine F^(n,i) and G^(n,i) by maximizing the γ_k ^(n,i); this yields:
G ^(n,i)=(F ^(n,i) H−I _P)R ^(n,i−1)

$F^{\left(n,i\right)}=\left[\begin{array}{ccc}1/\left[h_1^H{\Psi }^{\left(n,i\right)}{h}_1\right]& \cdots & 0\\ ⋮& ⋰& ⋮\\ 0& \cdots & 1/\left[h_P^H{\Psi }^{\left(n,i\right)}{h}_P\right]\end{array}\right]H^H{\Psi }^{\left(n,i\right)}$
where Ψ^(n,i)=[H(I_P−R ^(n,i−1)2)ΛH^H+σ²I_Q]⁻¹ is a Q×Q equalization matrix.

Applying the matrix inversion lemma gives the alternative formulation used in the foregoing listing of the detection steps:

$F^{\left(n,i\right)}=\left[\begin{array}{ccc}1/{\left[{\Phi }^{\left(n,i\right)}{H}^{H}H\right]}_{11}& \cdots & 0\\ ⋮& ⋰& ⋮\\ 0& \cdots & 1/{\left[{\Phi }^{\left(n,i\right)}{H}^{H}H\right]}_{\mathrm{PP}}\end{array}\right]{\Phi }^{\left(n,i\right)}{H}^H$
where Φ^(n,i)=[H^HH(I_P−R^(n,i−1)2)+σ²Λ⁻¹]⁻¹ is a P×P equalization matrix.

The transformation G^(n,i) controls the amount of interference cancellation at every iteration. When hard decisions are highly reliable (i.e., R^(n,i−1)≅I_P), full subtraction is performed; otherwise, partial subtraction according to the correlation.

The SINRs γ_k ^(n,i) expressed in terms of Φ^(n,i) and prior correlations ρ_k ^(n,i−1) are:
γ_k ^(n,i)={1/(1−ρ_k ^(n,i−1)2)}{1/[σ²Λ⁻¹Φ^(n,i)]_kk−1}
This follows from the maximization finding F and G.

The correlations ρ_p ^(n,i) can be assessed by either (i) ensemble averaging estimation or (ii) an analytical model using γ_p ^(n,i).

In particular, for ensemble averaging, a pilot channel or a set of training symbols that undergoes the same MIMO channel can be used. For instance, in WCDMA, the common pilot channel (CPlCh) can be used for this purpose. First, the received signal corresponding to the known training symbols are passed through the detector. Since the training symbols are foreknown, one can estimate the correlations from the ensemble average of the sample correlation between the training symbols and the resulting detected symbols.

In contrast, an analytic model allows direct computation without the delay of ensemble averaging. Indeed, the definition of correlation implies:
ρ_p ^(n,i)=(1/E[|s _p|²])Σ_jk c _j c _k*Prob[ŝ _p ^(n,i) =c _j ,s _p =c _k]
where c_j and c_k are symbols in the symbol constellation (which may depend upon p because the constellations may differ among the symbol streams), and
E[|s _p|²])=Σ_k |c _k|²Prob[s _p =c _k]
Now the hard decision ŝ_p ^(n,i)=D(z_p ^(n,i)) is a maximum likelihood decision (see section 6 for alternatives) for z_p ^(n,i)=s_p+η_p ^(n,i) where η_p ^(n,i) is the residual interference plus noise and may be approximated as (complex) normally distributed with mean 0 and variance λ_p/γ_p ^(n,i). (Exact expressions which avoid this approximation are more complicated with little accuracy improvement.) Thus the correlations ρ_p ^(n,i) can be evaluated as functions of the SINRs γ_p ^(n,i) and the symbol constellations.

As example, for the QPSK symbol constellation with all four symbols equi-likely, the maximum likelihood decision reduces to the signum function applied to the real and imaginary parts and scaled for energy:
Re{ŝ _p ^(n,i)}=sgn(Re{z _p ^(n,i)})√(λ_p/2)
Im{ŝ _p ^(n,i)}=sgn(Im{z _p ^(n,i)})√(λ_p/2).

Similarly for the 16-QAM constellation (scaled to λ_p=1 for notational convenience):

$\begin{array}{ccc}\mathrm{Re}\left\{{\hat{s}}_p^{\left(n,i\right)}\right\}=-3/\surd 10& \mathrm{if}& \mathrm{Re}\left\{z_p^{\left(n,i\right)}\right\}\le -2\surd 10\\ =-1/\surd 10& \mathrm{if}& -2\surd 10<\mathrm{Re}\left\{z_p^{\left(n,i\right)}\right\}\le 0\\ =1/\surd 10& \mathrm{if}& 0<\mathrm{Re}\left\{z_p^{\left(n,i\right)}\right\}\le 2\surd 10\\ =3/\surd 10& \mathrm{if}& 2\surd 10<\mathrm{Re}\left\{z_p^{\left(n,i\right)}\right\}\end{array}$
And analogously for the imaginary parts.

The normal distribution of η_p ^(n,i) implies for the QPSK constellation:
ρ_p=1−2Q(√γ_p)
where Q(x) is the x-tail of the normal distribution (Gaussian Q function). A similar result holds for the 16-QAM constellation:
ρ_p=1−(1/5){3Q(√γ_p/5)+4Q(3√γ_p/5)+3Q(5√γ_p/5)}
Section 5 describes preferred embodiment implementation details.
4. Illustrative Example of Preferred Embodiment Methods

To illustrate the refining of the interference cancellation in the foregoing preferred embodiment methods, consider the low-complexity case of P=Q=2 and a simple channel matrix

$\hspace{1em}\left[\begin{array}{cc}h& h\\ h& ih\end{array}\right]$
where h is an attenuation and phase factor. (A physical case of parallel antenna pairs with an offset could have such a channel matrix.) Now F is of order 1/h and Φ is of order 1/h², so simplify the notation by combining the channel coefficient h, symbol energy λ (take both symbols of equal energy), and the channel noise power σ² into a single parameter: S=σ²/(|h|²λ), and ignore the h scale factors of F and Φ. That is, take:

$H=\left[\begin{array}{cc}1& 1\\ 1& i\end{array}\right],\mathrm{so}{H}^H=\left[\begin{array}{cc}1& 1\\ 1& -i\end{array}\right],\mathrm{and}{H}^{H}H=\left[\begin{array}{cc}2& 1+i\\ 1-i& 2\end{array}\right]$
First compute the equalization matrix:

$\begin{array}{c}{\Phi }^{\left(n,i\right)}={\left[H^{H}H\left(I-R^{\left(n,i-1\right)2}\right)+SI\right]}^{-1}\\ ={\left[\begin{array}{cc}2\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)+S& \left(1+i\right)\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)\\ \left(1-i\right)\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)& 2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+S\end{array}\right]}^{-1}\\ \mathrm{where}{\rho }_k={\rho }_k^{\left(n,i-1\right)}\mathrm{for}k=1,2\\ =\left(1/\Delta \right)\left[\begin{array}{cc}2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+S& -\left(1+i\right)\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)\\ -\left(1-i\right)\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)& 2\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)+S\end{array}\right]\end{array}$
where the determinant is Δ=[2(1−ρ₁ ²)+S][2(1−ρ₂ ²)+S]−2(1−ρ₁ ²)(1−ρ₂ ²).

Next, compute the forward detection: F^(n,i)=D^(n,i)Φ^(n,i)H^H where D^(n,i) is the 2×2 diagonal matrix which makes the diagonal elements of F^(n,i)H equal 1 to avoid biasing the estimates. Of course, F^(n,i)H equals D^(n,i)Φ^(n,i)H^HH, which is used to evaluate D^(n,i) as follows:

$\begin{array}{c}{\Phi }^{\left(n,i\right)}{H}^{H}H=\left(1/\Delta \right)\left[\begin{array}{cc}2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+S& -\left(1+i\right)\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)\\ -\left(1-i\right)\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)& 2\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)+S\end{array}\right]\left[\begin{array}{cc}2& 1+i\\ 1-i& 2\end{array}\right]\\ =\left(1/\Delta \right)\left[\begin{array}{cc}2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+2S& \left(1+i\right)S\\ \left(1-i\right)\mathrm{<figure-callout\; id="2"\; label="S"\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">S</figure-callout>}& 2\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)+2S\end{array}\right]\\ \mathrm{so}{D}^{\left(n,i\right)}=\Delta \left[\begin{array}{cc}1/\left[2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+2S\right]& 0\\ 0& 1/\left[2\left(1-{\mathrm{<figure-callout\; id="1"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_1^2\right)+2S\right]\end{array}\right].\\ \mathrm{Hence},F^{\left(n,i\right)}=D^{\left(n,i\right)}{\Phi }^{\left(n,i\right)}{H}^H\\ =\Delta \left[\begin{array}{cc}1/\left[2\left(1-{\rho }_2^2\right)+2S\right]& 0\\ 0& 1/\left[2\left(1-{\rho }_1^2\right)+2S\right]\end{array}\right]\times \\ \left[\begin{array}{cc}2\left(1-{\rho }_2^2\right)+S& -\left(1+i\right)\left(1-{\rho }_2^2\right)\\ -\left(1-i\right)\left(1-{\rho }_1^2\right)& 2\left(1-{\rho }_1^2\right)+S\end{array}\right]\left[\begin{array}{cc}1& 1\\ 1& -i\end{array}\right]\\ =\Delta \left[\begin{array}{cc}1/\left[2\left(1-{\rho }_2^2\right)+2S\right]& 0\\ 0& 1/\left[2\left(1-{\rho }_1^2\right)+2S\right]\end{array}\right]\\ \left[\begin{array}{cc}\left(1-i\right)\left(1-{\rho }_2^2\right)+S& \left(1+i\right)\left(1-{\rho }_2^2\right)+S\\ \left(1+i\right)\left(1-{\rho }_1^2\right)+S& \left(-1-i\right)\left(1-{\rho }_1^2\right)-iS\end{array}\right]\\ =\Delta \left[\begin{array}{cc}\left[\left(1-i\right)\left(1-{\rho }_2^2\right)+S\right]/& \left[\left(1+i\right)\left(1-{\rho }_2^2\right)+S\right]/\\ \left[2\left(1-{\rho }_2^2\right)+2S\right]& \left[2\left(1-{\rho }_2^2\right)+2S\right]\\ \left[\left(1+i\right)\left(1-{\rho }_1^2\right)+S\right]/& \left[\left(-1-i\right)\left(1-{\rho }_1^2\right)-iS\right]/\\ \left[2\left(1-{\rho }_1^2\right)+2S\right]& \left[2\left(1-{\rho }_1^2\right)+2S\right]\end{array}\right]\end{array}$

Next, compute the cancellation matrix: G^(n,i)=(F^(n,i)H−I)R^(n,i−1). The term F^(n,i)H equals D^(n,i)Φ^(n,i)H^HH and the two factors D^(n,i) and Φ^(n,i)H^HH appear above, thus

$\begin{array}{cc}{F}^{\left(n,i\right)}H=\left[\begin{array}{cc}1/\left[2\left(1-{\rho }_2^2\right)+2S\right]& 0\\ 0& 1/\left[2\left(1-{\rho }_1^2\right)+2S\right]\end{array}\right]& \\ \mathrm{ }\left[\begin{array}{cc}2\left(1-{\rho }_2^2\right)+2S& \left(1+i\right)S\\ \left(1-i\right)S& 2\left(1-{\rho }_1^2\right)+2S\end{array}\right]& \\ =\left[\begin{array}{cc}1& \left(1+i\right)S/\left[2\left(1-{\rho }_2^2\right)+2S\right]\\ \left(1-i\right)S/\left[2\left(1-{\rho }_1^2\right)+2S\right]& 1\end{array}\right]& \\ \mathrm{Hence},G^{\left(n,i\right)}=\left(F^{\left(n,i\right)}H-I\right)R^{\left(n,i-1\right)}& \\ =\left[\begin{array}{cc}0& \left(1+i\right)S/\left[2\left(1-{\rho }_2^2\right)+2S\right]\\ \left(1-i\right)S/\left[2\left(1-{\rho }_1^2\right)+2S\right]& 0\end{array}\right]& \\ \left[\begin{array}{cc}{\rho }_1& 0\\ 0& {\rho }_2\end{array}\right]& \\ =\left[\begin{array}{cc}0& {\rho }_2\left(1+i\right)S/\left[2\left(1-{\rho }_2^2\right)+2S\right]\\ {\rho }_1\left(1-i\right)S/\left[2\left(1-{\rho }_1^2\right)+2S\right]& 0\end{array}\right]& \end{array}$
Now the soft estimates are:
z ^(n,i) =F ^(n,i) r−G ^(n,i) ŝ ^(n,i−1)

$\begin{array}{c}\mathrm{So},z_1^{\left(n,i\right)}=\Delta \left[\left(1-i\right)\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+S\right]/\left[2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+2S\right]{\mathrm{<figure-callout\; id="1"\; label="r"\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">r</figure-callout>}}_1+\\ \Delta \left[\left(1+i\right)\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+S\right]/\left[2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+2S\right]r_2-\\ [{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2\left(1+i\right)S/\left[2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+2S\right]{\hat{s}}_2^{\left(n,i-1\right)}\end{array}$
And similarly for z₂ ^(n,i). The symbol selected for detection is the symbol with greater SINR.

To see the cancellation in the z^(n,i) estimates, substitute r=Hs+w, thus z^(n,i)=F^(n,i)Hs+F^(n,i)w−G^(n,i)ŝ^(n,i−1). And so:
z ₁ ^(n,i) =s ₁+(1+i)S/[2(1−ρ₂ ^(n,i−1)2)+2S](s ₂−ρ₂ ^(n,i−1) ŝ ₂ ^(n,i−1))+[Fw] ₁
where

$\begin{array}{c}{\left[Fw\right]}_1=\Delta \left\{\left[\left(1-i\right)\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+S\right]/\left[2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+2S\right]\right\}{\mathrm{<figure-callout\; id="1"\; label="w"\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00002.png"\; state="\{\{state\}\}">w</figure-callout>}}_1+\\ \Delta \left\{\left[\left(1+i\right)\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+S\right]/\left[2\left(1-{\mathrm{<figure-callout\; id="2"\; label="\rho "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\rho </figure-callout>}}_2^2\right)+2S\right]\right\}{\mathrm{<figure-callout\; id="2"\; label="w"\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">w</figure-callout>}}_2\end{array}$
And analogously for z₂ ^(n,i).

Without loss of generality, consider z₁ ^(n,i); then the s₂ interference term is (with explicit superscripts)
{(1+i)S/[2(1−ρ₂ ^(n,i−1)2)+2S]}(s ₂−ρ₂ ^(n,i− 1)ŝ ₂ ^(n,i−1))
Now as ρ₂ increases from 0 to 1, the s₂ interference term from F increases because the denominator [2(1−ρ₂ ^(n,i−1)2)+2S] decreases, but this is offset by the increasing cancellation from G. The net effect can be quantified by the relation between ρ_k and S. Thus consider the asymptotic case when n becomes large and ρ₂ ^(n,i) and ŝ₂ ^(n,i−1) converge to limiting values independent of n and i, and ŝ₂ (hopefully) equals s₂. In this case the s₂ interference term simplifies to
{S(1−ρ₂)/[(1−ρ₂ ²)+S]}(1+i)s ₂/2.
To evaluate the factor in braces, apply the equations for ρ_k and γ_k. For simplicity take the QPSK symbol constellation, thus:
ρ_k=1−2Q(√γ_k).

$\begin{array}{c}\mathrm{So}1-{\rho }_k=2Q\left(\surd {\gamma }_k\right)\\ =\surd \left(2/\pi {\gamma }_k\right)e^{-\gamma k/2}\left[1-1/\left({\gamma }_k+2\right)+1/\left({\gamma }_k+2\right)\left({\gamma }_k+4\right)-+\dots \right]\\ \cong \surd \left(2/\pi {\gamma }_k\right)e^{-\gamma k/2}\mathrm{if}{\gamma }_k\mathrm{is}\mathrm{large}\\ \mathrm{And}{\gamma }_k=\left(1/\left[1-{\rho }_k^2\right]\right)\left(1/{\left[S\Phi \right]}_{kk}-1\right)\\ =\left\{2\left[\left(1-{\rho }_j^2\right)+S\right]/\left[2\left(1-{\rho }_j^2\right)+S\right]\right\}/S\mathrm{where}j\ne k\end{array}$
Because S is positive and the correlations are between −1 and 1, the factor 2[(1−ρ_j ²)+S]/[2(1−ρ_j ²)+S] is bounded between 1 and 2 with the upper bound tight for ρ_j close to 1. Hence,
1/S≦γ _k≦2/S
Note that 1/S is the received symbol SNR per channel and in this example there are 4 channels (combinations of 2 receive antennas and 2 transmit antennas); whereas, the post-detection SINR is the total SINR after detection and thus includes the combining across 2 receive antennas.

A small S implies a large γ_k and so 1−ρ_k is easily bounded by e^−1/2S and thus very small compared to S. Consequently, the asymptotic s₂ interference term for small S is:

$\left\{S\left(1-{\rho }_2\right)/\left[\left(1-{\mathrm{<figure-callout\; id="2"\; label="p"\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">p</figure-callout>}}_2^2\right)+S\right]\right\}\left(1+i\right)s_2/2\cong \left(1-{\rho }_2\right)\left(1+i\right)s_2/2\sim e^{-1/2S}\left(1+i\right)s_2/2$
Hence the overall effect of the preferred embodiment staged interference cancellation is a decrease of the s₂ interference from size on the order of S towards an asymptotic size bounded by e^−1/2S, and so the channel noise term dominates the detection error.

The asymptotic case with not-small S can be numerically or graphically estimated. Indeed, presume (due to symmetry) that ρ₁=ρ₂ and γ₁=γ₂ and take a channel with SNR=2(3 dB) as a semi-practical illustration. Thus S=1/2, γ=4[2(1−ρ²)+1]/[4(1−ρ²)+1], and ρ=1−2Q(√γ); this yields ρ≅0.93 and γ≅3.4. And the asymptotic s₂ interference term for S=1/2 is:
{S(1−ρ)/[(1−ρ²)+S]}(1+i)s ₂/2≅(0.07)(1+i)s ₂/2

To see the speed of convergence to the asymptotic state, consider the second stage: z₁ ^(2,1) and presume s₁ was detected prior to s₂ in the first stage so that ρ₂ ^(1,1)=0:
z ₁ ^(2,1) =s ₁+(1+i)S/[2(1−ρ₂ ^(2,0)2)+2S](s ₂−ρ₂ ^(2,0) ŝ ₂ ^(2,0))+[Fw] ₁
where ρ₂ ^(2,0)=ρ₂ ^(1,2) from the first stage. As previously, quantify by using the relations (for QPSK):
ρ₂ ^(1,2)=1−2Q(√γ₂ ^(1,2))
γ₂ ^(1,2)={[(1−ρ₁ ^(1,1)2)+S]/[(1−ρ₁ ^(1,1)2)+S/2]}/S
Again there is a bound, 1/S≦γ₂ ^(1,2)≦2/S, and the upper bound is tight for ρ₁ ^(1,1) close to 1. Now ρ^(1,1) comes from the first stage as
ρ₁ ^(1,1)=1−2Q(√γ₁ ^(1,1))
And ρ₁ ^(1,0)=0 implies γ₁ ^(1,1)={[1+S]/[1+S/2]}/S

Thus for small S, γ₁ ^(1,1) is large, ρ₁ ^(1,1) is close to 1, γ₂ ^(1,2) is large, 1−ρ₂ ^(1,2) is very small and roughly bounded by e^−1/2S, and the s₂ interference becomes {(1+i)/2}(s₂−ρ₂ ^(1,2)ŝ₂ ^(2,0)). Thus the second stage is close to full cancellation because ρ₂ ^(1,2) is very close to 1.

Again for not-small S, consider S=1/2; then γ₁ ^(1,1)=12/5, ρ₁ ^(1,1)≅0.88, γ₂ ^(1,2)≅3.0, ρ₂ ^(2,0)=ρ₂ ^(1,2)≅0.92, and the s₂ interference term is:
(1+i)S/[2(1−ρ₂ ^(2,0)2)+2S](s ₂−ρ₂ ^(2,0) ŝ ₂ ^(2,0))≅(0.1)(1+i)s ₂/2
which is roughly 1.5 times the asymptotic result.
5. Preferred Embodiment Implementations

The filtering with F^(n,i) and G^(n,i) can be combined with the symbol selector (see FIG. 2 a) to reduce the number of computations in the filtering process. This reduces P-dimensional to 1-dimensional filtering. That is, evaluate the SINRs, which uses the equalization matrix Φ (or Ψ), to pick the symbol to detect, and then only form the needed portions of F^(n,i) and G^(n,i) and filter.

Two alternative formulations of the derivations of F^(n,i) appear in the foregoing; namely, use of the P×P equalization matrix Φ or the Q×Q equalization matrix Ψ. Because Q≧P (in order to have the rank of H^HH equal to P), the use of Φ would seem to have lower complexity. However, only for Q much greater than P is Φ preferred. Rather, for Q=P, which MIMO systems typically use, preferred embodiment implementations avoid apparent matrix inversions using the Ψ approach as follows. First, notice that both Φ and Ψ change every iteration, and thus they need to be updated NP times for each complete detection. This seemingly implies NP matrix inversions. However, updating R^(n,i−1) to R^(n,i) changes only the (π(i),π(i)) diagonal element from ρ_π(i) ^(n,i−1) to ρ_π(i) ^(n,i) where π(i) denotes the symbol s_π(i) selected for detection during the ith iteration. Hence,

$\begin{array}{c}{\Psi }^{\left(n,i+1\right)}={\left[H\left(I_P-R^{\left(n,i\right)2}\right)\Lambda H^H+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{I}_Q\right]}^{-1}\\ ={\left[H\left(I_P-R^{\left(n,i-1\right)2}\right)\Lambda H^H+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{I}_Q+H\left(-\Delta R^{\left(n,i\right)2}\right)\Lambda H^H\right]}^{-1}\end{array}$
where the change from R^(n,i−1)2 to R^(n,i)2 has been separated out as the matrix ΔR^(n,i)2 which has a single (π(i),π(i)) diagonal element ρ_π(i) ^(n,i)2−ρ_π(i) ^{(n,i−1) 2}. Further, the first terms on the right side are the inverse of Ψ^(n,i), so
Ψ^(n,i+1)=[[Ψ^(n,i)]⁻¹+H(−ΔR^(n,i)2)ΛH^H]⁻¹
And the single nonzero element structure of ΔR^(n,i)2 allows the matrix products of the second term on the right to be expressed as a tensor product of two vectors:

$\begin{array}{c}H\left(-\Delta R^{\left(n,i\right)2}\right)\Lambda H^H=\left({{\rho }_{\pi \left(i\right)}}^{\left(n,i-1\right)2}-{{\rho }_{\pi \left(i\right)}}^{\left(n,i\right)2}\right){\lambda }_{\pi \left(i\right)}{h}_{\pi \left(i\right)}{h}_{\pi \left(i\right)}^H\\ ={\mathrm{\Delta \rho }}_{\pi \left(i\right)}{h}_{\pi \left(i\right)}{h}_{\pi \left(i\right)}^H\end{array}$
where (ρ_π(i) ^(n,i−1)2−ρ_π(i) ^(n,i)2)λ_π(i) has been abbreviated as Δρ_π(i). Thus
Ψ^(n,i+1)=[[Ψ^(n,i)]⁻¹+Δρ_π(i)h_π(i)h_π(i) ^H]⁻¹
Now applying the matrix inversion lemma yields:
Ψ^(n,i+1)=Ψ^(n,i)−Ψ^(n,i) h _π(i) [h _π(i) ^HΨ^(n,i) h _π(i)+1/Δρ_π(i)]⁻¹ h _π(i) ^HΨ^(n,i)
which has converted the Q×Q matrix inversion into a 1×1 matrix inversion (i.e., just numerical reciprocal) due to the structure of h_π(i)h_π(i) ^H as a tensor product of Q×1 matrices (vectors). In short, h_π(i) ^HΨ^(n,i)h_π(i)+1/Δρ_π(i) is a scalar. Consequently, the recursion:
Ψ^(n,i+1)=Ψ^(n,i)−Ψ^(n,i) h _π(i) h _π(i) ^HΨ^(n,i) /[h _π(i) ^HΨ^(n,i) h _π(i)+1/Δρ_π(i)]
has no matrix inversions beyond the initial Ψ^(1,0)=[HΛH^H+σ²I_Q]⁻¹.

And the SINR can also be expressed in terms of Ψ^(n,i) as
γ_p ^(n,i+1)=λ_p h _p ^H(Ψ^(n,i)−Ψ^(n,i) h _π(i) h _π(i) ^HΨ^(n,i) /D)h _p
where D=λ_p ⁻¹(1−ρ_π(i) ^(n,i)2)⁻¹−h_p ^HΨ^(n,i)h _p. Hence, the foregoing preferred embodiment could be cast in terms of Ψ^(n,i) with the foregoing updating.

The analogous analysis for Φ^(n,i) follows similar lines:

$\begin{array}{c}{\Phi }^{\left(n,i+1\right)}={\left[H^{H}H\left(I_P-R^{\left(n,i\right)2}\right)+{\mathrm{<figure-callout\; id="2"\; label="\sigma "\; filenames="US07302018-20071127-D00000.png,US07302018-20071127-D00003.png"\; state="\{\{state\}\}">\sigma </figure-callout>}}^2{\Lambda }^{-1}\right]}^{-1}\\ ={\left[{\left[{\Phi }^{\left(n,i\right)}\right]}^{-1}+H^{H}H\Delta R^{\left(n,i\right)2}\right]}^{-1}\end{array}$
and again the single nonzero element structure of ΔR^(n,i)2 allows the matrix product H^HH ΔR^(n,i)2 to be expressed as a constant times the tensor product of two vectors: Δρ_π(i)(H^HH)_π(i) 1_π(i) ^H where (H^HH)_π(i) is the π(i) column P-vector of H^HH and 1_π(i) is the column P-vector with all components equal to 0 except the π(i)th component. Note that here Δρ_π(i)=ρ_π(i) ^(n,i−1)2−ρ_π(i) ^(n,i)2).

Then the matrix inversion lemma yields the recursion:
Φ^(n,i+1)=Φ^(n,i)−Φ^(n,i)(H ^H H) _π(i)1_π(i) ^HΦ^(n,i)/[1_π(i) ^HΦ^(n,i)(H ^H H)_π(i)+1/Δρ_π(i)]
and again there are no matrix inversions beyond the initial [H^HH+σ²Λ⁻¹]⁻¹. Further, this updating of Φ^(n,i) could be used in the SINR computations.

Note that the foregoing updating method without matrix inversions would also apply to methods in which changes occur in more than one value of the correlation coefficients, such as in various iterated-decision LMMSE methods. Multiple single value recursions can be performed to achieve to update multiple changes, where each recursion is used to update one change.

6. Decision Alternatives

The foregoing preferred embodiments hard decision method D (which converts z_π(i) ^(n,i) into ŝ_π(i) ^(n,i)) uses a maximum likelihood approach that applies a hard-limiter (hard-clipping) non-linearity. Alternative preferred embodiment non-linearities, such as hyperbolic tangent and soft linear clipping, may result in better performance. They can be described mathematically as follows for BPSK:

• • Hyperbolic tangent: D(z)=tanh(Re{z}/ξ)
  • Soft linear clipped:

$\begin{array}{cc}\mathrm{??}\left(z\right)=+1& \mathrm{if}\mathrm{Re}\left\{z\right\}>\xi \\ =-1& \mathrm{if}\mathrm{Re}\left\{z\right\}<-\xi \\ =\mathrm{Re}\left\{z\right\}& \mathrm{if}\mathrm{Re}\left\{z\right\}\le \xi \end{array}$
where ξ is a non-negative constant. For higher order modulation, the decision is done in the bit level following the bit log likelihood ratio computation (as described in Section 7). Alternatively, for QPSK, the real and imaginary parts of: D(z) can be taken as tanh(Re{z}/ξ) and tanh(Im{z}/ξ), respectively; or the corresponding soft linear with Re{z} and Im{z}, respectively. Further, for higher order (larger constellations) the soft linear extends by selecting the constant to make D(z)=z for z in the constellation.

In these cases the decision statistics that are close to zero are effectively scaled instead of clipped; the scaling depends upon the variance of the residual interference plus noise. Use the same cancellation (interference subtraction) with these hard decisions defining the ŝ_k ^(n,i).

7. Decoding

The demodulator converts the soft symbol estimates z_π(1) ^(N,1), z_π(2) ^(n,2), . . . , z_π(P) ^(N,P) output by an N-stage detector into conditional probabilities; and the conditional probabilities translate into bit-level log likelihood ratios (LLRs) for (sequence) decoding. In more detail, the LLRs in terms of the bits u_pk which define the constellation symbol s_p (e.g., two bits for a QPSK symbol, four bits for a 16QAM symbol, etc.) are defined (dropping the superscripts on the z_p) as

$\begin{array}{c}\mathrm{LLR}\left(u_{\mathrm{pk}}\right)=\log \left\{P\left[u_{\mathrm{pk}}=1z_p]/P[u_{\mathrm{pk}}=0z_p\right]\right\}\\ =\log \left\{P\left[u_{\mathrm{pk}}=1z_p]\right\}-\log \left\{P[u_{\mathrm{pk}}=0z_p\right]\right\}\end{array}$
The LLRs can be computed using a channel model. For example,

$\begin{array}{c}\mathrm{LLR}\left(u_{\mathrm{pk}}\right)=\log \left\{P\left[u_{\mathrm{pk}}=1z_p]/P[u_{\mathrm{pk}}=0z_p\right]\right\}\\ =\log \left\{p\left(z_pu_{\mathrm{pk}}=1)/p(z_pu_{\mathrm{pk}}=0\right)\right\}+\\ \log \left\{P\left[u_{\mathrm{pk}}=1\right]/P\left[u_{\mathrm{nk}}=0\right]\right\}\end{array}$
where the first log term includes the probability distribution of the demodulated symbol z_p which can be computed using the channel model. The second log term is the log of the ratio of a priori probabilities of the bit values and typically equals 0. So for an AWGN channel where the residual interference (interference which is not cancelled) is also a zero-mean, normally-distributed, independent random variable, the channel model gives:
p(z _p |s _p =c)˜exp(−|z _p −c| ²/γ_p)
where c is a symbol in the symbol constellation and γ_p is a normalization typically derived from the channel characteristics and the detector type. Of course, γ_p is just twice the variance of the estimation error random variable.

Then compute LLRs by using an approximation which allows direct application of the channel model. Take p(z_p|u_pk=1)=p(z_p|s_p=c_pk=1) where c_pk=1 is the symbol in the sub-constellation of symbols with kth bit equal 1 and which is the closest to z_p; that is, c_pk=1 minimizes |z_p−c_k=1|² for c_k=1 a symbol in the sub-constellation with kth bit equal to 1. Analogously for p(z_p|u_pk=0) using the sub-constellation of symbols with kth bit equal 0. Then with equal a priori probabilities of the bit values and the notation subscript k=1 and k=0 indicating symbols with kth bit 1 and 0, respectively, the approximation yields

$\begin{array}{c}\mathrm{LLR}\left(u_{\mathrm{pk}}\right)=\log \left\{p\left(z_pu_{\mathrm{pk}}=1)/p(z_pu_{\mathrm{pk}}=0\right)\right\}\\ \cong 1/{\gamma }_p\left\{{\min }_{k=0}{z_p-c_{k=0}}^2-{\min }_{k=1}{z_p-c_{k=1}}^2\right\}\end{array}$
Thus the LLR computation just searches over the two symbol sub-constellations for the minima.

The LLRs are used in decoders for error correcting codes such as Turbo codes (e.g., iterative interleaved MAP decoders with BCJR or SOVA algorithm using LLRs for each MAP) and convolutional codes (e.g. Viterbi decoders). Such decoders require soft bit statistics (in terms of LLR) from the detector to achieve their maximum performance (hard bit statistics with Hamming instead of Euclidean metrics can also be used but result in approximately 3 dB loss). Alternatively, direct symbol decoding with LLRs as the conditional probability minus the a priori probability could be used

8. Simulation Results

Simulations for 4×4 MIMO QPSK and 16QAM systems in uncorrelated (IID) MIMO channels were performed to compare a preferred embodiment with other detection methods. In particular, the comparison included a preferred embodiment method, a linear MMSE, an iterative MMSE with optimal SINR-based detection ordering (similar to the first stage of a preferred embodiment with correlations either 0 or 1), and an iterated-decision linear MMSE as in Chen and Wornell (see background) using five stages to achieve convergence.

FIGS. 3 a-3 b illustrate the raw BER versus E_b/N₀ results. The following table expresses the performance of a 2-stage preferred embodiment relative to the 5-stage iterated-decision linear MMSE and the iterative MMSE with optimal SINR detection in terms of E_b/N₀ (dB) required to achieve a BER of 0.01:

	Method	(4,4) QPSK	(4,4) 16QAM
	Iterative MMSE/optimal SINR	3.3	10/0
	5-stage iterative decision	4.2	10.7
	2-stage preferred embodiment	3.1	9.5

The preferred embodiments show considerable gain with 2 stages relative to the iterative MMSE with optimal SINR which is roughly comparable to a preferred embodiment with only 1 stage. The implementation preferred embodiments of section 5 allow inclusion of additional stages with small complexity increase. Note that the implementations of section 5 may also be applied to updates in the iterative-decision methods.
9. Modifications

The preferred embodiments may be varied while retaining one or more of the features of stages of detections with each detection using optimal SINR and cancellation with correlation-weighted current hard decisions.

For example, a sub-optimal alternative to the correlation weighting is to choose a set of pre-determined weighting coefficients, which do not depend on the channel, but increase as the number of iterations and/or stages increase. Also, a stage may be terminated prior to completion and the then-current decoding used, so some symbol estimates would be from the current stage and some from the prior stage.

Claims (20)

1. A multistage iterative minimum mean square error detection method for a multiple-input, multiple-output (MIMO) communications system, comprising the steps of:

(a) receiving Q samples at Q receive antennas of P symbols transmitted from P transmit antennas of a remote transmitter, where P is an integer greater than 1 and Q is an integer not less than P;

(b) iteratively first detecting each of said P symbols using said Q samples wherein each of said first detections is a minimum mean square error detection including cancellation of symbols detected in prior iterations of said first detections; and

(c) after step (b), iteratively second detecting ones of said P symbols using said Q samples wherein each of said second detections is a minimum mean square error detection including cancellation of said symbols detected in either prior iterations of said second detections or in said first detections.

2. The method of claim 1, wherein:

(a) said second detections have the form F⁽ⁱ⁾ r−G⁽ⁱ⁾ ŝ⁽ⁱ⁻¹⁾ where i is an iteration index, r is a Q-vector of received samples, F⁽ⁱ⁾ is a normalized product of an equalization matrix and a matched filter matrix, G⁽ⁱ⁾ is a cancellation matrix, and ŝ⁽ⁱ⁻¹⁾ is a P-vector of hard decisions on said symbols detected in either prior iterations of said second detections or in said first detections.

3. The method of claim 2, wherein:

(a) said F⁽ⁱ⁾ and G⁽ⁱ⁾ are found by maximizing the signal-to-interference-plus-noise ratios of said symbols.

4. The method of claim 2, wherein:

(a) said matrix F⁽ⁱ⁾ has the form DΦ⁽ⁱ⁾ H^H where D is a diagonal matrix normalizing the diagonal terms of F⁽ⁱ⁾, Φ⁽ⁱ⁾ is a P×P equalization matrix, and H^H is the Hermitian conjugate of the Q×P channel matrix for the channel from the P sources to the Q sinks.

5. The method of claim 4, wherein:

(a) said P×P equalization matrix Φⁱ is found by a matrix-inversion-free recursion using Φ^(i−l).

6. The method of claim 2, wherein:

(a) the matrix F⁽ⁱ⁾ has the form D H^H Ψⁱ where D is a diagonal matrix normalizing the diagonal terms of F⁽ⁱ⁾, Ψ⁽ⁱ⁾ is a Q×Q equalization matrix, and H^H is the Hermitian conjugate of the Q×P channel matrix for the channel from the P sources to the Q sinks.

7. The method of claim 6, wherein:

(a) said P×P equalization matrix Ψ⁽ⁱ⁾ is found by a matrix-inversion-free recursion using Ψ⁽ⁱ⁻¹⁾.

8. The method of claim 3, wherein:

(a) the matrix G⁽ⁱ⁾ has the form G⁽ⁱ⁾=(F⁽ⁱ⁾H-I_F)R⁽ⁱ⁻¹⁾ where H is the Q×P channel matrix for the channel from the P sources to the Q sinks and R⁽ⁱ⁻¹⁾ is a correlation matrix of correlations of symbols and their hard decision estimates available in iteration i−1.

9. The method of claim 8, wherein

(a) the elements of said correlation matrix R⁽ⁱ⁻¹⁾ of claim 8 are found from said signal-to-interference-plus-noise ratio together with said hard decision method and a normal approximation for residual interference plus noise.

10. The method of claim 8, wherein

(a) the elements of said correlation matrix R⁽ⁱ⁻¹⁾ of claim 8 are found from pilot symbols and ensemble averaging.

11. The method of claim 2, wherein:

(a) said hard decision vector ŝ⁽ⁱ⁾ of claim 2 is computed by ŝ_k ⁽ⁱ⁾=D(z_k ⁽ⁱ⁾) where z=F⁽ⁱ⁾ r−G⁽ⁱ⁾ŝ⁽ⁱ⁻¹⁾ and D(z) is defined as the symbol of a constellation which minimizes the distance from D(z) to z.

12. The method of claim 2, wherein:

(a) said hard decision vector ŝ⁽ⁱ⁾of claim 2 is computed by ŝ_k ⁽ⁱ⁾=D(z_k ⁽ⁱ⁾) where z=F⁽ⁱ⁾r−G^{(i) ŝ (i−1)} and D(z) is defined as tanh(Re{z}/c) for symbols z with ±1 values where c is a non-negative constant and an extension on the bit level for symbols represented with multiple bits.

13. The method of claim 1, further comprising:

(a) after step (c) of claim 1, iteratively third detecting each of said P symbols wherein each of said third detections is a minimum mean square error detection and includes cancellation of previously-detected ones of said symbols including results of step (c) of claim 1.

14. The method of claim 1, wherein:

(a) said P symbols are transmitted from P antennas; and

(b) said Q samples are received from Q antennas.

15. An iterative minimum mean square error detection method for a multiple-input, multiple-output (MIMO) communications system, comprising the steps of:

(a) receiving Q samples at Q receive antennas of P symbols transmitted from P transmit antennas of a remote transmitter, where P is an integer greater than 1 and Q is an integer not less than P;

(b) iteratively detecting each of said P symbols wherein each of said detections is a minimum mean square error detection using said Q samples including cancellation of previously-detected ones of said symbols; and

(c) wherein said detections have the form F⁽ⁱ⁾ r−G⁽ⁱ⁾ ŝ⁽ⁱ⁻¹⁾ where i is an iteration index, r is a Q-vector of received samples, F⁽ⁱ⁾is a normalized product of an equalization matrix and a matched filter H^H where H is a Q×P channel matrix, G⁽ⁱ⁾ is a cancellation matrix, and ŝ⁽ⁱ⁻¹⁾ is a P-vector of hard decisions on previously-detected symbols; and

(d) wherein said equalization matrix for iteration i is found by a matrix-inversion-free recursion using said equalization for iteration i−1.

16. The method of claim 15, wherein:

(a) said F⁽ⁱ⁾ and G⁽ⁱ⁾ are found by maximizing the signal-to-interference-plus-noise ratios of the symbols.

17. The method of claim 15, wherein:

(a) the matrix F⁽ⁱ⁾ has the form DΦ⁽ⁱ⁾ H^H where D is a diagonal matrix normalizing the diagonal terms of F⁽ⁱ⁾ and Φ⁽ⁱ⁾ is a P×P equalization matrix.

18. The method of claim 17, wherein:

(a) said P×P equalization matrix Φ⁽ⁱ⁾ is found by a matrix-inversion-free recursion using Φ⁽ⁱ⁻¹⁾ with Φ^(σ)=[H^HH+σ² Λ⁻¹]⁻¹ where σ²is the variance of the channel noise and Λis a diagonal matrix of transmitted symbol energies.

19. The method of claim 15, wherein:

(a) the matrix F⁽ⁱ⁾ has the form D H^H Ψ⁽ⁱ⁾ where D is a diagonal matrix normalizing the diagonal terms of F⁽ⁱ⁾ and Ψⁱ is a Q×Q equalization matrix.

20. The method of claim 19, wherein:

(a) said P×P equalization matrix Ψ⁽ⁱ⁾ is found by a matrix-inversion-free recursion using Ψ⁽ⁱ⁻¹⁾ with Ψ⁽⁰⁾=[H Λ H^H+σ² I_Q]⁻¹ where σ² is the variance of the channel noise and Λ is a diagonal matrix of transmitted symbol energies.

Images