This chapter extends the basic material on Monte Carlo (MC) techniques explored in the preceding chapter. In Section 10.1 we consider two example simulations of communications systems. The first of these systems, a phase-shift key (PSK) digital communications system, although very simple, serves as a building block for simulations developed later in this book. This is followed by a more complicated simulation of a differential QPSK system in which the effects of phase and symbol synchronization errors are considered. In Section 10.2 we turn our attention to the semianalytic (SA) technique, which combines MC simulation and analysis.
The two methodologies explored in this chapter are quite different. Monte Carlo simulations require very little mathematical analysis and can be applied to any communication system for which the signal-processing algorithm required to represent each functional block in the block diagram of the system is known. Monte Carlo simulation is therefore a very general tool, but is applied at the expense of very long simulation run times, since, as we saw in the preceding chapter, a basic tradeoff exists between simulation accuracy and the time required to execute the simulation. Semianalytic simulation requires a higher level of analysis, but the payoff is a significantly reduced run time. In addition, execution of an MC simulation yields an estimate of the bit error rate (BER) at a single value of Eb/N0, while an SA simulation provides a complete curve of BER as a function of Eb/N0. We will see, however, that SA simulation is not a methodology that can be universally applied, since it is applicable to a restricted class of systems. For most applications, an SA simulation consumes a trivial amount of computer time and, therefore, is the preferred methodology when it can be applied.
As we saw in the previous chapter, the Monte Carlo technique, applied to the estimation of the BER of a digital communication system, is implemented by passing N data symbols through a simulation model of the system and counting the number of errors that occur. Assuming that passing N symbols through the simulation model results in Ne errors, the estimate of the BER is
We learned in the previous chapter that 
is a random variable, and accurate estimation of the BER requires that the estimator be unbiased and have small variance. Small variance requires that N be large and this in turn results in long computer run times. In the work to follow, the Monte Carlo technique is illustrated by giving two simple examples. Other examples are contained in the remainder of this book.
Example 10.1. (PSK).
For our first example consider the basic system illustrated in Figure 10.1. We assume binary PSK modulation with both signal points in the signal constellation lying in the direct (in-phase) channel. (Recall Example 9.3.) With this assumption, we can eliminate the quadrature channel from the simulation. The filter at the output of the modulator, which is assumed to be a third-order Butterworth filter with a bandwidth equal to the bit rate (BW = rb), leads to intersymbol interference (ISI). The purpose of the simulation is to determine the increase in BER resulting from the filter-induced ISI. The program for simulating the system is given in Appendix A. A block-serial approach is used in which blocks of 1,000 symbols are processed iteratively until N total symbols are processed. This was primarily done so that the MATLAB routine filter, which is a built-in MATLAB function implementing a time-domain convolution, could be used. As a built-in function it is very efficient and results in a significant reduction in the simulation run time. Note that one must ensure that the filter output is continuous from block to block. This is accomplished by using the initial condition parameter provided in filter.
The first problem is to determine the value of delay. There are a number of ways in which this can be accomplished. The most elegant way is to crosscorrelate the modulator output and the receiver output, as was done in Chapter 8 in the signal-to-noise ratio estimator. This will be the method used when we consider SA simulation. In order to illustrate the importance of correctly choosing the value of delay, we will use a different technique in this example. Specifically, we will choose a value of Eb/N0, simulate the system using different values of delay, and observe the results. The MATLAB routine for accomplishing this follows:
% File: c10_MCBPSKdelay.m
EbNodB = 6; % Eb/No (dB) value
z = 10.^(EbNodB/10); % convert to linear scale
delay = 0:8; % delay vector
BER = zeros(1,length(delay)); % initialize BER vector
Errors = zeros(1,length(delay)); % initialize Errors vector
BER_T = q(sqrt(2*z))*ones(1,length(delay));% theoretical BER vector
N = round(100./BER_T); % 100 errors for ideal (zero ISI) system
FilterSwitch = 1; % set filter switch (in=1 or out=0)
for k=1:length(delay)
[BER(k),Errors(k)] = c10_MCBPSKrun(N(k),z,delay(k),FilterSwitch)
end
semilogy(delay,BER,'o',delay,BER_T,'-'), grid;
xlabel('Delay'), ylabel('Bit Error Rate'),
% End of script file.Note that the preceding MATLAB script is essentially a combined preprocessor and postprocessor. (The simulation engine is the MATLAB function given in Appendix A.) The assumed value of Eb/N0 is 6 dB and delay is iterated from 0 to 8 samples. Since the sampling frequency is 10 samples per symbol, the step size of delay is 0.1Ts, where Ts is the symbol duration. The value of N is chosen so that a sufficient number of errors occur to ensure that the estimator variance is suitably small. In this case, we set N to 100/PT, where PT is the theoretical error probability for the additive, white, Gaussian noise (AWGN) case. The presence of ISI and other disturbances will of course increase the number of errors that occur for a given Eb/N0 over the average value of 100. Note that for each value of delay, both the value of the BER and the number of errors used to compute the BER are displayed. This allows the assumption of “a sufficient number of errors to form a reliable BER estimate” to be verified.
Executing the simulation yields the result illustrated in Figure 10.2. The various simulation results are indicated by the small circles, and the performance of the ideal (zero ISI) system operating in an AWGN environment with Eb/N0 = 6 dB is given by the solid line for reference. An incorrectly chosen value of delay clearly results in a value of BER that is too large. Since the value of the BER is a minimum at a delay of 5 samples, one might assume that 5 samples is the appropriate value of delay. However, since the delay is quantized to an integer number of sample periods, the value of 5 may not be precisely correct. Observation of Figure 10.2 implies that the corect delay value is most likely between 5 and 6 sampling periods. A more precise estimate of the correct value of delay can be determined by executing the simulation again with a higher sampling frequency (smaller sampling periods). (See Problem 10.3.) It should also be remembered that the estimator defined in (10.1) is a random variable and, as a result, any given value of BER may be too high or too low. The transmitter filter causes ISI, and the effect of ISI will prevent the BER from achieving the zero-ISI limit for any value of delay. If the transmitter filter is removed (FilterSwitch=0), the zero-ISI limit can be achieved.
Now that we know the appropriate value of delay, we can execute the simulation and determine the value of as a function of Eb/N0. The MATLAB script follows:
% File: c10_MCBPSKber.m
EbNodB = 0:8; % vector of Eb/No (dB) values
z = 10.^(EbNodB/10); % convert to linear scale
delay = 5; % enter delay value (samples)
BER = zeros(1,length(z)); % initialize BER vector
Errors = zeros(1,length(z)); % initialize Errors vector
BER_T = q(sqrt(2*z)); % theoretical (AWGN) BER vector
N = round(20./BER_T); % 20 errors for ideal (zero ISI) system
FilterSwitch = 1; % Tx filter out=0 or in=1
for k=1:length(z)
N(k) = max(1000,N(k)); % ensure at least one block processed
[BER(k),Errors(k)] = c10_MCBPSKrun(N(k),z(k),delay,FilterSwitch)
end
semilogy(EbNodB,BER,'o',EbNodB,BER_T)
xlabel('E_b/N_o - dB'), ylabel('Bit Error Rate'), grid
legend('System Under Study','AWGN Reference',0)
% End of script file.Note that Eb/N0 is stepped in 1 dB steps from 0 dB to 8 dB.
When executing a Monte Carlo simulation over a range of Eb/N0 values, using the same value of N for each value of Eb/N0 will result in an estimated value of BER that is based on a decreasing number of observed errors as Eb/N0 increases. As a result, the BER estimate at large values of Eb/N0 will be less reliable than the BER estimate at smaller values of Eb/N0. This problem can partially be overcome by setting N, the number of samples processed, to K/PT, where PT is the error probability for the AWGN case. Since system impairments such as ISI and synchronization errors will result in simulated values of that exceed PT, one will typically observe more than K errors in a given simulation run. The preceding MATLAB code uses K = 20. If N is determined in this fashion, values of N less than 1,000 are possible for sufficiently small values of Eb/N0. Since the simulation is based on the processing of sequential blocks of samples, with a block size of 1,000 symbols (10,000 samples), we must ensure that N > 1,000 so that at least one complete block is processed by the simulation. If N < 1,000, incorrect results will occur.
Executing the simulation provides the results illustrated in Figure 10.3. The simulation results are given by the small circles, and the BER of the ideal (zero ISI) results are given by the solid line. The increased BER resulting from the ISI caused by the filter is evident.
The block-serial technique used in this example will be used again in Chapter 18 when we examine a simulation of a simplified CDMA system. This simulation, although simple, will serve as a building block for simulations presented later in this book.
Example 10.2. (QPSK).
The previous example on BPSK modulation made a number of simplifying assumptions to keep the simulation code compact and the analysis tractable. This example of a QPSK system models some new sources of error, and includes some new parameters that may make it easier to relate the simulation results to those obtained from a physical communications system. Note for example that the channel attenuation, accounting for the propagation loss between transmitter and receiver, is included as a simulation parameter. Real (nonscaled) values for the symbol rate and the sampling frequency are also included in this simulation. The block diagram for this system is shown in Figure 10.4, and the code is given in Appendix B. Note that the transmitter filter, although illustrated in Figure 10.4, is not used in the simulations presented here, in order to reduce simulation execution time. Provision for including the transmitter filter is included in the simulation code given in Appendix B.
In coherent radio frequency (RF) systems, the receiver must provide carrier and symbol synchronization capabilities. Noise and distortion in the channel will make it impossible for the carrier and symbol synchronizers to operate perfectly. Incorrect carrier synchronization will result in a phase error, or phase rotation, of the received signal relative to the transmitted signal. The simulation provided in this example allows the user to simulate this phase error as a stochastic process. Symbol synchronization errors will result in the integrate-and-dump detector processing the received signal over the incorrect time interval. The simulation in this example also allows the user to examine the errors resulting from this effect.
As we know from basic communication theory, QPSK systems suffer from a problem called phase ambiguity. Since the channel introduces an unknown time delay to the signal, it will be impossible for the receiver to determine the absolute phase of the transmitted signal. For example, in a QPSK system, a transmitter may send the phase sequence 45°, 135°, 45°, and -45°. Suppose the channel introduces a time delay equal to 100.75 cycles of the RF carrier. The receiver will now mistakenly detect the original 45 degree signal to be -45 degrees, and make similar errors on the remaining symbols to produce the received sequence of -45°, -135°, 135°. If the information bits are contained in the absolute phase of the transmitted signal, the receiver will make a large number of errors. The solution to this problem involves encoding the information not in the absolute phase, but in the phase difference between symbols. For example, if the transmitter phase increases 90°, from 45° to 135° between the first and second symbols, the receiver will detect these two signals as -135° and 45°, which still shows a phase increase of 90°. Differential encoding is implemented in the MATLAB code given in Appendix B, which is the main simulation code for the QPSK system. All simulations presented in this example use repeated calls to this code.
As in the previous example, the time delay through the system must be determined. The following MATLAB program determines the optimal time delay to be used at the receiver to account for the signal propagation delay through the system:
% File: c10_MCQPSKdelay.m Eb = 23; No = -50; % Eb (dBm) and No (dBm/Hz) ChannelAttenuation = 70; % channel attenuation in dB N = 1000; delay = -0.1:0.1:0.5; EbNo = 10.^(((Eb-ChannelAttenuation)-No)/10); BER_MC = zeros(size(delay));
for k=1:length(delay)
BER_MC(k) = c10_MCQPSKrun(N,Eb,No,ChannelAttenuation,...
delay(k),0,0,0);
disp(['Simulation ',num2str(k*100/length(delay)),'% Complete']);
end
BER_T = 0.5*erfc(sqrt(EbNo))*ones(size(delay)); % Theoretical BER
semilogy(delay,BER_MC,'o',delay,2*BER_T,'-') % Plot BER vs Delay
xlabel('Delay (symbols)'), ylabel('Bit Error Rate'),
legend('MC BER Estimate','Theoretical BER'), grid;
% End of script file.Since no channel filter was used, the optimal delay is zero symbols, as shown in Figure 10.5. Note that we have measured delay with respect to the symbol period rather than in samples, as was done in the previous example.
Now that we know the delay, we will measure the sensitivity of the BER to static synchronization phase error. The phase error is measured from 0 to 90 degrees in 10-degree increments. The code for accomplishing this follows:
% File: c10_MCQPSKphasesync.m
PhaseError = 0:10:90; % Phase Error at Receiver
Eb = 24; No = -50; % Eb (dBm) and No (dBm/Hz)
ChannelAttenuation = 70; % dB
EbNo = 10.^((Eb-ChannelAttenuation-No)/10);
BER_T = 0.5*erfc(sqrt(EbNo)*ones(size(PhaseError)));
N = round(100./BER_T);
BER_MC = zeros(size(PhaseError));
for k=1:length(PhaseError)
BER_MC(k) = c10_MCQPSKrun(N(k),Eb,No,ChannelAttenuation,0,0,...
PhaseError(k),0);
disp(['Simulation ',num2str(k*100/length(PhaseError)),'...
% Complete']);
end
semilogy(PhaseError,BER_MC,'o',PhaseError,2*BER_T,'-')
xlabel('Phase Error (Degrees)'),
ylabel('Bit Error Rate'),
legend('MC BER Estimate','Theoretical BER'), grid;
% End of script file.Executing the simulation yields the result illustrated in Figure 10.6. Figure 10.6 shows that the BER, as determined by the simulation, reaches a maximum at a phase error of 45 degrees and decreases back to the optimal value (value for zero synchronization phase error) for a phase error of 0 or 90 degrees. This behavior is due to the differential encoder.
Now that we know the optimal phase shift and time delay for the channel, we can measure the BER as a function of the signal-to-noise ratio (SNR). The MATLAB code for accomplishing this follows:
% File: c10_MCQPSKber.m
Eb = 22:0.5:26; No = -50; % Eb (dBm) and No (dBm/Hz)
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 = round(100./BER_T); % Symbols to transmit
BER_MC = zeros(size(Eb)); % Initialize BER vector
for k=1:length(Eb) % Main Loop
BER_MC(k) = c10_MCQPSKrun(N(k),Eb(k),No,ChannelAttenuation,...
0,0,0,0);
disp(['Simulation ',num2str(k*100/length(Eb)),'% Complete']);
end
semilogy(EbNodB,BER_MC,'o',EbNodB,2*BER_T,'-')
xlabel('Eb/No (dB)'), ylabel('Bit Error Rate'),
legend('MC BER Estimate','Theoretical BER'), grid;
% End of script file.Executing the simulation yields the plot shown in Figure 10.7. We see that the simulated result is very close to the theoretical AWGN result. This provides at least a partial sanity check on the simulation.
Next we examine the impact of phase jitter on the system BER. The phase error process is modeled as white Gaussian noise. The MATLAB code for the simulation follows:
% File: c10_MCQPSKPhaseJitter.m
PhaseBias = 0; PhaseJitter = 0:2:30;
Eb = 24; No = -50; % Eb (dBm) and No (dBm/Hz)
ChannelAttenuation = 70; % dB
EbNo = 10.^((Eb-ChannelAttenuation-No)/10);
BER_T = 0.5*erfc(sqrt(EbNo)*ones(size(PhaseJitter)));
N=round(100./BER_T);
BER_MC = zeros(size(PhaseJitter));
for k=1:length(PhaseJitter)
BER_MC(k) = c10_MCQPSKrun(N(k),Eb,No,ChannelAttenuation,0,0,...
PhaseBias,PhaseJitter(k));
disp(['Simulation ',num2str(k*100/length(PhaseJitter)),'...
% Complete']);
end
semilogy(PhaseJitter,BER_MC,'o',PhaseJitter,2*BER_T,'-')
xlabel('Phase Error Std. Dev. (Degrees)'),
ylabel('Bit Error Rate'),
legend('MC BER Estimate','Theoretical BER'), grid;
% End of script file.Executing the simulation yields the result illustrated in Figure 10.8. As expected, the BER increases as the standard deviation of the phase jitter increases. In many system simulations it is not appropriate to model phase jitter as a white-noise process. Should this be the case, a finite impulse response (FIR) filter can be designed to realize the required power spectral density (PSD) of the phase jitter process.
The final simulation in this sequence examines the sensitivity of BER to symbol timing error. The MATLAB code follows:
% File: c10_MCQPSKSymJitter.m
SymJitter = 0:0.02:0.2;
Eb = 24; No = -50; % Eb (dBm) and No (dBm/Hz)
ChannelAttenuation = 70; % channel attenuation in dB
EbNo = 10.^((Eb-ChannelAttenuation-No)/10);
BER_T = 0.5*erfc(sqrt(EbNo)*ones(size(SymJitter)));
N=round(100./BER_T);
BER_MC = zeros(size(SymJitter));
for k=1:length(SymJitter)
BER_MC(k) = c10_MCQPSKrun(N(k),Eb,...
No,ChannelAttenuation,0,SymJitter(k),0,0);
disp(['Simulation ',num2str(k*100/length(SymJitter)),...
'% Complete']);
end
semilogy(SymJitter,BER_MC,'o',SymJitter,2*BER_T,'-')
xlabel('Symbol Timing Error Std. Dev. (Symbols)'),
ylabel('Bit Error Rate'),
legend('MC BER Estimate','Theoretical BER'), grid;
% End of script file.The result of executing the simulation is shown in Figure 10.9. Just as in the phase jitter case, the symbol synchronization error is modeled as a white Gaussian stochastic process. Once again, if memory effects in the symbol jitter process must be accurately modeled, an FIR filter can be designed to realize the required PSD.
In addition, in this simulation the transmitted symbols are crosscorrelated (note the use of the function vxcorr) in order to ensure that the transmitted and received symbols are properly aligned so that the BER is correctly determined. In future simulations, including Examples 10.3 and 10.4 to follow, the crosscorrelation technique will be used to calculate the appropriate value of delay and a separate simulation to determine this parameter will not be required. A separate simulation program was used here to illustrate the sensitivity of the simulation results to this important parameter. This simulation will be encountered again in Chapter 16 when we consider importance sampling.
As we have seen, the Monte Carlo simulation method is completely general and may be applied to any system for which simulation models of the various system building blocks can be defined, or at least approximated, in terms of a numerical (digital signal processing, DSP) algorithm. No analytical knowledge, outside of that required to implement the subsystem models, is required. The price paid for using Monte Carlo methods is the run time required for executing the simulation. If the system and the channel model are complicated, and the BER is low, the required run time is sometimes so long that the use of Monte Carlo techniques becomes impractical for all but the most important simulations.
In the work to follow we stress the estimation of the BER, since the BER is the most common measure used to evaluate the performance of digital communication systems. In executing a simulation to estimate the BER, information is collected that allows other items of interest to be determined. These include waveforms, eye diagrams, signal constellations, and PSD estimates at various points in the system of interest.
Fortunately there are alternatives to the pure Monte Carlo method. One of the most powerful of these alternatives is the semianalytic method, in which analytical and simulation techniques are used together in a way that yields very rapid estimation of the BER. The SA simulation method, like all rapid simulation techniques, allows analytical knowledge to be traded off against simulation run time.
The block diagram of a simple system, to which the SA simulation technique is applicable, is illustrated in Figure 10.10. The kth transmitted symbol is denoted dk, and the corresponding symbol at the receiver output is 
. The transmitted symbol is received correctly if 
, and an error is made if 
. As we know from basic communication theory, the quantity Vk is the decision statistic for the kth transmitted symbol, and the receiver makes a decision by comparing the value of Vk with a threshold T.
The decision statistic Vk, shown in Figure 10.10, is a function of three components:
where Sk is the component of Vk due to the transmitted signal; Dk results from system-induced distortion such as ISI due to filtering or multipath, and Nk is due to the channel disturbances such as noise and interference. In applying semianalytic simulation, the combined effects of Sk and Dk are determined by an MC simulation, and the effects of noise, represented by Nk, are treated analytically. The SA simulation technique is applicable to any system in which the probability density function of the noise component of Vk can be analytically determined. A simple case is for the additive Gaussian-noise channel in which the system is linear from the point where noise is injected to the point at which the decision statistic Vk is defined. This follows since any linear transformation of a Gaussian random process yields a Gaussian process. Referring to Figure 10.10, we see that if the channel noise is Gaussian, the decision statistic Vk is a Gaussian random variable with the mean determined by Sk and Dk. Thus, semianalytic simulation is a combination of MC simulation and analysis.
As an example, consider a binary PSK (BPSK) system operating an AWGN (additive, white, Gaussian noise) environment. For the moment we neglect the transmitter filter and assume full-response signaling.[1] The pdf of Vk, conditioned on transmission of a binary 1 or binary 0, is Gaussian, as illustrated in Figure 10.11. The probability density functions (pdfs) of the decision statistic, Vk, conditioned on dk = 0 and dk = 1, are given by
and
where v1 and v2 are the means of the random variable Vk conditioned on dk = 0 and dk = 1, respectively. The probability of error conditioned on dk = 0 is
where T is the decision threshold, and the probability of error conditioned on dk = 1 is
If the binary symbols dk = 0 and dk = 1 are transmitted with equal probability, the optimum threshold T is the point at which the two conditional pdfs are equal. For this case the two conditional error probabilities are equal and the overall error probability is
which is
The probability for an AWGN environment is usually expressed in terms of the Gaussian Q function. This gives
where
defines the Gaussian Q function. Note that with equal energy signals the threshold T is zero and therefore v2 = −v1. Thus, determination of v1 and σn completely determine system BER. In order to determine the error probability we need only to develop a simulation for estimating v1 and σn. The Monte Carlo technique of counting errors is not required.
The value of v1 can be determined by executing a noiseless simulation. If the channel noise is removed, the two conditional pfds shown in Figure 10.11(a) collapse to impulse functions (σn = 0) as shown in Figure 10.11(b). Each of these impulse functions has unity area, and the locations of the impulse functions define v1 and v2.
The value of σn is determined by executing a simple simulation of that portion of the system through which the noise passes. For the system being considered this consists of the receiver, which is modeled as an integrate-and-dump symbol detector. Assume that this portion of the system has transfer function H(f). If white (delta-correlated) noise having a two-sided power spectral density of N0/2 is input to the matched-filter receiver, the variance of the random variable Vk is
As we saw in Chapter 7, the equivalent noise bandwidth is defined as
and is the equivalent noise bandwidth of the receiver. It follows that the probability of error is given by
Note that even though the system is AWGN, the nonlinear amplifier may affect system performance, since the nonlinear amplifier will affect the shape of the transmitted signals and this, in turn, will affect the value of v1.
We now consider the presence of the transmitter filter. The effect of this filter is to spread, in time, the energy associated with the transmitted symbols beyond the symbol period giving rise to intersymbol interference. If the memory length of this filter is two symbols, the error probability associated with the transmission symbol will depend not only on the transmitted symbol but also on the previously transmitted symbol. As a result, the calculation of the probability of error will involve four conditional probability density functions rather than two conditional probability density functions as shown in Figure 10.11. This is illustrated in Figure 10.12(a).
As before, execution of a noiseless simulation will yield the values of v1, v2, v3, and v4. The system probability of error becomes
The extension to a memory length of M symbols is obvious.
In applying the semianalytic technique we make use of the idea of an equivalent noise source. We have seen that the decision statistic Vk is a function of three components. In other words, Vk = f(Sk,Dk,Nk) where Sk is due to the signal; Dk results from system-induced distortion such as ISI, and Nk is due to noise. The effects of Sk and Dk are determined by an MC simulation, and the effects of noise, represented by Nk, are treated, as we have seen, analytically. If a noise-free simulation is executed, the resulting sufficient statistic, which is denoted Vk,nf, will be a function of only Sk and Dk. To this statistic is added a random variable Nk having the variance defined by (10.11). Thus
The random variable Nk may be viewed as a sample from an equivalent noise source ne(t) as shown in Figure 10.13. This equivalent noise source contains the combined effects of themal noise, interference, and other channel impairments reflected to the integrator output of the integrate-and-dump detector. If the channel noise is white, the impulse response, or equivalently the transfer function, defined in (10.11) is used to transform the channel noise to the integrator output.
We now briefly consider the development of an algorithm for the determination of the BER in a binary PSK system using semianalytic simulation. We do this in a way that is easily extendable to QPSK. Consider the signal constellation illustrated in Figure 10.14. The transmitted signal points are denoted S1 and S2 and the corresponding decision regions are denoted D1 and D2. A correct decision is made at the receiver if Si is transmitted and the received signal falls in region Di; otherwise an error occurs. In Figure 10.14 we assume that S1 is transmitted and 
is received. As discussed in the previous section, S1 and differ because of intersymbol interference, nonlinear distortion, of other signal-degrading effects. The difference between S1 and is denoted dx. The conditional error probability, conditioned on the transmission of S1 is
which is
In terms of the Gaussian Q-function, the preceding equation becomes
Thus, knowledge of , determined using MC simulation, and σn, allows the conditional BER to be determined. In determining σn the value of BN is found from the simulated impulse response h[n].
Assume that Sk is the kth transmitted bit in a simulated sequence of N bits. For each value of k, 1 ≤ k ≤ N, Sk will be S1 or S2. The conditional BER is
The overall BER, obtained by averaging over the entire sequence of N bits, is given by
Example 10.3. (PSK).
The MATLAB code for executing a semianalytic simulation of a PSK system is given in Appendix C. The methodology used is that presented in the preceding paragraphs. Due to symmetry, the received symbols are rotated to positive values. The bandwidth of the transmitter filter, which gives rise to ISI, is equal to the bit rate. The increase in the BER resulting from ISI is clearly seen in Figure 10.15.
We now consider a semianalytic estimator for the symbol error probability PS in a QPSK system.[2] Since a QPSK signal constellation has four signal points rather than two, and since the signal space has two dimensions rather than one, the semianalytic estimator for QPSK is different from the estimator for PSK in that a dimension must be added for the quadrature channel.
Consider the signal constellation illustrated in Figure 10.16. The transmitted signal points are denoted Si, i = 1, 2, 3, 4, and the decision regions are denoted Di, i = 1, 2, 3, 4. As in the preceding section, a correct decision is made at the receiver if Si is transmitted and the received signal falls in decision region Di; otherwise an error occurs. In Figure 10.16 it is assumed that S1 is transmitted, and the noiseless received signal is denoted . As a result of intersymbol interference and distortion, 
. It is rather than S1 that is determined by the semianalytic simulation, since the simulation will account for the effects of intersymbol interference but not the effects of noise. The direct and quadrature components of are denoted 
and 
, respectively, where 
and 
. When noise is considered, by adding nx and ny to and , respectively, a correct decision is made, conditioned on S1 transmitted, if 
. An error is made if 
. Keep in mind that since we are developing a semianalytic estimator, the impact of noise is treated analytically and does not appear in Figure 10.16.
The problem is to determine the noise components nx and ny that will result in an error given the received (noiseless) point in signal space . The problem is very similar to the PSK example just considered. The essential difference is that we are working in two dimensions rather than one. We assume that the direct and quadrature additive noise components are uncorrelated and jointly Gaussian. Thus, given that S1 is transmitted and is received, an error is made if
where nx and ny are the direct and quadrature noise components, and 
represents the noise variance. In order to simplify the notation let
and
With these changes (10.21) becomes
This can be bounded by the expression
where the bound occurs since the decision region D3 appears twice in (10.25). From the definition of the decision regions we can write
Recognizing that two of the four integrals in (10.26) are equal to one yields
Substituting (10.22) and (10.23) in the preceding expression, and using the definitions of and , yields the bound on the conditional error probability. This conditional error probability bound is
where, as always, Q(·) is the Gaussian Q function. By symmetry, the conditional error probability is the same for any of the four possible transmitted symbols.
As with PSK assume that Sk is the kth transmitted symbol in a simulated sequence of N symbols. For each value of k, 1 ≤ k ≤ N, Sk will be S1, S2, S3, or S4. The bound on the conditional symbol error rate is, from (10.28):
The overall symbol error rate, obtained by averaging the conditional symbol error probability over the entire sequence of N symbols, is given by
The bit error rate, PE, is PS/2. Note that in the PSK case we obtained an exact solution, whereas in the case of QPSK we have a bound. The technique used here to develop a semianalytic estimator is easily extended to MPSK and QAM [1].
The estimator developed here will be used throughout the remainder of this book for evaluating the performance of a number of systems. Included will be examples illustrating the effect of multipath and fading in a wireless system and the effect of nonlinear distortion in a frequency multiplexed satellite communications system.
Example 10.4. (QPSK).
The MATLAB code for executing a semianalytic simulation of a QPSK system is given in Appendix D. The simulation is run to examine the effect of the ISI resulting from transmitter filtering. The filter bandwidth is set equal to the symbol rate (one-half the bit rate, i.e., BW = rb/2). Since the signal constellation is symmetric, all received signal points are rotated to the first quadrant as discussed in the preceding paragraphs.
Executing the simulation yields the signal constellation and BER illustrated in Figure 10.17. The received signal constellation is shown in the left-hand pane of Figure 10.17. Note that the received signal constellation no longer consists of 4 points, as is the case for ideal QPSK, but now consists of 16 points. In order to understand the reason for this, assume that the signal point in the first quadrant represents the data bit 00 and that the system memory, as a result of ISI, is two symbols (the current and the previous transmitted symbols). As a result, four signal points will result from the transmission of 00. These four signal points correspond to 00|00, 00|01, 00|10, and 00|11, where the vertical bar delineates the current symbol (00) and the previously transmitted symbols. Note also that each of the four points in the first quadrant are composed of points that are slightly scattered. This scattering results from the fact that the system exhibits a memory length that exceeds two symbols, although the effect of this additional memory is small. The left-hand pane of Figure 10.17 illustrates the BER of the system with transmitter filtering. The AWGN result is also illustrated for reference. The increase in the BER resulting from the ISI is clearly seen.
In applying semianalytic techniques to a system having memory, it is important to use a data source that generates sequences exhibiting all possible combinations of data symbols for the given memory length of the system. For example, if the memory length is three (current symbol plus the two preceding symbols) the symbol error probability is given by
The error probability will, in general, be different for each (Sk,Sk−1,Sk−2) sequence. Thus, all combinations of Sk, Sk−1, and Sk−2, must appear an equal number of times to properly account for the memory efficts. In general, if a binary system exhibits significant memory spanning N symbols, then all binary sequences of length N should be generated by the data source an equal number of times in the simulation. For a binary system there are 2N sequences of length N. There are three popular ways to accomplish, or at leat approximate, this requirement as follows:
If N is reasonably large one may wish to use a PN Sequence for the data source. As discussed in Chapter 7, the number of sequences generated will not be 2N as desired but will be L = 2N − 1, since the sequence of N consecutive zeros will not occur. The result will be an unbalanced sequence having
ones and 
zeros. If N is large, this effect is negligible. Note that one is free to select N greater than the memory length in order to mitigate this effect. Choosing N larger than necessary will, however, result in a longer simulation execution time.If a perfectly balanced sequence is desired, a deBruijn sequence [2] may be used. As briefly discussed in Chapter 7, a deBruijn sequence is formed by adding an extra zero at the point where there are N − 1 zeros in the output of a PN sequence generator.
One may of course simply execute the semianalytic sequence using random data. If this sequence is sufficiently long, all data symbol combinations will occur approximately the same number of times. This is the approach taken in Examples 10.3 and 10.4.
In this chapter, simulation examples of binary PSK and differential QPSK communication systems were presented. Strict Monte Carlo simulations were first developed. These simulations were easily developed using the concepts presented in the previous chapter. The PSK system was very simple and served to illustrate the basic concepts. The only degrading effects were intersymbol interference and additive channel noise. The differential QPSK example considered the simulation of a much more realistic system.
We next considered semianalytic simulation. Both PSK and QPSK illustrated that the semianalytic estimators for BER are different for the PSK and QPSK cases. Thus, a unique procedure for conducting a semianalytic simulation does not exist. While the estimators are quite different, the methodologies are the same, in that the semianalytic simulation captures all deterministic system perturbations, such as intersymbol interference and distortion due to nonlinearities through a conventional Monte Carlo simulation. The effects of noise and other stochastic effects are dealt with analytically. This requires that the pdf of the samples upon which a bit or symbol decision is made is known. The simplest case, and the case most often used, assumes that the noise is Gaussian and that the system is linear from the point at which the noise enters the system to the point where bit or symbol decisions are made is linear. In this case, the pdf of the decision statistic is Gaussian, and the Monte Carlo simulation is executed to establish the mean of the decision statistic. We saw that for those cases in which the semianalytic method can be used, very fast simulations result.
10.1 | Verify that the MATLAB code segment used in Appendix A [Btr,Atr]=butter(5,0.2) with a sampling frequency of 10 samples/symbol yields a filter bandwidth equal to the symbol rate. |
10.2 | Rerun the simulation described in Example 10.1 using filter bandwidths of one-half the symbol rate and double the symbol rate. Compare the result of these two simulations with the result given in Example 10.1. |
10.3 | Rerun the simulation described in Example 10.1 using a sampling rate of 20 samples per symbol. Does this result in an improved estimate of delay? Explain. Estimate the BER using 20 samples/symbol and compare the result with that obtained in Example 10.1 using 10 samples/symbol. |
10.4 | Using the appropriate MATLAB routines, compare the required execution times for the Monte Carlo and semianalytic simulations of the binary PSK system given in Examples 10.1 and 10.3, respectively. |
10.5 | The semianalytic BER estimator for the PSK system (Appendix C) contains the line of code nbwideal=1/(2*tb) Explain the purpose of this line of code and verify that it is correct. The equivalent line of code for the BER estimator for the QPSK system (Appendix D) is nbwideal=1/(2*tb*2) Verify the correctness of this line of code. |
10.6 | Modify the simulation given in Example 10.1 so that the PSK system is simulated using a symbol-by-symbol approach rather than by using a block serial approach. In other words, the random binary data source will generate binary bits (0 or 1), and the waveform samples corresponding to these binary symbols will be repeated the required number of times to satisfy a given samples/symbol specification. (You may wish to review Chapter 5 in order to derive an efficient filter simulation using MATLAB to establish the necessary numerator and denominator polynomials for the filter transfer function and use them in a sample-by-sample simulation.) |
10.7 | Example 10.2 examines the Monte Carlo simulation of a differential QPSK system. Rewrite this simulation for QPSK rather than for differential QPSK. Sanity check the simulation result by comparing it with the theoretical result for QPSK. |
Appendix A: Simulation Code for Example 10.1
% File: c10_MCBPSKrun.m
function [BER,Errors]=MCBPSKrun(N,EbNo,delay,FilterSwitch)
SamplesPerSymbol = 10; % samples per symbol
BlockSize = 1000; % block size
NoiseSigma = sqrt(SamplesPerSymbol/(2*EbNo)); % scale noise level
DetectedSymbols = zeros(1,BlockSize); % initialize vector
NumberOfBlocks = floor(N/BlockSize); % number of blocks
% processed
[BTx,ATx] = butter(5,2/SamplesPerSymbol); % compute filter
% parameters
[TxOutput,TxFilterState] = filter(BTx,ATx,0); % initialize state
% vector
BRx = ones(1,SamplesPerSymbol); ARx=1; % matched filter
% parameters
Errors = 0; % initialize error
% counter
%
% Simulation loop begins here.
%
for Block=1:NumberOfBlocks
%
% Generate transmitted symbols
%
[SymbolSamples,TxSymbols] = random_binary(BlockSize,...
SamplesPerSymbol);
%
% Transmitter filter if desired.
%
if FilterSwitch==0
TxOutput = SymbolSamples;
else
[TxOutput,TxFilterState] = filter(BTx,ATx,SymbolSamples,...
TxFilterState);
end
%
% Generate channel noise.
%
NoiseSamples = NoiseSigma*randn(size(TxOutput));
%
% Add signal and noise.
%
RxInput = TxOutput + NoiseSamples;
%
% Pass Received signal through matched filter.
%
IntegratorOutput = filter(BRx,ARx,RxInput);
%
% Sample matched filter output every SamplesPerSymbol samples,
% compare to transmitted bit, and count errors.
%
for k=1:BlockSize,
m = k*SamplesPerSymbol+delay;
if (m < length(IntegratorOutput))
DetectedSymbols(k)=(1-sign(IntegratorOutput(m)))/2;
if (DetectedSymbols(k) ~= TxSymbols(k))
Errors = Errors + 1;
end
end
end
end
BER = Errors/(BlockSize*NumberOfBlocks); % calculate BER
% End of function file.% file: random_binary.m
function [x, bits] = random_binary(nbits,nsamples)
% This function genrates a random binary waveform of length nbits
% sampled at a rate of nsamples/bit.
x = zeros(1,nbits*nsamples);
bits = round(rand(1,nbits));
for m=1:nbits
for n=1:nsamples
index = (m-1)*nsamples + n;
x(1,index) = (-1)^bits(m);
end
end
% End of function file.Appendix B: Simulation Code for Example 10.2
% file c10_MCQPSKrun.m
function BER_MC=MCQPSKrun(N,Eb,No,ChanAtt,...
TimingBias,TimingJitter,PhaseBias,PhaseJitter)
fs = 1e+6; % sampling Rate (samples/second)
SymRate = 1e+5; % symbol rate (symbols/second)
Ts = 1/fs; % sampling period
TSym = 1/SymRate; % symbol period
SymToSend = N; % symbols to be transmitted
ChanBW = 4.99e+5; % bandwidth of channel (Hz)
MeanCarrierPhaseError = PhaseBias; % mean of carrier phase
StdCarrierPhaseError = PhaseJitter; % stdev of phase error
MeanSymbolSyncError = TimingBias; % mean of symbol sync error
StdSymbolSyncError = TimingJitter; % stdev of symbol sync error
ChanGain = 10^(-ChanAtt/20); % channel gain (linear units)
TxBitClock = Ts/2; % transmitter bit clock
RxBitClock = Ts/2; % receiver bit clock
%
% Standard deviation of noise and signal amplitude at receiver input.
%
RxNoiseStd = sqrt((10^((No-30)/10))*(fs/2)); % stdev of noise
TxSigAmp = sqrt(10^((Eb-30)/10)*SymRate); % signal amplitude
%
% Allocate some memory for probes.
%
SampPerSym = fs/SymRate;
probe1 = zeros((SymToSend+1)*SampPerSym,1);
probe1counter = 1;
probe2 = zeros((SymToSend+1)*SampPerSym,1);
probe2counter = 1;
%
% Counters to keep track of how many symbols have have been sent.
%
TxSymSent = 1;
RxSymDemod = 0;
%
% 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; % set first I bit
TxBitsQ = SourceBitsQ*0; % set first Q bit
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));
%
RxIntegrator = 0; % initialize receiver integrator
TxBitClock = 2*TSym; % initialize transmitter
%
% Design the channel filter, and create the filter state array.
%
[b,a] = butter(2,ChanBW/(fs/2));
b = [1]; a = [1]; % filter bypassed
[junk,FilterState] = filter(b,a,0);
%
% Begin simulation loop.
%
while TxSymSent < SymToSend
%
% Update the transmitter's clock, and see
% if it is time to get new data bits
%
TxBitClock = TxBitClock+Ts;
if TxBitClock > TSym
%
% Time to get new bits
%
TxSymSent = TxSymSent+1;
%
% 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 scale it up appropriately.
%
TxOutput = TxBits(TxSymSent)*TxSigAmp;
end
%
% Pass the transmitted signal through the channel filter.
%
[Rx,FilterState] = filter(b,a,TxOutput,FilterState);
%
% Add white Gaussian noise to the signal.
%
Rx = (ChanGain*Rx)+(RxNoiseStd*(randn(1,1)+sqrt(-1)*randn(1,1)));
%
% 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;
xTSym = TSym*(1+MeanSymbolSyncError+...
(StdSymbolSyncError*randn(1,1)));
if RxBitClock > RxTSym % time to demodulate symbol
RxSymDemod = RxSymDemod+1;
RxBitsI(RxSymDemod) = round(sign(real(RxIntegrator))+1)/2;
RxBitsQ(RxSymDemod) = round(sign(imag(RxIntegrator))+1)/2;
RxBitClock = RxBitClock - TSym; % reset receive clock
RxIntegrator = 0; % reset integrator
end
end
%
% Differential decoder.
%
SinkBitsI = SourceBitsI*0; % set first I sink bit
SinkBitsQ = SourceBitsQ*0; % set first Q sink bit
%
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 for 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));
elseif BestLag < 0
SinkBitsI = SinkBitsI(-BestLag+1:length(SinkBitsI));
SinkBitsQ = SinkBitsQ(-BestLag+1:length(SinkBitsQ));
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);
%
% Find the number of errors and the BER.
%
Errors = sum(SourceBitsI ~= SinkBitsI) + sum(SourceBitsQ ~=...
SinkBitsQ);
BER_MC = Errors/(2*length(SourceBitsI));
% End of function file.Program random_binary.m is defined in Appendix A of this chapter.
% File: vxcorr.m
function [c,lags] = vxcorr(a,b)
% This function calculates the unscaled cross-correlation of 2
% vectors of the same length. The output length(c) is
% length(a)+length(b)-1. It is a simplified function of xcorr
% function in matlabR12 using the definition:
% c(m) = E[a(n+m)*conj(b(n))] = E[a(n)*conj(b(n-m))]
%
a = a(:); % convert a to column vector
b = b(:); % convert b to column vector
M = length(a); % same as length(b)
maxlag = M-1; % maximum value of lag
lags = [-maxlag:maxlag]'; % vector of lags
A = fft(a,2^nextpow2(2*M-1)); % fft of A
B = fft(b,2^nextpow2(2*M-1)); % fft of B
c = ifft(A.*conj(B)); % crosscorrelation
%
% Move negative lags before positive lags.
%
c = [c(end-maxlag+1:end,1);c(1:maxlag+1,1)];
%
% Return row vector if a, b are row vectors.
%
[nr nc]=size(a);
if(nr>nc)
c=c.';
lags=lags.';
end
% End of function file.Appendix C: Simulation Code for Example 10.3
% File: c10_PSKSA.m
NN = 256; % number of symbols
tb = 1; % bit file
p0 = 1; % power
fs = 16; % samples/symbol
ebn0db = [0:1:8]; % Eb/No vector in dB
[bt,at] = butter(5,2/fs); % transmitter filter parameters
x = random_binary(NN,fs); % establish PSK signal
y1 = x; % save signal
y2a = y1*sqrt(p0); % scale amplitude
y2 = filter(bt,at,y2a); % transmitter output
br = ones(1,fs); br = br/fs; ar = 1; % matched filter parameters
y = filter(br,ar,y2); % matched filter output
%
% End of simulation.
%
% The following code sets up the semianalytic estimator. Find the
% max. magnitude of the cross correlation and the corresponding lag.
%
[cor lags] = vxcorr(x,y); % compute crosscorrelation
[cmax nmax] = max(abs(cor)); % maximum of crosscorrelation
timelag = lags(nmax); % lag at max crosscorrelation
theta = angle(cor(nmax)); % determine angle
y = y*exp(-i*theta); % derotate
%
% Noise BW calibration.
%
hh = impz(br,ar); % receiver impulse response
nbw = (fs/2)*sum(hh.^2); % noise bandwidth
%
% Delay the input and do BER estimation on the NN-20+1 128 bits.
% Use middle sample. Make sure the index does not exceed number
% of input points. Eb should be computed at the receiver input.
%
index = (10*fs+8:fs:(NN-10)*fs+8);
xx = x(index);
yy = y(index-timelag+1);
eb = tb*sum(abs(y2).^2)/length(y2);
eb = eb/2;
[peideal,pesystem] = psk_berest(xx,yy,ebn0db,eb,tb,nbw);
semilogy(ebn0db,pesystem,'ro-',ebn0db,peideal); grid;
xlabel('E_b/N_0 (dB)'), ylabel('Bit Error Rate')
legend('System Under Study','AWGN Reference',0)
% End of script file.Program random_binary.m is defined in Appendix A of this chapter. Program vxcorr.m is defined in Appendix B of this chapter.
% File: psk_berest.m
function [peideal,pesystem] = psk_berest(xx,yy,ebn0db,eb,tb,nbw)
% ebn0db is an array of Eb/No values in db (specified at the
% receiver input); tb is the bit duration and nbw is the noise BW
% xx is the reference (ideal) input; yy is the filtered output;
%
nx = length(xx);
%
% For comparision purposes, set the noise BW of the ideal
% receiver (integrate and dump) to be equal to rs/2.
%
nbwideal = 1/(2*tb); % noise bandwidth
for m=1:length(ebn0db)
peideal(m) = 0.0; pesystem(m) = 0.0; % initialize
%
% Find n0 and the variance of the noise.
%
ebn0(m) = 10^(ebn0db(m)/10); % dB to linear
n0 = eb/ebn0(m); % noise power
sigma = sqrt(n0*nbw*2); % variance
sigma1 = sqrt(n0*nbwideal*2); % variance of ideal
%
% Multiply the input constellation/signal by a scale factor so
% that input constellation and the constellations/signal at the
% input to receive filter have the same ave power
% a = sqrt(2*eb/(2*tb)).
%
b = sqrt(2*eb/tb)/sqrt(sum(abs(xx).^2)/nx);
d1 = b*abs(xx);
d3 = abs(yy);
peideal(m) = sum(q(d1/sigma1));
pesystem(m) = sum(q(d3/sigma));
end
peideal = peideal/nx;
pesystem = pesystem/nx;
% End of function file.Appendix D: Simulation Code for Example 10.4
% File: c14_QPSKSA.m
%
% Default parameters
%
NN = 256; % number of symbols
tb = 0.5; % bit time
p0 = 1; % power
fs = 16; % samples/symbol
ebn0db = [0:1:10]; % Eb/N0 vector
[b,a] = butter(5,1/16); % transmitter filter parameters
%
% Establish QPSK signals
%
x = random_binary(NN,fs)+i*random_binary(NN,fs); % QPSK signal
y1 = x; % save signal
y2a = y1*sqrt(p0); % scale amplitude
%
% Transmitter filter
%
y2 = filter(b,a,y2a); % filtered signal
%
% Matched filter
%
b = ones(1,fs); b = b/fs; a = 1; % matched filter parameters
y = filter(b,a,y2); % matched filter output
%
% End of simulation
%
% Use the semianalytic BER estimator. The following sets
% up the semi analytic estimator. Find the maximum magnitude
% of the cross correlation and the corresponding lag.
%
[cor lags] = vxcorr(x,y);
cmax = max(abs(cor));
nmax = find(abs(cor)==cmax);
timelag = lags(nmax);
theta = angle(cor(nmax));
y = y*exp(-i*theta); % derotate
%
% Noise BW calibration
%
hh = impz(b,a); % receiver impulse response
nbw = (fs/2)*sum(hh.^2); % noise bandwidth
%
% Delay the input, and do BER estimation on the last 128 bits.
% Use middle sample. Make sure the index does not exceed number
% of input points. Eb should be computed at the receiver input.
%
index = (10*fs+8:fs:(NN-10)*fs+8);
xx = x(index);
yy = y(index-timelag+1);
[n1 n2] = size(y2); ny2 = n1*n2;
eb = tb*sum(sum(abs(y2).^2))/ny2;
eb = eb/2;
[peideal,pesystem] = qpsk_berest(xx,yy,ebn0db,eb,tb,nbw);
subplot(1,2,1)
yscale = 1.5*max(real(yy));
plot(yy,'+')
xlabel('Direct Sample'), ylabel('Quadrature Sample'), grid;
axis([-yscale yscale -yscale yscale])
subplot(1,2,2)
semilogy(ebn0db,peideal,ebn0db,pesystem,'ro-'), grid;
xlabel('E_b/N_0 (dB)'), ylabel('Bit Error Rate')
legend('AWGN Reference','System Under Study')
% End of script file.Program random_binary.m is defined in Appendix A of this chapter. Program vxcorr.m is defined in Appendix B of this chapter. Program q.m is defined in Appendix C of this chapter.
% File: qpsk_berest.m
function [peideal,pesystem] = qpsk_berest(xx,yy,ebn0db,eb,tb,nbw)
% ebn0db is an array of Eb/No values in db (specified at the
% receiver input); tb is the bit duration and nbw is the noise BW
% xx is the reference (ideal) input; yy is the distorted output;
%
[n1 n2] = size(xx); nx = n1*n2;
[n3 n4] = size(yy); ny = n3*n4;
[n5 n6] = size(ebn0db); neb = n5*n6;
%
% For comparision purposes, set the noise BW of the ideal
% receiver (integrate and dump) to be equal to rs/2.
%
nbwideal = 1/(2*tb*2);
for m=1:neb
peideal(m) = 0.0; pesystem(m) = 0.0; % initialize
%
% Find n0 and the variance of the noise.
%
string1 = ['Eb/No = ',num2str(ebn0db(m))];
disp(string1) % track execution
ebn0(m) = 10^(ebn0db(m)/10); % dB to linear
n0 = eb/ebn0(m); % noise power
sigma = sqrt(n0*nbw*2); % variance
sigma1 = sqrt(n0*nbwideal*2); % variance of ideal
%
% Multiply the input constellation/signal by a scale factor so
% that input constellation and the constellations/signal at the
% input to receive filter have the same ave power
% a=sqrt(2*eb/(2*tb)).
%
b = sqrt(2*eb/tb)/sqrt(sum(abs(xx).^2)/nx);
for n=1:nx
theta = angle(xx(n));
if (theta<0)
theta = theta+2*pi;
end
%
% Rotate x and y to the first quadrant and compute BER.
%
xxx(n) = b*xx(n)*exp(-i*(theta-(pi/4)));
yyy(n) = yy(n)*exp(-i*(theta-(pi/4)));
d1 = real(xxx(n)); d2 = imag(xxx(n)); % reference
d3 = real(yyy(n)); d4 = imag(yyy(n)); % system
pe1 = q(d1/sigma1)+q(d2/sigma1); % reference
pe2 = q(d3/sigma)+q(d4/sigma); % system
peideal(m) = peideal(m)+pe1; % SER of reference
pesystem(m) = pesystem(m)+pe2; % SER of system
end
end
peideal = (1/2)*peideal./nx; % convert to BER
pesystem = (1/2)*pesystem./nx; % convert to BER
% End of function file.[1] Recall that a full-response system is one in which the signal energy at the receiver is constrained to the symbol interval so that the matched-filter receiver, which integrates the received signal over a symbol period, captures all of the transmitted symbol energy.
[2] Note that for QPSK we use semianalytic simulation to compute the symbol error probability, since QPSK points in signal space are defined by symbols rather than bits. Once the symbol error probability is determined, the symbol error probability can be converted to the bit error probability using analysis. For binary PSK, the symbol error probability and the bit error probability are, of course, equivalent.