The basic Monte Carlo (MC) technique was examined in Chapter 9 and, in Chapter 10, the Monte Carlo technique was applied to several simple communications systems. In both Chapters 9 and 10 we saw that the Monte Carlo method, while applicable to all systems without regard to architecture or complexity, has one very significant drawback. The fundamental problem with the Monte Carlo method is that the time required to execute a simulation and obtain a reliable estimate of system performance is often very long. In some cases the required run time may be so long that use of the Monte Carlo method is not practical.
The semianalytic (SA) technique was introduced in Chapter 10. While the SA technique yields simulations that execute very rapidly, application of the SA technique is restricted to situations in which the pdf of the sufficient statistic, upon which symbol decisions are based, is known. In many simulations this statistic is not known, and the SA technique cannot be applied.
In this chapter we take a very brief look at simulation techniques aimed at overcoming the lengthy run-time requirements of the Monte Carlo method. The three methods considered here are quite different. The first of these methods, tail extrapolation, involves curve fitting to MC simulation results. The second technique is based on estimating the pdf of a decision metric through the application of moment methods. The third method to be addressed, importance sampling, involves biasing the channel noise in a way that forces more decision errors to be made. Importance sampling is a variance reduction technique that provides an estimate of the BER to be obtained that has a smaller variance than the estimate provided by an MC simulation of equivalent execution time. The application of variance reduction techniques essentially involves a tradeoff between analysis and computer run time. A number of these techniques have been investigated and the search for efficient simulation techniques remains an area of active research.
Of the three techniques presented in this chapter, importance sampling has received the most attention and appears to be the most generally applicable. We therefore consider importance sampling in more detail than the other two simulation methods. The effective application of any of these methods (especially importance sampling) requires considerable analytical skill. This chapter should be considered only a brief introduction to these topics. The interested student should consult the research literature on this subject. A few key references are given at the end of this chapter.
The first technique we consider is tail extrapolation [1, 2]. This method is applied by executing a number of Monte Carlo simulations using Eb/N0 values for which reliable simulation results can be obtained with reasonable run times. These results are then extrapolated to values of Eb/N0 where Monte Carlo simulations are not practical. Of course, one must be extremely careful not to extrapolate into regions of Eb/N0 where extrapolation is not justified. An example for an assumed system is illustrated in Figure 16.1. For A < Eb/N0 < B, reliable values of the error probability can be computed with reasonable run times using Monte Carlo simulation. This yields points 1, 2, 3, and 4. For our system of interest, extrapolation is valid for B < Eb/N0 < C, and therefore point 5, although extrapolated, is valid. Continuing extrapolation gives point 6. However, our example system exhibits a floor on the probability of error and, as a result, extrapolation of the simulated values into the region for which Eb/N0 > C is not valid. Therefore, point 6 is not valid. One must be very careful about extrapolating results into regions where experemental or simulated results are not available.
Suppose we know, or are willing to assume, that the pdf of the decision metric is a Gaussian random variable. For this case the probability of symbol error will be given by (
, where Q(x) is the Gaussian Q-function and g(Eb/N0) is a function of Eb/N0 determined by system parameters such as the modulation format. We see that fitting simulated points to a Gaussian Q-function (such as points 1, 2, 3, and 4 in Figure 16.1) may allow for accurate performance extrapolation to higher values of Eb/N0.
As a more general example, assume that the decision metric is not strictly Gaussian but may be approximated by a pdf of a generalized exponential class.
This class of pdfs is defined by the expression
where Γ(·) denotes the gamma function, m is the mean of X, v is the mode, and σ is a parameter related to the spreading of the pdf about the mean. The parameter σ is, of course, the standard deviation of X for ν = 2, since a Gaussian pdf results for v = 2. A Laplacian pdf results for v = 1. If we assume that the mean is zero, one can show that at for large values of T (error threshold) [1]
The left-hand side of (16.2) is the BER of the system when we use a threshold of T and, therefore we write the left side as 
. Taking the logarithm of (16.2) twice, we see that
or
Now suppose we run M MC simulations using a variety of thresholds, T1,T2...TM. The BER for each of these thresholds is estimated yielding 
. The plot of ln 
as a function of ln 
should be a straight line, determined using linear regression, with y intercept −v ln (σ) and slope v. This is illustrated in Figure 16.2. The mode v is determined from the slope. Once v is known, σ is determined from the y intercept yi according to
Once v and σ have been determined, (16.2) can be used as an estimator of the BER of the system. Tail extrapolation can be extended to other pdfs. This technique will be accurate if the analyst is skilled, or lucky, enough to select an appropriate family of pdfs for a given system.
We know that the error probability of a digital communications system can be written in the form
where T is the error threshold (values of X > T result in errors) and fV(v) is the pdf of the decision statistic. One technique for estimating the BER is to first estimate the pdf and then evaluate the integral defined in (16.6) using numerical integration. An intuitive way of accomplishing this is to form a histogram of the data representing the values of V. In order to be useful the histogram must have a sufficient number of bins for V > T to allow the integration defined by (16.6) to be carried out with the required accuracy. In addition, as we saw in Chapter 8, the histogram is a biased estimator for finite data size N. Fortunately there are alternatives to the histogram method. The two most popular techniques for estimating the pdf from a data set x[n] are the Gram-Charlier series [3, 4] and the Parzen series [4].
The Gram-Charlier series approximates a pdf in the form
where Hk(·) represents the Chebyshev-Hermite polynomials and Ck denotes the series coefficients. In writing (16.7) we have made, for mathematical simplicity, the assumption that X is a zero-mean unit-variance random variable. The Chebyshev-Hermite polynomials are defined by the recursion relationship
in which H0(y) = 1 and H1(y) = y. Given the Chebyshev-Hermite polynomials, the coefficients can be calculated from
Note that even though fY(y) is unknown, (16.9) can be evaluated in terms of the moments of Y. If the random variable of interest Y is not a zero-mean unit-variance variable, we can apply the transformation
where µy and σy represent the mean and standard deviation of Y, respectively. Under this transformation Z becomes a zero-mean unit variance random variable so that (16.7) can be applied directly with Y replaced by Z.
The Gram-Charlier series provides a good approximation to the target pdf in the neighborhood of the mean but usually provides a poor approximation in the tails of the pdf. Unfortunately, for BER estimation it is the tails of the pdf that are of interest. The Gram-Charlier series has a number of other difficulties. The Gram-Charlier approximation is not asymptotically unbiased and does not uniformly converge to the target pdf as more terms are added. Also, for a finite number of terms, N, the approximation defined by (16.7) is not a true pdf and may even be negative for particular values of y. Despite these shortcomings, the Gram-Charlier series is a useful pdf estimator for many applications.
The natural estimator for the moments of a random variable producing a data sequence x[n] is given by
In practice of course, the data size N must be finite. For finite N, this estimator, unfortunately, has no optimality properties. It is, however, usually a consistent estimator and therefore satisfactory performance can be expected for sufficiently large N [5].
The shortcomings of the Gram-Charlier series can at least be partially overcome through the use of the Parzen pdf estimator. The Parzen estimator takes the form
The choices for g(x) and h(N) are somewhat arbitrary but reasonable choices that work well in many applications are [4]
and
In (16.12) N is the size of the data set, the function g(·) is a weighting function, and h(N) is a smoothing factor. Although the Parzen estimator is biased for finite N, it can be shown that if h(N) → 0 as N →∞ and Nh(N) →∞ as N →∞, the Parzen estimator is asymptotically unbiased. Clearly, the smoothing function defined by (16.14) satisfies these requirements. It can also be shown that the Parzen estimator is consistent.
Both the Gram-Charlier and the Parzen techniques have been successfully applied to a number of problems. A recent paper explores the use of both of these techniques for estimating the BER of a digital communications system operating in a variety of environments. These include the AWGN channel, AWGN with cochannel interference, and AWGN plus multipath [6].
As we saw in Chapters 9 and 10, MC simulations perform a random experiment a large number of times, N, and count the number of outcomes, NA, denoting some event, A, of interest. One can estimate the probability of event A using the equation
Assuming that NA denotes the number of errors in the transmission of N symbols, high-performance communication systems yield values of NA many orders of magnitude smaller than N. The brief discussion of confidence intervals in Chapter 9 showed that to generate a reliable estimator, we require NA to be at least 10, and preferably 100. In many situations the requirement for large NA can lead to very long simulation execution times. This is especially true when the BER of a communication system is low. If the BER is small, a simulated demodulation error is a rare (low probability) but important event. The goal of importance sampling is to alter the simulation in a controlled way, in order to increase the number of these important or rare events while still enabling the determination of the correct probability of demodulation error. This requires a change in (16.15). Importance sampling falls into a class of simulation methodologies known as variance reduction techniques. The goal of variance reduction techniques is to develop an unbiased estimator which exhibits a reduced variance and/or a reduced simulation execution time as compared to a Monte Carlo simulation.
Before beginning our discussion of importance sampling, it is important to point out that a detailed description of importance sampling is well beyond the scope of this introductory text. Rather than treating this subject in depth, we instead consider the estimation of the area of a geometrical shape in a way that makes use of the fundamental concepts of importance sampling. This is followed by an example of importance sampling applied to a communications system.
In Chapter 9 the subject of Monte Carlo simulation was introduced by developing an estimator for the area of a two-dimensional shape. (Recall that this led to an estimator for the value of π.) In this chapter we again consider an estimator for the area of a two-dimensional shape (an ellipse in this case) in a way that provides insight into the importance sampling simulation technique. At the end of this section a simple communications system simulation is presented. This introduction hopefully makes the basic literature on importance sampling more understandable to the reader. One wishing to apply importance sampling to a practical communications system will need to become familiar with the basic literature. A suggested list of papers is given at the end of this chapter.
Suppose one wishes to use Monte Carlo simulations to estimate the area of an ellipse in the (x, y) plane. In particular, let us examine the problem of finding the area of an ellipse which has a major axis in the x direction of length 
and a minor axis of length 2.
The specified ellipse contains all points (x, y) such that
The area of the ellipse, Ae, is
By rewriting the integral, we can put the problem in a form better suited for Monte Carlo integration and for demonstrating importance sampling. The first step is to represent the limits of the integral by a function in the integrand. This is accomplished by defining an indicator function, he(x, y), such that
This gives
We next define a bounding area that contains the ellipse and has an area easy to calculate. The area of this bounding area is denoted Abound. Multiplying the integrand by unity in the form Abound/Abound gives
The next step is to define a pair of random variables (X, Y) that are uniformly distributed over Abound. The pdf for this pair of random variables is
Combining (16.20) and (16.21) gives
As in Chapter 9, the estimator for the area of the ellipse is
where N denotes the number of points generated within the bounding area, and Ne denotes the number of points falling inside the ellipse.
One can show that the estimator defined by (16.23) is consistent and unbiased provided the bounding area fully encloses the ellipse. This is illustrated in the following example.
Example 16.1.
As a simple example, N = 500 points in a square box, centered on the origin, are generated. The area of the ellipse is estimated using (16.23) with three different box sizes; Abound = 4, 8, and 100. The upper half of Figure 16.3 shows the sample points generated in this simulation with Abound = 100 (10 units on a side), and the lower half of Figure 16.3 shows the running estimate of the area, along with the theoretically calculated value. It is clear from Figure 16.3 that we need more than 500 samples in this simulation to generate an accurate estimate of the ellipse area.
Next, let us examine the effect of changing the area of the bounding box. Setting Aboound = 8 dramatically improves the convergence properties of the estimator. As a result, estimated area is within a few percent of the correct value after only 100 samples. Based on this observation we are tempted to conclude that the smaller the bounding box, the faster the estimator will converge. This is correct, to a point. Note, however, that when Abound = 4, the estimator quickly converges, but to an incorrect value. The estimator is biased, since the bounding box no longer includes the entire ellipse. If the bounding box were to continue to shrink, the rate of convergence would continue to improve at the cost of increased bias and ultimately we would simply be finding the area of the bounding box and not the area of the ellipse. The convergence properties for Abound = 100, 8, and 4 are illustrated Figure 16.4.
The effect of the bounding box size can be predicted from (16.23). The relative precision of the ratio Ne/N increases with Ne. In other words, we wish to have as many points as possible inside the ellipse, leading us to select small boxes. However, the MC estimator will become biased if the box fails to fully enclose the entire ellipse, leading to an optimal box size of Abound = 8, which is the smallest square that encloses the ellipse. The challenge is to find estimators that both converge quickly and are unbiased.
We saw in the previous section that a two-dimensional MC simulation will generate random points over some region of the (x, y) plane. To make the estimator converge quickly, we would like the region used by the MC estimator to be as small as possible, while still ensuring that the area of interest is fully enclosed. For simplicity, we initially assume that the bounding region is a rectangle. However, there is nothing in the development of the estimation algorithm that requires this. The algorithm is applicable for any bounding region of known area. Consider, for example, a bounding rectangle. The smallest rectangle that includes the ellipse has an area of 
Since this is less than the area of the smallest bounding square, the estimator based on a rectangular bounding region will converge faster than the estimator based on a square bounding region.
We now extend the search for a superior bounding box to higher-order polygons. We can find polygons that enclose the ellipse and have areas even smaller than It will become increasingly difficult to calculate the area of the bounding polygon analytically and will also become increasingly difficult to generate points uniformly distributed within the polygon. There is typically a point of diminishing returns, where the saving in computer execution time does not justify the additional analytical difficulties. The ultimate bounding region is simply the ellipse itself. In this case, Abound = Ae, and every randomly generated point will fall inside the ellipse. Thus, Ne = N. Equation (16.23) indicates that such a simulation will be extremely efficient and will produce the exact answer after a single random point has been generated. Of course, in this limiting case there is no need to perform the simulation, since Ae must to be determined analytically before the simulation code can be developed. This, of course, eliminates the need for the simulation.
One of the dangers of making highly efficient simulations is that a small calculation error may cause part of the ellipse to fall outside the bounding region, generating a biased estimator. If one did not know the correct answer, it would be impossible to detect the bias by simply observing the simulation results.
In summary, when using a uniform pdf in MC simulations, one should:
Determine the area of the bounding region analytically.
Prove analytically that the bounding region fully encloses the area of interest.
Find an algorithm for efficiently generating uniformly distributed points within the bounding region.
Minimize the area of the bounding region, given constraints 1, 2, and 3.
Run the simulation long enough to observe a large number of samples in the ellipse.
These concerns and requirements are quite easy to understand for this simple problem. They will reappear in the following sections in more complicated simulation strategies, and the mathematical complexities will tend to make them more difficult to visualize and appreciate.
To this point we have considered only the basic Monte Carlo algorithm. Recall from Chapter 9 that in Monte Carlo simulations we have two counters, a sample counter and an event counter. The sample counter is incremented each time the underlying random experiment, such as generating a sample and testing to see if it falls in an area of interest, is performed. The event counter is incremented each time an event of interest, such as a sample falling within a given area occurs. Counters are usually incremented by adding the number one to the counter contents. We now show that, by adding a weight wi to the event counter rather than the number one, the necessity for having a uniform distribution across the area of interest is removed. The concept of incrementing by weights rather than by integers is central to the application of importance sampling methodology to BER estimation.
Previously in (16.19) we multiplied the integrand by unity in the form Abound/ Abound to yield (16.20). A more general estimator is formed by multiplying the integrand in (16.19) by a two-dimensional weighting function divided by itself, w(x, y)/w(x, y). As before, we replace 1/w(x, y) by a pdf. Specifically
This gives
This is the more general form of (16.22), where fXY(x, y) can be any pdf that is nonzero within the ellipse. The area estimator corresponding to (16.24) is
where random vectors (xi, yi) are IID (independent and identically distributed) with pdf fXY(x, y). This is a generalization of the previous result, since for the uniform case the reciprocal of the pdf, that is, the weighting function, was the constant Abound for all points inside the ellipse. Equation (16.25) allows us to use nearly any pdf in our simulation and includes the uniform pdf as a special case. The variance and bias of the estimator will be a function of the pdf selected for the simulation.
To avoid estimator bias, the uniform pdf had to include the entire ellipse. In other words, the uniform pdf could not be zero over any finite area inside the ellipse. This same condition applies here. One can show that the area estimator will be unbiased provided fXY(x, y) > 0 for all points inside the ellipse. Since many common pdfs are nonzero over the entire (x, y) plane, it is not difficult to generate unbiased estimators. Minimizing the variance is a more difficult problem since the convergence rate of the simulation is a rather complex function of the pdf. The following example provides insight into this problem.
Example 16.2.
In this example we restrict our attention to Gaussian pdfs, and investigate the convergence rate for three different standard deviations, namely, σ = 10, σ = 1, and σ = 0.2. Figure 16.5 shows the result with σ = 10. (Note the original 5 × 5 box. The ellipse is obscured by the sampling points.) When σ is large, the pdf inside the ellipse is small but nearly constant. Thus, relatively large weight is assigned to each sample falling in the ellipse. Several large jumps followed by intervals with a hyperbolic (1/N) decay rate can be seen. The area estimate slowly converges to the correct value, in much the same way as a uniform pdf and a large bounding area.
We next consider the case in which σ = 1. For this case the standard deviation is close to the dimensions of the ellipse and, as a result, the convergence rate is much more rapid. A large number of samples fall inside the ellipse, and each is assigned a relatively small weight. This case roughly corresponds to the uniform pdf when the bounding area is just slightly larger than the ellipse. This result is shown in Figure 16.6.
Perhaps the most interesting case is when σ is small compared to the dimensions of the ellipse. Nearly all of the points in the sample simulation fall inside the ellipse, most very close to the origin. However, at sample 153, the Gaussian random number generator produces a sample near the edge of the ellipse. Because the pdf for this sample is so small, it is assigned a phenomenally large weight, 1.3874 × 1016 in this case, and the large weight causes the area estimate to suddenly increase by many orders of magnitude. The resulting behavior is shown in Figure 16.7. We call this point an extreme event because it is an extremely rare, but extremely important, event. Since the estimate is unbiased, the extreme event is required to offset the effect an initially low estimate. The only way that the estimator can be unbiased is to run the simulation long enough to see many extreme events. This will typically require far more processor time than any of the methods mentioned to date. This simulation is somewhat similar to the uniform case when the bounding area was smaller than the ellipse. In the uniform sampling case a biased estimate was produced. In this new case in which Gaussian sampling is used the estimator is unbiased, but only because of the presence of extreme events. If one terminates the simulation prior to the occurrence of the first extreme event or shortly after the occurrence of an early (small N) extreme event, very misleading results will be obtained from the simulation. Note that in all three cases the estimator Ne/N considered here is unbiased. However, it may require an unreasonable amount of computer time to average out the extreme events.
The effect of extreme event behavior is familiar to us from past studies. As a simple example, assume that a digital communications system has a BER of 10−6. Also assume that an MC simulation is performed and the first 99 transmitted symbols are received without error. The BER estimator 
= Ne/N will produce zero for 1 ≤ N ≤ 99, a value that is clearly too small. If symbol number 100 in the simulation is received in error (an extreme event), the value of jumps from zero to 10−2, which is a value four orders of magnitude too large. If a long span of correct decisions follows this error, will decrease as 1/N for N > 100. This behavior will continue until the next error occurs. This is essentially the same behavior observed in Figure 16.7, with the only difference being that the error in Figure 16.7 is multiplied by a large weight wi. Note that extreme events occurring early in a simulation run have a more pronounced effect than extreme events occurring later in a simulation run.
All three simulations used the same random number seed. All three estimators are unbiased and consistent, but the mode of convergence varies widely. These modes are summarized in Table 16.1.
In this very simple example, there are a variety of methods for predicting the mode of convergence. The scatter diagrams indicate when the pdf covers too little of the ellipse or too much of the plane. Another clear indicator is the weighting function. Ideally the weighting function should be small over the entire area of the ellipse. If enough data is collected, the estimates and samples of the error functions also clearly indicate the modes. Ideally, one would like to establish confidence intervals for these estimators. This is difficult because the estimates no longer have a binomial distribution. The Gaussian approximation can be used if a sufficiently large number of symbols are processed. The number of samples required for the approximation to be accurate will depend on the convergence mode. When the convergence is rapid, only a few hundred samples may be necessary. The slowly converging estimators are similar to the uniform case, and a few hundred points inside the ellipse are necessary. In the last case, many hundred extreme events would need to be observed before the Gaussian approximation could be applied.
Finding the area of an ellipse is a deterministic problem. The random variables used in the MC simulation are strictly artificial, and do not model any physical device or signal. Most physical communication systems are subject to random perturbations and process random data. The first inclination is to develop a simulation that uses the same noise source pdf as the physical system. However, as we will show, this is not necessary. The simulation can use a biased noise source, having a different pdf than the one used in the physical system. For example, suppose a physical system generates a jointly Gaussian pair of random variables (x, y) having the pdf
which σx and σy are the standard deviations of X and Y, respectively, and r represents the correlation coefficient. The probability that a pair of random variables falls inside the ellipse described earlier is
This expression is difficult to evaluate. A simple MC solution to this problem would be to use the sum
where the random vectors (xi, yi) are IID with pdf fphy(x, y). Suppose we did not have a random number generator that would produce fphy(x, y), but we did have a random number generator to produce samples with the distribution fsim(x, y). We can still solve the problem using Monte Carlo techniques. As in the previous section, we begin by multiplying the integrand of (16.27) by unity in the form of w(x, y)/w(x, y). Unlike the earlier development, the weighting function is now defined as
Substituting (16.29) into (16.27) gives
which can be approximated by the sum
The random vectors (xi, yi) are IID with pdf fsim(x, y). As with the previous problem, one must be concerned about the bias, consistency, and convergence mode of this estimator. It is easy to show that when fsim(x, y) > 0 for all points where he(x, y)fphy(x, y) > 0, the estimate will be unbiased. One can also show that this ensures that the estimate is consistent. As in the simpler problem, the mode of convergence may be slow, fast, or very slow because of the existence of exteme events. The mode will depend on the weighting function w(x, y) = fphy(x, y)/fsim(x, y) inside the ellipse. Convergence rates are summarized in Table 16.2.
Before applying the ellipse area estimation problem to a communications system, we need to add one final twist. The channel generates noise, but this noise may be altered by the receiver into a new distribution with an unknown pdf. We need to find the probability that samples generated by this new pdf fall in the region of interest. Working only with two-dimensional pdfs for now, assume that the channel generates the physical noise, (x, y), with pdf fphy(x, y). These random variables are fed to a receiver, which applies operator g(x, y) = (α, β). This new set of random variables has pdf frec(α, β). The probability the random variables (α, β) falls into some region of interest, such as the ellipse, is
The complexity of g may make it impractical to calculate frec(α, β) from fphy(x, y). The MC simulation of this problem will be
where (αi,βi) has pdf frec(α, β). It is possible to run this simulation even though we do not know how to calculate frec(α, β). If we generate random vectors (xi, yi) with pdf fphy(x, y) and then pass each sample vector through function g, the resulting vector (αi,βi) will have pdf frec(α, β). A more descriptive way of writing (16.33) is
Suppose we now alter fphy(x, y) to fbias(x, y). This will cause frec(α, β) to change to frec−bias(α, β). To remove the effects of this bias from the MC estimate, we need to calculate a weighting function
where (αi,βi) = g(xi, yi). This weighting function can be used to alter (16.34) to create
which can be shown to be a consistent and unbiased estimate.
The geometric problem of finding the area of an ellipse can be mapped to the problem of finding the BER of a communication system. The two sets of random variables (x, y) and (α, β) represent the channel noise waveform and the decision metric, respectively. In the communications system problem, the channel noise often has far more than two dimensions, even though the decision metric may have only one dimension. These two pairs of random variables were related by operator g(·), which represents the receiver. The ellipse of the geometric problem now corresponds to the values of the decision metric which causes the receiver to have a demodulation error. Finding the BER of the system is equivalent to finding the probability of a randomly generated sample falling inside the ellipse in the geometric problem. We will bias the channel noise in order to increase the frequency of error events, and calculate a weighting function
The BER estimator is
In the geometric problem, we ensured an unbiased estimator by making sure that the bounding region included the entire ellipse. The equivalent requirement here is simply that all values of the decision metric can be generated. It is not difficult to satisfy this requirement. The rate of convergence in the geometric problem depends on how we select the variance of the noise source that establishes the sampling points. If a large number of samples fall outside the ellipse, the convergence rate is slow. If we generate almost no samples outside the ellipse, the convergence is again slow because of extreme events. We ideally desire a mix of samples both inside and outside the ellipse. In the communication system problem, we ideally want a mix of decision metric values both in the error region and in the error-free region. Too many samples in either region will cause the convergence rate to be very slow. Figure 16.8 illustrates the system with the additional functions required to implement importance sampling represented by heavy lines.
There are a variety of methods for biasing the noise when applying importance sampling. The two the most common methods are called Conventional Importance Sampling (CIS) and Improved Importance Sampling (IIS). In CIS one simply increases the variance of the channel noise, which is equivalent to operating the system at a lower signal-to-noise ratio (SNR). In IIS the mean of the noise is altered rather than the variance. For a comparison of these two techniques and examples illustrating their application to communications systems, the interested student is referred to the literature. (See the Further Reading section.)
Example 16.3.
In this example we consider CIS applied to the differential QPSK system previously considered in Chapter 10. The MATLAB script for running a CIS simulation follows:
% File: c16_CISQPSK.m
Eb = 20:2:32; No = -50; % Eb and Noin dBm
ChannelAttenuation = 70; % channel attenuation in dB
EbNodB = (Eb-ChannelAttenuation)-No; % Eb/No in dB
EbNo = 10.(EbNodB./10); % Eb/No in linear units
BER_T = 0.5*erfc(sqrt(EbNo)); % BER (theoretical)
N = ones(size(BER_T))*2000; % set N=2000 for all runs
CISBias = 1+(EbNo/20); % set CIS bias
BER_CIS = zeros(size(Eb)); % initialize BER vector
for k=1:length(Eb) % main loop starts here
BER_CIS(k) = c16_CISQPSKrun(N(k),Eb(k),...
No,ChannelAttenuation,0,0,0,0,CISBias(k));
disp(['Simulation ',num2str(k*100/length(Eb)),'% Complete']);
end
semilogy(EbNodB,BER CIS,'o',EbNodB,2*BER_T,'-')
xlabel('Eb/No (dB)'), ylabel('BER'), grid;
% End of script file.This script sets the bias for the CIS simulation and then calls the function listed in Appendix A. This function is nearly identical to the Monte Carlo simulation of the QPSK system considered in Chapter 10. The only difference is that the channel noise is biased by increasing the noise variance.
The BER estimator in the example is defined by (16.38). A sample output from this simulation is shown in Figure 16.9. The circles in this figure represent the CIS estimates of the BER, while the solid curve is the theoretical performance. At first glance there appears to be some significant errors in the measurements at high SNR. However, these results are based only on 1,000 demodulated symbols. Using conventional MC simulations it would be impractical to measure any BER less than approximately 10−3. However, from the results of the CIS simulation we are able to get some idea of the BER for error rates approaching 10−8.
This chapter focused on alternatives to the Monte Carlo simulation method for the determination of the BER in a digital communications system. The search for alternative simulation methodologies is motivated by the lengthy execution times often encountered with Monte Carlo methods, especially when the BER is small. Three techniques were briefly discussed. These were tail extrapolation, pdf estimation, and importance sampling. The successful application of any of these techniques requires considerable analytical skill. The material presented in this chapter should be considered only a high-level overview and the student wishing to pursue an in-depth understanding of these methods will wish to consult the literature.
The topics covered in this chapter are also covered in
M. C. Jeruchim, P. Balaban, and K. S. Shanmugan, Simulation of Communication Systems, 2nd ed., New York: Kluwer Academic/Plenum Press, 2000.
A very large number of papers have been written on the subject of importance sampling. A partial reading list consists of the following:
P. Smith, M. Shafi, and G. Hongshen, “Quick Simulation: A Review of Importance Sampling Techniques in Communications Systems,” IEEE Journal on Selected Areas in Communications, Vol. 15, No. 4, May 1997, pp. 597–613.
M. C. Jeruchim, “Techniques for Estimating the Bit Error Rate in the Simulation of Digital Communication Systems,” IEEE Journal on Selected Areas Communications, Vol. 2, No. 1, January 1984, pp. 153–170.
R. J. Wolfe, M. C. Jeruchim, and P. M. Hahn, “On Optimum and Suboptimum Biasing Procedures for Importance Sampling in Communication Simulation,” IEEE Transactions on Communications, Vol. 38, No. 5, May 1990, pp. 639–647.
K. S. Shanmugan and P. Balaban, “A Modified Monte Carlo Simulation Technique for the Evaluation of Error Rate in Digital Communication Systems,” IEEE Transactions on Communications, Vol. 28, No. 11, November 1980, pp. 1916–1924.
M. C. Jeruchim, “On the Application of Importance Sampling to the Simulation of Digital Satellite Multihop Links,” IEEE Transactions on Communications, Vol. 32, No. 10, October 1984, pp. 1088–1092.
B. R. Davis, “An Improved Importance Sampling Method for Digital Communication System Simulations,” IEEE Transactions on Communications, Vol. 34, No. 7, July 1986, pp. 715–719.
P. M. Hahn and M. C. Jeruchim, “Development in the Theory and Application of Importance Sampling,” IEEE Transactions on Communications, Vol. 35, No. 7, July 1987, pp. 706–714.
Q. Wang and V. Bhargava, “On the Application of Importance Sampling to BER Estimation in the Simulation of Digital Communication Systems,” IEEE Transactions on Communications, Vol. 35, No. 11, November 1987, pp. 1231–1233.
D. Lu and K. Yao, “Improved Importance Sampling Technique for Efficient Simulation of Digital Communication Systems,” IEEE Journal on Selected Areas in Communications, Vol. 6, No. 1, January 1988, pp. 67–75.
M. C. Jeruchim, P. M. Han, K. P. Smyntek, and R. T. Ray, “An Experimental Investigation of Conventional and Efficient Importance Sampling,” IEEE Transactions on Communications, Vol. 37, No. 6, June 1989, pp. 578–587.
D. Remondo, R. Srinivasan, V. Nicola, and W. Van Etten, “Adaptive Importance Sampling for Performance Evaluation and Parameter Optimization of Communication Systems,” IEEE Transactions on Communications, Vol. 48, No. 4, April 2000, pp. 557–565.
J. Porath and T. Aulin, “Improved Technique for Quick Error Rate Estimation of Multi-Dimensional Communication Schemes,” IEE Proceedings—Communications, Vol. 146, No. 6, December 1999, pp. 343–346.
16.1 | Assume that the conditional error probability, conditioned on the transmission of a binary 0, for a binary digital communications system can be expressed  where we assume that m is negative. Show that PE|0 can be expresssed in terms of a incomplete gamma function Γ(a, x) where  By choosing m appropriately, and assuming that PE|0 = PE|0, plot the system error probability for v = 1.9, 2, and 2.1. (Note: MATLAB contains a routine for computing the incomplete gamma function.) |
16.2 | Using the recursion relationship for the Chebyshev-Hermite polynomials, evaluate Hk(y) for 2 ≤ k ≤ 6. Using these results determine the Gram-Charlier series coefficients Ck for 0 ≤ k ≤ 6. |
16.3 | A three-dimensional surface is defined by the equation  find the volume enclosed by this surface using Monte Carlo simulation. For each iteration of the MC simulation, generate three independent, identically distributed, uniform random numbers, ranging from -10 to +10. How many points did you need to generate before you became 90% confident that your error was less than 1%? |
16.4 | Repeat Problem 16.3, but allow the x, y, and z points to cover different ranges. How many points did you need to generate before you became 90% confident that your error was less than 1%? |
16.5 | Repeat Problem 16.3, but use a three-dimensional, jointly Gaussian pdf for x, y, and z. What are the optimal values for mean and variance for each of the three marginal pdfs? |
16.6 | Suppose you wish to estimate the area of a two-dimensional ellipse, using Monte Carlo simulation and a jointly Gaussian pdf. Would you always select uncorrelated random variables, or would there be times when it would help to have them correlated? |
16.7 | Water dripping from a faucet hits a the bottom of a sink with a jointly Gaussian pdf (zero mean, standard deviation 0.1 in both directions, correlation coefficient is zero). There is a circular drain in the sink, centered at the origin, of radius 0.1. Using Monte Carlo simulation, estimate what fraction of the drops hit the drain. |
16.8 | Repeat Problem 16.7, but in the simulation use a jointly Gaussian pdf that has a standard deviation of 0.2 in both directions. The physical problem has not changed—you will be using a different pdf in the simulation than you observe in the physical world. |
16.9 | Write a simulation that implements an MSK transmitter, AWGN channel, and matched filter receiver. Perform a Monte Carlo simulation for this system to determine the SNR required to achieve a BER of 10−2. |
16.10 | Repeat Problem 16.9, but use importance sampling. Bias the noise by altering it’s variance. What is the optimal amount to alter the variance, and how much execution time does this save you? |
16.11 | Repeat Problem 16.10, but bias the noise by adding a constant. What is the best constant to add and how much simulation time does this save you? |
16.12 | Write a simulation that implements an OQPSK transmitter. Place a bandpass filter on the transmitter output to limit the transmitted signal to the main lobe. Create two adjacent channels, which are also OQPSK transmitters limited to the main lobe. Place the adjacent channels as close as possible to the desired channel, without allowing the spectra to overlap. Pass the signals through an AWGN channel, and a matched filter receiver. If the interfering signals are not present, set the BER of the system to 10−3. Use MC simulation to estimate the BER when the interfering signals are equal in power to the desired signal. |
Appendix A: MATLAB Code for Example 16.3
% file: c16_CISQPSKrun.m
function BER_CIS=CISQPSKrun(N,Eb,No,ChanAtt,...
TimingBias,TimingJitter,PhaseBias,PhaseJitter,CISBias)
fs = 1e+6; % sampling rate (samples/second)
SymRate = 1e+5; % symbol rate (symbols/second)
Ts = 1/fs; % sampling period
TSym = 1/SymRate; % symbol period
SampPerSym=fs/SymRate; % samples per symbol
SymToSend = N; % symbols to be transmitted
ChanBW = 4.99e+5; % bandwidth of channel (Hz)
CISWeightIntegrator = 1; % importance sampling weight
CISWeightIntegratorOld = 1; % importance sampling weight
MeanCarrierPhaseError = PhaseBias; % mean of carrier phase
StdCarrierPhaseError = PhaseJitter; % std dev of phase error
MeanSymbolSyncError = TimingBias; % mean symbol sync error
StdSymbolSyncError = TimingJitter; % std dev symbol sync error
ChanGain = 10^(-ChanAtt/20); % channel gain (linear units)
TxBitClock = Ts/2; % Tx clock period
RxBitClock = Ts/2; % Rx clock period
TxSymSent = 1; RxSymDemod = 0; % Tx and Rx symbol counters
%
RxNoiseStd = sqrt((10^((No-30)/10))*(fs/2)); % std dev of noise
TxSigAmp = sqrt(10^((Eb-30)/10)*SymRate); % signal amplitude
probe1 = zeros((SymToSend+1)*SampPerSym,1); % probe 1 memory
probe2 = zeros((SymToSend+1)*SampPerSym,1); % probe 2 memory
probe1counter = 1; probe2counter = 1; % initialize probes
%
% Buffers that contain the transmitted and received data.
%
[unused,SourceBitsI] = random_binary(SymToSend,1);
[unused,SourceBitsQ] = random_binary(SymToSend,1);
%
% Differentially encode the transmitted data.
%
TxBitsI = SourceBitsI*0; TxBitsQ = SourceBitsQ*0;
for k=2:length(TxBitsI)
TxBitsI(k) = or(and(not(xor(SourceBitsI(k),SourceBitsQ(k))),...
xor(SourceBitsI(k),TxBitsI(k-1))), ...
and(xor(SourceBitsI(k),SourceBitsQ(k)),...
xor(SourceBitsQ(k),TxBitsQ(k-1))));
TxBitsQ(k) = or(and(not(xor(SourceBitsI(k),SourceBitsQ(k))),...
xor(SourceBitsQ(k),TxBitsQ(k-1))), ...
and(xor(SourceBitsI(k),SourceBitsQ(k)),...
xor(SourceBitsI(k),TxBitsI(k-1))));
end;
%
% Make a complex data stream of the I and Q bits.
%
TxBits = ((TxBitsI*2)-1)+(sqrt(-1)*((TxBitsQ*2)-1));
%
% Initialize transmitter and the receiver integrate and dump filter.
%
RxIntegrator = 0; TxBitClock = 2*TSym;
%
% Design the channel filter and state array if needed.
%
[b,a] = butter(2,ChanBW/(fs/2));
b = [1]; a = [1]; % bypass filter
[junk,FilterState] = filter(b,a,0);
%
% Loop once for each sample.
%
while TxSymSent < SymToSend
%
% Update transmitter clock. Get new data bits if required.
%
TxBitClock = TxBitClock+Ts;
if TxBitClock > TSym
TxSymSent = TxSymSent+1; % get new bit
%
% We don't want the clock to increase to infinity so
% subtract off an integer number of Tb seconds.
%
TxBitClock = mod(TxBitClock,TSym);
%
% Get the new bit and appropriately.
%
TxOutput = TxBits(TxSymSent)*TxSigAmp;
end
[Rx,FilterState] = filter(b,a,TxOutput,FilterState);
%
% Add white Gaussian noise to the signal.
% First create unbiased (Monte Carlo) noise and then bias.
%
UnbiasedNoise = RxNoiseStd*(randn(1,1)+sqrt(-1)*randn(1,1));
BiasedNoise = CISBias*UnbiasedNoise;
%
% Calculate the CIS weight for this particular noise sample.
%
CISWeight = cgpdf(BiasedNoise,0,RxNoiseStd)./...
cgpdf(BiasedNoise,0,CISBias*RxNoiseStd);
%
% Since we are using white noise, the total CIS weight will be
% the product of the individuals CIS weights.
%
CISWeightIntegrator = CISWeightIntegrator*CISWeight;
Rx = (ChanGain*Rx)+BiasedNoise;
%
% Phase rotation due to receiver carrier synchronization error.
%
PhaseRotation = exp(sqrt(-1)*2*pi*(MeanCarrierPhaseError+...
(randn(1,1)*StdCarrierPhaseError))/360);
Rx = Rx*PhaseRotation;
probe1(probe1counter) = Rx; probe1counter = probe1counter+1;
%
% Update the Integrate and Dump Filter at the receiver.
%
RxIntegrator = RxIntegrator+Rx;
probe2(probe2counter) = RxIntegrator;...
probe2counter=probe2counter+1;
%
% Update the receiver clock, to see if it is time to
% sample and dump the integrator.
%
RxBitClock = RxBitClock+Ts;
RxTSym = TSym*(1+MeanSymbolSyncError+...
(StdSymbolSyncError*randn(1,1)));
if RxBitClock > RxTSym
RxSymDemod = RxSymDemod+1;
RxBitsI(RxSymDemod) = round(sign(real(RxIntegrator))+1)/2;
RxBitsQ(RxSymDemod) = round(sign(imag(RxIntegrator))+1)/2;
RxBitsCISWeight(RxSymDemod) = ...
CISWeightIntegrator*CISWeightIntegratorOld;
%
% Reset clock and dump the integrator.
%
RxBitClock = RxBitClock-TSym; RxIntegrator = 0;
CISWeightIntegratorOld = CISWeightIntegrator;
CISWeightIntegrator = 1;
end
end
%
% Implement differential decoder.
%
SinkBitsI = SourceBitsI*0;
SinkBitsQ = SourceBitsQ*0;
for k=2:RxSymDemod
SinkBitsI(k) = or(and(not(xor(RxBitsI(k),RxBitsQ(k))),...
xor(RxBitsI(k),RxBitsI(k-1))),...
and(xor(RxBitsI(k),RxBitsQ(k)),...
xor(RxBitsQ(k),RxBitsQ(k-1))));
SinkBitsQ(k) = or(and(not(xor(RxBitsI(k),RxBitsQ(k))),...
xor(RxBitsQ(k),RxBitsQ(k-1))),...
and(xor(RxBitsI(k),RxBitsQ(k)),...
xor(RxBitsI(k),RxBitsI(k-1))));
end;
%
% Look for best time delay between input and output, 100 bits.
%
[C,Lags] = vxcorr(SourceBitsI(10:110),SinkBitsI(10:110));
[MaxC,LocMaxC] = max(C);
BestLag = Lags(LocMaxC);
%
% Adjust time delay to match best lag.
%
if BestLag > 0
SourceBitsI = SourceBitsI(BestLag+1:length(SourceBitsI));
SourceBitsQ = SourceBitsQ(BestLag+1:length(SourceBitsQ));
RxBitsCISWeight = ...
RxBitsCISWeight(BestLag+1:length(RxBitsCISWeight));
elseif BestLag < 0
SinkBitsI = SinkBitsI(-BestLag+1:length(SinkBitsI));
SinkBitsQ = SinkBitsQ(-BestLag+1:length(SinkBitsQ));
RxBitsCISWeight = ...
RxBitsCISWeight(-BestLag+1:length(RxBitsCISWeight));
end
%
% Make all arrays the same length.
%
TotalBits = min(length(SourceBitsI),length(SinkBitsI));
TotalBits = TotalBits-20;
SourceBitsI = SourceBitsI(10:TotalBits);
SourceBitsQ = SourceBitsQ(10:TotalBits);
SinkBitsI = SinkBitsI(10:TotalBits);
SinkBitsQ = SinkBitsQ(10:TotalBits);
RxBitsCISWeight = RxBitsCISWeight(10:TotalBits);
%
% Find the number error events and the BER.
%
IErrors = SourceBitsI ^= SinkBitsI;
QErrors = SourceBitsQ ^= SinkBitsQ;
BER_CIS = ...
sum(IErrors.*RxBitsCISWeight)+sum(QErrors.*RxBitsCISWeight);
BER_CIS = BER_CIS/(2*length(SourceBitsI));
% End of function file.The program random_binary.m is given in Chapter 10, Appendix A. The program vxcorr.m is given in Chapter 10, Appendix B.