Chapter 3 Probabilistic Data Association Based Detectors
3.3 Simulation Results
Perfect Channel State Information (CSI) Perfect noise variance estimation
Number of subcarrier 64 Length of cyclic prefix 16
Channel Rayleigh Fading
Path 2
Relative power (dB) (0,0)
Modulation QPSK, 16QAM, 64QAM
Table 3.1 Simulation parameters for comparing PDA with GPDA.
In this Section, we use several computer simulation examples to show the performance and the computational cost of the PDA detector and the GPDA detector. We also compare the PDA detector and the GPDA detector with the V-BLAST ZF OSIC [9] and the optimal ML detector in the examples. The simulation parameters are shown in Table 3.1.
In Fig. 3.2, we compare the BER performance for the GPDA detector and the PDA detector with . The result shows that the BER performance of the GPDA detector has about 1 dB improvements as compared with the PDA detector. Since the dimension of the GPDA detector is two times as the dimension of the PDA detector making the PDF of the interference and noise closer to Gaussian distribution. In Fig. 3.3, it shows the performance of the aforementioned detectors for the case
q=4
q=16.
In Fig. 3.4, we compare the effect of sorting for the GPDA detector and the PDA detector with . The result shows that the BER performance of the PDA and the GPDA detectors has approximately more 2 dB than that of the unsorted PDA and the unsorted GPDA detectors. Thus, we can identify that sorting is an important factor for the GPDA and the PDA detectors.
q=4
In Fig. 3.5, we compare the number of iterations for the GPDA detector and the PDA
detector with . The result shows that both the GPDA detector and the PDA detector just
need 2-4 iterations to converge. In Fig. 3.6, it shows the BER performance of the aforementioned detectors for the case
q=4
q =16.
In Fig. 3.7, we compare the complexity for the GPDA detector and the PDA detector with . For the case of complexity of the system, the measurement was calculated using FLOPS function in MATLAB [14], which counts the approximated floating point operations
that the algorithm needs to complete decoding in one block of transmitted symbols. The result shows that there is a great gap of the original complexity between the GPDA detector and the PDA detector. After using the matrix speed-up (speed-up I) tactic, the gap between the GPDA detector and the PDA detector has reduced. If we use matrix speed-up tactic and successive skip tactic simultaneously (speed-up II), the complexity may reduce once more. Overall, the complexity of the GPDA detector is slight more than that of the PDA detector for
q=4
q=4.
In Fig. 3.8, we compare the complexity for the GPDA detector and the PDA detector with . As the figure suggests that, after using speed-up tactic, the complexity of the GPDA detector will be significantly less than that of the PDA detector.
q=16
In Fig. 3.9, we compare the complexity for the GPDA detector and the PDA detector with different modulation order, at SNR=0 dB. We can observe that with greater modulation order, the gap between the GPDA detector and the PDA detector gets wider.
Fig. 3.10 shows the noise variance estimation error for the GPDA detector and the PDA
detector with , at SNR=25 dB. The result shows that the GPDA detector and the PDA detector are almost free from the impact of the noise variance estimated error. Note that,
q =16
E[
2 2
[ ] [ ]
Eσ =Eσ + Δσ2] where E[ ]σ2 is estimated noise variance composed of real noise
variance E[σ2] and estimation error E[Δσ2].
After above comparisons, the GPDA detector shows improved BER performance over the PDA detector by reducing the number of probabilities associated with each transmit symbol. Moreover, the complexity of the GPDA detector is much less than that of the PDA detector, especially in high order modulation. Hence, that is the reason why we choose the GPDA detector rather than the PDA detector.
Fig. 3.2 The BER performance for PDA and GPDA with q=4.
Fig. 3.3 The BER performance for PDA and GPDA with q=16.
Fig. 3.4 The effect of sorting for PDA and GPDA with q=4.
Fig. 3.5 The number of iterations for PDA and GPDA with q=4.
Fig. 3.6 The number of iterations for PDA and GPDA with q=16.
Fig. 3.7 The complexity for PDA and GPDA with q=4.
Fig. 3.8 The complexity for PDA and GPDA with q=16.
Fig. 3.9 The complexity for PDA and GPDA with different modulation order.
Fig. 3.10 The noise variance estimation error for GPDA and PDA.
Chapter 4
Data Detection in MIMO-OFDM System Based on GPDA Detector
4.1 GPDA-MCPDA Detector
The basic idea of the GPDA-MCPDA detector is using the GPDA detector at the low SNR regions, and using parallel MCPDA method to generate numbers of random samples at
the high SNR regions. After generating samples, we will pick up a sample from the final iteration of parallels, and the sample which has minimum distance (i.e. arg min 2
a X∈ i r - Ha ).
Thus, we can get a solution from the GPDA-MCPDA detector.
4.1.1 Markov Chain Monte Carlo Method
A major limitation towards more widespread implementation of Bayesian approaches is that obtaining the posterior distribution often requires the integration of high-dimensional functions. This can be computationally very difficult, but several approaches short of direct integration have been proposed. The MCMC methods [6], which attempt to simulate direct draws from some complex distribution of interest. MCMC approaches are so-named because one uses the previous to randomly generate the next sample value, generating a Markov Chain (as the transition probabilities between sample values are only a function of the most recent
sample value).
The realization in the early 1990’s that one particular MCMC method, the Gibbs sampler, is widely applied to a broad class of Bayesian problems has sparked a major increase in the application of Bayesian analysis.
4.1.1.1 Monte Carlo Integration
The original Monte Carlo approach was a method developed by physicists to use random number generation to compute integrals. Suppose we wish to compute a complex integral
b ( )
This is referred to as Monte Carlo integration.
Monte Carlo integration can be used to approximate posterior (or marginal posterior) distributions required for a Bayesian analysis. Consider the integral ( )I y =
∫
f y x p x d( | ) ( ) x , which we approximate by1
It was observed in the preceding Section 4.1.1.1 that the integral can be approximate by Monte Carlo integration. However, not every density p x( ) can be drawn directly. Now, we suppose the density q x( ) roughly approximates the density p x( ), then
This forms the basis for the method of importance sampling, with
1
4.1.1.3 Introduction to Markov Chains
Before introducing the Gibbs sampler, a few introductory comments on Markov Chains are in order. Let denote the value of a random variable at time , and let the state space refer to the range of possible
Xt t
X values. The random variable is a Markov process if the
transition probabilities between different values in the state space depend only on the random variable’s current state, i.e.,
1 0 1
Pr(Xt+ =s Xj| =sk, ,Xt =si)=Pr(Xt+ =s Xj| t =si) (4.7)
Thus for a Markov random variable the only information about the past needed to predict the future is the current stage of the random variable, knowledge of the values of earlier states do not change the transition probability. A Markov chain refers to a sequence of random variables generated by a Markov process. A particular chain is defined most critically by its transition probabilities, vector of the state space probabilities at step . We start the chain by specifying a starting vector . Often all the elements of are zero except for a single element of 1, corresponding to the process starting in that particular state. As the chain progresses, the probability values get spread out over the possible state space.
t
π ( ) t
t
)
π (0) π (0
The probability that the chain has state value si at time (or step) is given by the Chapman-Kolomogrov equation, which sum over the probability of being in a particular state at the current step and the transition probability from that state into state
t+1
si,
1
Successive iteration of the Chapman-Kolomogrov equation describes the evolution of the chain.
We can more compactly write the Chapman-Kolomogrov equation in matrix form as follows. Define the probability transition matrix P as the matrix whose i j, th element is
( , )
P i j , the probability of moving from state to state i j , P i( → j) . The
Chapman-Kolomogrov equation becomes
( 1) t + = ( ) t
π π P
(4.11)Using the matrix form, we immediately see how to quickly iterate the Chapman-Kolomogrov
equation, as
( )t = (t−1) =( (t−2) ) = (t−2) 2
π π P π P P π P (4.12)
Continuing in this fashion show that
( )t = (0) t
π π P (4.13)
Defining the n-step transition probability p as the probability that the process is in state i j( ),n j given that it started in state , n step ago, i.e., i
Finally, a Markov chain is said to be irreducible if there exists a positive integer such that
( ) ,n 0
pi j > for all . That is, all states communicate with each other, as one can always go from any state to any other state (although it may take more than one step). Likewise, a chain is said to be aperiodic when the number of steps required to move between two states (say
i j,
x
and ) is not required to be multiple of some integer. Put another way, the chain is not forced into some cycle of fixed length between certain states.
y
π*
* *
A Markov chain may reach a stationary distribution , where the vector of probabilities of being in any particular given state is independent of the initial condition. The
stationary distribution satisfies
= π π P
*
(4.15) The condition for a stationary distribution is that the chain is irreducible and aperiodic. When a chain is periodic, it can cycle in a deterministic fashion between states and hence never settles down to a stationary distribution.
A sufficient condition for a unique stationary distribution is that the detailed balance
equation holds,
( ) *j ( ) k
P j→k π =P k→ j π (4.16)
If equation (4.16) holds for all , the Markov chain is said to be reversible, and hence equation (4.16) is also called the reversibility condition. Note that this condition implies
, as the
with the last step following since rows sum to one.
4.1.1.4 Gibbs Sampler
One problem with applying Monte Carlo integration is in obtaining samples from some complex probability distribution. Attempts to solve this problem are the roots of MCMC methods. In particular, they trace to attempts by mathematical physicists to integrate very complex functions by random sampling, and resulting Metropolis-Hastings sampling [6]. The Gibbs sampler [6] [15] (introduced in the context of image processing by Geman 1984), is a special case of Metropolis-Hastings sampling wherein the random value is always accepted (i.e. α =1). The task remains to specify how to construct a Markov Chain with values converged to the target distribution. The key to the Gibbs sampler is that we only consider the univariate conditional distributions (the distribution when all of the random variables but one is assigned fixed value). Such conditional distributions are far easier to simulate than complex joint distributions and usually have simpler forms. Thus, we simulate n random variables sequentially from the n univariate conditions rather than generating a single n-dimensional vector in a single pass using the full joint distribution.
To introduce the Gibbs sampler, consider a bivariate random variable ( , )x y , and suppose we want to compute one or both marginal, p x( ) and p y( ). The idea behind the sampler is that it is far easier to consider a sequence of conditional distributions, p x y( | ) and
( | )
p y x , than to obtain the marginal by integral of the joint density p x y( , ), e.g.,
( ) ( , )
p x =
∫
p x y dy. The sampler starts with some initial value y0 for y and obtains x0 by generating a random variable from the conditional distribution, p x y( | = y0). We use x0 to generate a new value of from the sampler. To do so, we draw from the conditional distribution based on the valuey1 Repeating this process for k times, then it generates a Gibbs sequence of length k, where a subset of points ( ,x y for j j) 1≤ ≤j m k< is taken as our simulated draws from the full joint
distribution. The first m times of the process, called burn-in period, can make the Markov Chain converge to the distribution that near its stationary one.
When more than two variables are involved, the sampler is extended in the obvious fashion. In particular, the value of the kth variable is drawn from the distribution
where denotes a vector containing all of the variables but k. Thus, during the ith iteration of the sample, to obtain the value of i( )k
For example, if there are four variab y z)
|
Now, we consider the equation (3.2) and we use Gibbs sampler to generate samples of
generate initial samples randomly
sample from ( | , , , , )
i( )m
The first Nb iterations of the loop, called burn-in period, are to let the Markov Chain converge to near its stationary distribution. During the next Ns iterations, the Gibbs sampler generate the Ns samples, i.e.,
t sample can be istribution should be converged
after
las the solution of the distribution. Since the d
b s
N +N times iteration.
There are two problems of Gibbs sampler:
1) How do you choose an initial point? A poor choice of initial point will greatly increase the required burn-in time, and an area of much current research is whether an optimal initial point can be found.
2) How many iterates are needed? This question does not have exact answer, the majority
.1.2 GPDA-MCPDA Detector
at higher values of SNR, some of the transition
tioned answers are obtained from the experience.
4
In [12], the author mentioned that
probabilities in the underlying Markov Chain may become very small. As a result, the Markov Chain may be effectively divided into a number of nearly disjoint chains. The term nearly disjoint here means the transition probabilities that allow movement between the disjoint chains are very low. Therefore, a Gibbs sampler that starts from a random point will remain within the set of points surrounding the initial point and thus may not get a chance of visiting sufficient points to find the global solution. In [13] two solutions for solving this problem were proposed: (i) run a number of parallel Gibbs sampler with different initial points; (ii) while running the Gibbs sampler, we assume a noise variance which is higher than it actually is. These two methods turned out to be effective for low and medium SNRs.
In the parallel MCPDA detector, we will focus on these two methods which men
above to improve the performance of MCMC method. First, we use parallel Gibbs samplers with the initial point generated randomly. Second, we compute covariance according to (3.6) rather than (4.25), so named MCPDA. Since we take the variance of residual interference caused by the random samples into account in the equation (3.6), the covariance will increase.
Furthermore, with few times of iteration, the covariance will be gradually narrowing. This may be regarded as automatic Simulated Annealing method [6]. Finally, we will pick up a
sample from the final iteration of parallels, and the sample which has minimum distance (i.e.
arg min 2
a X∈ i r - Ha ). Thus, we can get a solution from the parallel MCPDA detector.
r 3, we have mentioned that the GPDA detector performs well at th
In Chapte e low SNR
regions, so we can only use the GPDA detector at the low SNR regions; however, with the SNR increasing (exceed M dB), the performance of the GPDA detector will get worse gradually. Thus, we need to use parallel MCPDA method to assist the GPDA detector in order to reach better performance, so named GPDA-MCPDA detector. Moreover, MCPDA is similar to GPDA, we only need to add few blocks, and then the GPDA detector will become the MCPDA detector. Therefore, it may be quite simple to implement. The block diagram of the GPDA-MCPDA detector is shown in Fig. 4.1 and the discrepancy between the GPDA detector and the MCPDA detector is shown in Fig. 4.2.
−1
Fig. 4.1 Block diagram of the GPDA-MCPDA detector.
, ( ) 1
Fig. 4.2 The discrepancy between the PDA detector and the MCPDA detector.
.2 GPDA-SD Detector
not have good performance at the high SNR
.2.1 Sphere Decoding
e lgorithm is a quasi-ML detection technique. It promises to find
4
As we know, the GPDA detector does
regions. In order to solve this problem, we try to find a solution which is better than the GPDA solution. Therefore, we use the Sphere Decoding (SD) algorithm to do this work, so named GPDA-SD, which can attain a better performance and lower complexity in the MIMO-OFDM spatial multiplexing system.
4
Th sphere decoding [3] a
the optimal solution with low computational costs under some conditions. The SD algorithm
was first introduced by Finke and Pohst [16] in the context of the closest point search in lattices but it has become very popular in digital communication literature. Its various applications include lattice codes, CDMA systems, MIMO systems, global positioning system (GPS), etc.
4.2.1.1 Real Sphere Decoding
Fig. 4.3 A sphere of radius and centered at
In communication system, the SD algorithm is used to solve the ML problem as follows:
d
r
.2
ˆML arg min
a X∈ i
=
a r - Ha (4.27) The computational complexity of above exhaustive search
m
method is really high. Therefore, the SD algorithm is brought up to avoid the exhaustive search and search only over the possible a which lie in a certain sphere of radius d around the given vector
r
, thereby reducing the search space and, hence, the required co putations (see Fig. 4.3).It is clear that the closest point inside the sphere will also be the closest point for the who
ny points,
2) oints are inside the sphere? If this requires testing the distance of
The SD algorithm does not really address the first question, but usually uses ZF solution to be le point. However, close scrutiny of this basic ideal leads to two key questions.
1) How do you choose radius d? Clearly, if radius is too large, we obtain too ma
and the search remains exponential in size, whereas if radius is too small, we obtain no points inside the sphere.
How can we tell which p
each point form
r
, then there is no point in SD, as we will still need an exhaustive search.the radius of the sphere. However it does propose an efficient way to answer the second. The basic observation is the following. Although it is difficult to determine the points inside a general NT-dimensional sphere, it is trivial to do so in the one-dimensional case. The reason is that a one dimensional sphere reduces to the endpoints of an interval, and so the desired points will be the integer values that lie in this interval. We can use this observation to go from dimension k to dimension k +1. Suppose that we have determined allk-dimensional points that lie in sphere of radius Then, for any such k -dimensional point, the set of admissible values of the
( k + 1)
th d ensional coordinate th lie in the higher dimensional sphere of the same radius ms an interval.The above means that we can determine all
a d.
im at
d for
points in a sphere of dimension NT and
radius d by successively determining all points in spheres of lower dimensions 1, , NT and the same radius d . Such an algorithm for determining the points in an NT-di
sphere essentially constructs a tree where the branches in the kth level of the tree correspond to the points inside the sphere of radius d and dimension k (see Fig. 4.4).
Moreover, the complexity of such an algorithm will depend on the size of the tree, i.e., on the number of points visited by the algorithm in different dimensions.
2, mensional
at hand. To this Fig. 4.4 Search tree of sphere decoding.
With this brief discussion, we can now be mo e specific about the problemr
end, we shall assume that NR≥NT. Note that, the point Ha lies inside a sphere of radius d centered at
r
if and only if2 2
d ≥ r - Ha (4.28)
In order to break the problem into the sub-problems described above, it is useful to consider
In order to break the problem into the sub-problems described above, it is useful to consider