UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

On the computation of the probability of undetected error for linear block codes on the Gilbert channel Wong, Brenden 1991

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


831-UBC_1991_A7 W66.pdf [ 4.68MB ]
JSON: 831-1.0098608.json
JSON-LD: 831-1.0098608-ld.json
RDF/XML (Pretty): 831-1.0098608-rdf.xml
RDF/JSON: 831-1.0098608-rdf.json
Turtle: 831-1.0098608-turtle.txt
N-Triples: 831-1.0098608-rdf-ntriples.txt
Original Record: 831-1.0098608-source.json
Full Text

Full Text

O N THE COMPUTATION OF THE PROBABILITY OF UNDETECTED ERROR FOR LINEAR BLOCK CODES O N THE GILBERT C H A N N E L by BRENDEN WONG B. Sc. University of British Columbia A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF ELECTRICAL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September, 1991 © Brenden Wong, 1991 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. 1 further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of E l e c t r i c a l Engineering The University of British Columbia Vancouver, Canada Date Oct. 15, 1991  DE-6 (2/88) Abstract An important measure of the performance of error detecting codes is the probability of undetected error. Extensive study on the subject has yielded results which allow for the computation of the probability of undetected error for many codes on the binary symmetric channel (BSC). However, little is known about code performance in more complicated channel models. The Gilbert channel is a two-state, three-parameter model with memory which simu-lates the effects of burst noise. In this thesis, we investigate methods to compute the probability of undetected error of binary linear block codes on this channel. We examine an approach to approximate code performance based on the P(m,n) distribution which is the probability of m errors in a block of n bits and the weight distribution of the code. For the Gilbert channel, P(m,n) can in prin-ciple be calculated from the channel parameters. In practice however, existing methodologies suffer from rather excessive computational requirements, particu-larly when n is larger than one thousand or so. We have developed an efficient method to calculate P(m,n) for reasonable channel parameters. This allows the probability of undetected error for many codes to be readily estimated. For certain channel and code parameters, the approximation method described above may not be sufficiently accurate. Exact analytical results are difficult to obtain, however; because unlike the BSC, the probability of a particular error pattern on the Gilbert channel depends not just on the number of i's in the pattern. Nevertheless, by appropriately exploiting certain symmetries present on the Gilbert channel, we can acquire some useful results. We have derived the probability of undetected error for the single parity check code. We ii have also obtained a formula for summing over a cyclic class of vectors and shown that reciprocal generator polynomials generate cyclic codes which have the same probability of undetected error on the Gilbert channel. The Monte Carlo simulation technique is often used when exact analysis is difficult. In a simulation study of CRC codes, we are able to observe several interesting qualitative results with just a reasonable amount of computational effort. We find that as on the BSC, on the Gilbert channel the probability of undetected error does not always increase with worsening channel conditions. Also, the CRC-CCITT code appears to maintain its superiority in terms of error detection performance over the CRC-ANSI code on the Gilbert channel, and perhaps most significantly, for some ranges of channel parameters, the probability of undetected error estimated using BSC results with the effective bit error rate can be quite inaccurate. iii Contents Abstract ii List of Tables v List of Figures vi Acknowledgement viii 1. Introduction 1 2. Approximation By Averaging 11 2.1 Using P(m,n) to Approximate P u 11 2.2 A Series Expansion for P(m,n) 13 2.3 Comparison with Other Methods of Computing P(m,n) 24 2.4 Some Observations on Short BCH Codes 29 3. Some Exact Results 34 3.1 No Coding 34 3.2 Single Parity Check Encoding 34 3.3 Probability of a Vector and its Cyclic Shifts 40 3.4 Reciprocal Generator Polynomials 56 4. Monte Carlo Simulation 58 4.1 The Monte Carlo Method 58 4.2 Simulation Study of CRC-16 Codes 59 5. Conclusions 73 References 75 Appendix A 77 Appendix B 80 Appendix C 83 iv List of Tables Table 2.1 Computation of P(20,n) using (2.6) and double precision C language routines 16 Table 2.2 Efficiency of series expansion technique for computing P0(m,n) ...22 Table 2.3 Comparison with Cuperman's expression for P0(m,n) 25 Table 2.4 Comparison with Cuperman's expression for P(m,n) 26 Table 2.5 Comparison of computational times of P(m,n) for the recursion and series expansion methods 28 v List of Figures Figure 1.1 The Gilbert channel 5 Figure 3.1 P u for single parity check code, n=8, h=.5 41 Figure 3.2 P u for single parity check code, n=8, h=.7 42 Figure 3.3 P u for single parity check code, n=8, h=.9 43 Figure 3.4 P u for single parity check code, n=16, h=.5 44 Figure 3.5 P u for single parity check code, n=16, h=.7 45 Figure 3.6 P u for single parity check code, n=16, h=.9 46 Figure 3.7 P u for single parity check code, n=32, h=.5 47 Figure 3.8 P u for single parity check code, n=32, h=.7 48 Figure 3.9 P u for single parity check code, n=32, h=.9 49 Figure 3.10 P u for single parity check code, n=64, h=.5 50 Figure 3.11 P u for single parity check code, n=64, h=.7 51 Figure 3.12 P u for single parity check code, n=64, h=.9 52 Figure 4.1 Simulation of CRC-16 codes, p=.01, h=.5 62 Figure 4.2 Simulation of CRC-16 codes, p=.01, h=.7 63 Figure 4.3 Simulation of CRC-16 codes, p=.01, h=.9 64 Figure 4.4 Simulation of CRC-16 codes, p=.l, h=.5 65 Figure 4.5 Simulation of CRC-16 codes, p=.l, h=.7 66 Figure 4.6 Simulation of CRC-16 codes, p=. 1, h=.9 67 vi Figure 4.7 Simulation of CRC-16 codes, p=.3, h=.5 68 Figure 4.8 Simulation of CRC-16 codes, p=.3, h=.7 69 Figure 4.9 Simulation of CRC-16 codes, p=.3, h=.9 70 vii Acknowledgement First, I must thank Prof. C. Leung for his very excellent supervision throughout the course of this work and also for his help in obtaining financial support. I must also thank my peers from whom I have benefitted greatly by asso-ciation. Space restrictions prevent me from naming them all so I should just mention that several discussions with C. Lai and V. Wong on the subject on this thesis have been most helpful. Finally, I gratefully acknowledge the financial support of NSERC through Grant #OGP0001731 and The Department of Electrical Engineering which were in the form of partial research assistantship and partial scholarship respectively. viii 1. Introduction The ability to exchange or retrieve information accurately and efficiently is of paramount importance in everyday life. The speed and accessibility requirements of today have led to complex communications architectures, enough so that it is customary to consider a layered structure where each layer serves as a virtual communications link. At the most basic level (the "physical layer"), the link is physically real and consists of a data source, an encoder, a channel over which the data is transferred, a decoder, and the final output as the receiver's interpretation of the original message. Clearly, correct and efficient operation at this level is essential and ultimately related to the performance of a complex hierarchical system. The physical communications channel may be thought of as being described by the medium over which data is transmitted or stored and the method by which such data is represented (modulation). The problem is that in general, due to various imperfections, the channel distorts the original message. It is the purpose of coding to add redundancy in a well defined manner to the source data in order that discrepancies from the original message may be detected and/or corrected by the receiver. Coding is a mathematical scheme which maps each possible message con-sisting of a sequence of symbols into a distinct sequence of symbols known as a codeword [1-3]. In practice, the sets of message and codeword symbols, that is, the message and codeword alphabets, are often identical. Naturally, the most common alphabets used for digital communications consist of powers of two symbols so the binary alphabet {0,1} is one of particular significance. In a block code, the number of symbols in a message or codeword is fixed. If the messages 1 7. Introduction 2 are of length k and the alphabet is of size s, there are sk k-tuples representing all possible messages. The k-tuples spanning the space are mapped one to one into the space of n-tuples V n . A linear mapping (with respect to addition and multiplication in the finite field GF(s)) is particularly simple to visualize and implement as it can be accomplished by matrix multiplication. A linear block code C(n,k) is then defined to be a k-dimensional linear subspace of V n . For a message m (k-dimensional row vector), the corresponding codeword c (n-dimensional row vector) is given by c=mG where G (kxn matrix with rank k) is known as the generator matrix. The extra n-k symbols are called parity symbols. It is often convenient to represent a codeword by a polynomial. That is, we identify the coordinate representation (co,Ci,...,cn.i) of a codeword c with the polynomial c(x)=c0+Cix+...+c11.1xn"1. One particular class of linear block code which has been much studied and often used in practice is the class of cyclic codes. A cyclic code has the property that if a given vector is a codeword, then so are all vectors which arise from cyclic permutations of this vector. A cyclic code can be characterized by a generating polynomial g(x) of degree n-k. The codewords of a cyclic code are then given by the set of all polynomials (g(x)d(x)} where d(x) is a polynomial with degree at most k-1. Let cGC(n,k) be a transmitted codeword. The received vector v € V n may be written as v=c+e where e G V n is the error vector representing the distortions or noise introduced by the channel. Hence, the ith component of e is a 0 or 1 depending on whether the ith bit transmission is error free or in error respectively. In general, e^O and v£C(n ,k ) ; this being the case the receiver has several recourses. If the code is used for error detection, v is discarded and a retransmission is requested (the procedure by which this occurs is the function 7. Introduction 3 of the second layer in the comrnunications hierarchy, the "data link layer", with which we will not be concerned here). It is also possible that the received vector v £ C ( n , k ) but v^c in which case the code has failed to detect the error. The probability of undetected error is then a measure of reliability of an error-detecting code. If the code is used for error correction, an estimate c' 6 C(n,k) is found based on v and the maximum likelihood principle. An error in decoding results if the estimate is incorrect, that is if cVc . Thus the probability of decoding error is a measure of the ability of an error-correcting code to identify correctly a codeword that is corrupted by noise. A code can also be used for simultaneous error detection and correction, that is to decode received vectors which are reasonably similar to corresponding codewords and reject those which are too dissimilar. In general, error correction is much more difficult to analyze and implement than error detection. In the following we consider only binary linear block codes used for error detection. One of the simplest noise models is the binary symmetric channel (BSC). In this channel, two possible symbols, 0 and 1, or bits in this case, are transmit-ted and received. There is a probability e that a transmitted bit will be received incorrectly. Successive errors are independent. The bit error rate e is assumed to be constant and representative of the overall channel conditions. For a linear code, the probability of undetected error is independent of the particular code-word transmitted and is just the probability that the error vector is itself a code-word. The Hamming weight (or simply weight) of a vector v, denoted wt(v), is the number of l's that it contains. If the number of codewords with Hamming weight i is Aj (i=0,...,n), the Hamming weight enumerator is the polynomial defined by 7. Introduction 4 A(z)= I Aiz1 . (1.1) i=0 For the BSC, all error vectors of weight i have probability e^l-e)11'1 so that knowledge of the weight distribution of a code is sufficient for computing its probability of undectected error since P u = prob(e*0, eeC(n,k)) = I AieKl-e)11-1 i=l = (l-e)n[A(e/(l-e))- 1] . (1.2) There is a large body of literature describing many results and techniques used in the study of weight enumerators [1-3]. Regarding P u itself, an important characteristic is whether it is monotonic with respect to e. A code where P u is monotonically increasing with e is called proper. Properness and other aspects of P u are discussed for various codes on the BSC in [4-11]. While the BSC is useful for describing channels afflicted by random or shot noise errors, its memoryless nature renders it inappropriate for channels where errors tend to occur in clusters or bursts. Burst noise effects are typically associated with multipath fading in radio channels, switching transients and crosstalk on wired channels, and surface defects and foreign particles on mag-netic storage media. The Gilbert channel, illustrated in Fig. 1.1, is a model devised to simulate such burst effects [12]. The channel consists of two states, a "good" state G where transmission is completely error-free and a "bad" state B where the probability of correct transmission is h; a state transition from G to B occurs with probability P and from B to G with probability p. It is convenient to define Q=l-P, q=l-p, and h'=l-h so Q is the probability of remaining in state G, 7. Introduction 5 P Figure 1.1 The Gilbert channel J. Introduction 6 q is the probability of remaining in state B, and h ' is the bit error rate while in state B . As well, we denote the average or unconditional probability of the system being in state B by PB=P/(P+P), the corresponding probability for state G by p G = p / ( P + p ) , and the effective bit error rate of the channel by e e f f = p R h ' . For suitable (small) values of P and p, the model produces the type of burst error statistics that are of interest. The idea is that most of the time the channel is in the good state but occasionally, due to a change in the transmission characteristics, the channel lapses into the bad state which tends to persist. As an example, Gilbert found an adequate fit using h=0.84, P=0.003, p=0.034 for Call #1296 which was over a 500-mile radio path with loaded cables at the ends from the telephone circuit measurements of Alexander, Gryb, and Nast [13]. An extension to the Gilbert channel where G is not completely error-free but has a probability of correct transmission g ( g >h ) was considered by Elliott [14] and hence is known as the Gilbert-Elliott channel. As the Gilbert channel is a system with memory, the probability of the next error vector bit being 0 or 1 depends on the current state of the channel and therefore on the past history of bits. However, this dependence can extend back only as far as when the last 1 was received since that pins down the system in the B state at that instant. Alternatively, the conditional independence can be expressed by prob(e I ...20i)=prob(e 110) where e=0,l and the superscript j=0,l,... denotes the number of repetitions. If w{j)=prob($i), v(j)=prob(0>2 11), and u(j)=prob($ 11), then the probability of a particular error vector e can be decomposed [14] as l e l - l prob(e) = w(a) j [~[ v(bO \ u(c) (1.3) i=l 7. Introduction 7 where a is the number of O's before the first 1, b\ is the number of O's between the ith and (i+l)th 1, and c is the number of O's following the last 1 in e. The above "building block" probabilities in the Gilbert channel can be drawn directly from the study of recurrent events in probability theory [15]. The results are u(j) = [(J+p-Q)Jj - (L+p-Q)U] v(j) = j £ [(qJ+p-Q)JJ - (qL+p-Q)U] w(j) = pBu(j) (1.4) where J = |{Q+hqW(Q+hq)2+4h(p-Q)} L = |{Q+hq-V(Q+hq)2+4h(p-Q)} (1.5) are the roots of the quadratic D(t)=t2-(Q+hq)t-h(p-Q). Generally, in order to find the probability of undetected error of a C(n,k) code in the Gilbert channel, we need to sum over 2 k -l expressions of the form given in (1.3). Consider the computational effort required for such an exhaus-tive evaluation technique. To generate a codeword e (using matrix multiplica-tion) requires 2nk additions and multiplications. Assuming that the u(-), v(-), and w(-) are precomputed, to compute prob(e) requires wt(e)+l multiplications. Then overall, at least approximately ( l + 2 n k + n / 2 ) 2 k « n k 2 k + 1 arithmetic operations are required to obtain P u for each set of parameters P, p, h. We see that the number of such operations and therefore the computational complexity grows rapidly with k. For even a modest size code, say the Hamming (31,26) code, the actual computation using a Sun SPARC® station 1 workstation takes approximately 4 days. As a matter of comparison, the MATLAB® benchmark /. Introduction 8 tests rate the SPARC® station 1 at 1.4 Mflops (millions floating point instructions per second) and the Cray XMP® supercomputer at 33 Mflops. The codewords themselves can be precomputed and stored to effect a small savings of factor 4k, but storing all codewords requires 2kn bits of media space which for the above example of the Hamming (31,26) code is approximately 260 Mbyte. Clearly, determining P u by exhaustive computation for all but the shortest codes is impossible even with the most powerful computing facilities that are currently available. The main difficulty is the exponential increase in the number of codewords with k. For the BSC, the situation is alleviated by the natural divi-sion of the codewords into n weight classes since all codewords with a given weight have the same probability and that is why the weight distributions of block codes have been studied extensively. For the Gilbert channel however, knowledge of a code's weight distribution is of no direct value; the general prob-lem of determining P u is much more difficult. Nevertheless, we have already alluded to the importance of coding and the fact that channel conditions are often not representable by the BSC alone so this problem deserves investigation. Our objective here will be to find viable means of computing P u. Two binary linear block codes are said to be equivalent to each other if the generator matrix of one can be obtained by a column (coordinate label) permuta-tion of any generator matrix of the other [3]. Evidently, this definition is motivated by the fact that, on a BSC, equivalent codes have identical error detecting (and error correcting) properties since they have exactly the same weight enumerators. In general, the average probability of undetected error E[PU(C)] over all equivalent codes of a given code C(n,k) depends only on the weights Ai and the probability P(m,n) of m errors in a block of n bits. E[PJ may 7. Introduction 9 be considered an approximation to PU(C) [14]. This approach appears feasible for the Gilbert channel as P(m,n) may be obtained from recursion relations given by Elliott [14,16] or a closed form expression given by Cuperman [17]. In practice, however, the computation is problematical, especially if n is larger than one thousand or so. We address these problems in Chapter 2 by developing a technique that computes P(m,n) efficiently for reasonable channel parameters. As well, we make some observations regarding the validity of the E[PJ approximation and the performance of some short cyclic codes of the BCH (Bose-Chaudhuri-Hocquenghem) type relative to their equivalents. There has also been an approach to approximating P u using the burst length distribution Q(l,n) where 1 is the burst length defined as the number of bit positions from the first to last error inclusive although it does not appear to be applicable for many realistic channel conditions [18]. In general, it is difficult to find exact analytical expressions for P u . In fact, no nontrivial result of this type is known. A single parity check code is one of the simplest linear block codes; an extra bit which is the modulo 2 sum of the message bits is appended to the original message to form a codeword. Chapter 3 consists primarily of the derivation of P u for the single parity check code. We also give there an expression for the probability of a given vector and its cyclic permutations which may be useful in determining P u for cyclic codes as well as an identity regarding the P u of cyclic codes which have generator polynomials satisfying a particular relationship. A commonly used technique for determining error rates in communica-tions systems is the Monte Carlo Simulation [19]. Random noise effects are com-puter generated according to some fixed distribution corresponding to system 7. Introduction 10 parameters and error events are tabulated. The average number of such events over a large number of samples can then be taken as an estimate of the actual physical error rate. As an application of the Monte Carlo technique, we study the performance of two CRC (Cyclic Redundancy Check) codes in the Gilbert Channel. A CRC codeword consists of the message bits plus a number of parity bits which corresponds to the remainder of a modulo 2 division of the message expressed as a polynomial and a predetermined polynomial known as the CRC generator polynomial. A judiciously chosen generator polynomial endows the CRC code with some desirable error detection characteristics. In particular, CRC codes are known for their strength in detecting burst errors. Hence, it is reasonable and useful to examine CRC codes under burst channel conditions in order to compare with and augment the conventional BSC results. In Chapter 5, we give a summary of our findings and discuss topics which may deserve further study. Two appendices give some necessary details for the computations in Chapter 2. A third appendix contains the source code listings for the programs used in the numerical computations and simulations. Maple is a registered trademark of Waterloo Software Systems. M A T L A B is a registered trademark of The Math Works, Inc. Sun SPARC is a registered trademark of Sun Microsystems. Cray XMP is a registered trademark of Cray Research, Inc. 2. Approximation By Averaging 2.1 Using P(m,n) to Approximate Pu If P(m,n) is the probability of m 2's in an error vector of length n bits, then the average probability of a weight m error vector is P(m,n)/cm where cm=n!/[m!(n-m)!] is the binomial coefficient. Hence an approximation to the probability of undetected error of a C(n,k) code is given by The idea is that even though codewords within the same weight class may in general have rather different probabilities, collectively these differences tend to cancel so if A m is quite large, we can expect the total probability of all weight m codewords to be reasonably approximated by A m times the average probability of a weight m vector. Let TC be a permutation of n indices. Denote the equivalent code resulting from the application of K to the vectors of C(n,k) by 7tC(n,k). Then as noted by Elliott [14], 'm (2.1) m=l 11 2. Approximation By Averaging 12 E[PJ = ^ I Pu(7tO 7t = T7 £ E prob(x) * rt x^OxercC = ~7 E E prob(Trx) * X ^ O x S C TC = 37 E [wt(x)]![n-wt(x)]!P(wt(x),n) " x ? i O x € C = £ (Ai/c?)P(i,n) (2.2) i=l so the approximation for P u in (2.1) is in fact also the exact average over equiva-lent codes. Unfortunately, there is no a priori way of determining the effectiveness of the approximation (2.1). We can obtain some empirical results only by experimentation. However, because of the connection set forth by (2.2), at least we also get information regarding the relative performance of codes as com-pared to their equivalents. That is, if the exact P u for a code were known, comparison with (2.1) or (2.2) in the first context gives a measure of the accuracy of the approximation and in the second context a measure of how well the code performs in error detection relative to its equivalent codes. From a practical viewpoint though, we must first be able to compute P(m,n) efficiently. The P(m,n) distribution, also known as the counting distribution, has been studied for a variety of channel models [20] and is often parametrized by experimentally measured statistical quantities. For the Gilbert channel however, it can be expressed directly in terms of the fundamental channel parameters. 2. Approximation By Averaging 13 2.2 A Series Expansion for P(m,n) A renewal process is a process where the occurrence of an event resets the state of the system, or in our present terminology, one where an error deter-mines the state of the channel. The P(m,n) distribution for a renewal process is given quite generally by the recursion relations [16] P(0,n) = 1 - E w(i) i=l n-m+1 P(m,n)= E Piu(i-l)R(m,n-i+l); l < m < n (2.3) i=l where R(l,n) = u(n-1) n-m+1 R(m,n)= E v(i-l)R(m-l,n-i); 2<m<n , (2.4) i=l u(-), v(-), and w(-) have the same meaning as the definitions given before (1.4), and pi is the unconditional bit error rate (eeff for the Gilbert channel). Computing P(m,n) this way requires a number of summations increasing as n 4 . The main difficulty however is the amount of memory required to store R(m,n) in a typical computational algorithm. For instance, a 1000x1000 array of numbers of type double in the C programming language (providing approximately 15 decimal place accuracy) requires about 8 Mbyte. Hence, in practice it would be awkward to handle n much larger than 1000 or so. Moreover, it is difficult to infer directly any relationship between P(m,n) and the channel parameters using this recursive technique. This would be a disadvantage, for instance, in the reverse problem where channel parameters 2. Approximation By Averaging 14 must be extracted from experimental distributions for then an analytic form is almost certainly preferred. A set of recursion relations for computing P(m,n) specific to the Gilbert channel requiring order n 2 summations is given in [14]. See (3.3), (3.4). However, the memory requirement is double that of the general recursive technique. Also, for n small, the computational time is not much different and is in fact inferior to that of the method which we will propose. Hence we do not discuss the latter recursion method further in this chapter. A nonrecursive approach to computing P(m,n) for the Gilbert channel was given by Cuperman [17]. If P0(m,n) is the probability that the system is in the B state m times out of n, then P(m,n) is just P(m,n) = (l-h)m £ cih^PoG,!!) (2.5) It can be shown that Po(m,n) = <j P + p Q n l ; m=0 p m - l 2 B T r Z I Vm+tC? q n > J - l Q n - m . j . t ( p _ q ) i + t . l < M < N . 1 r ~ t * 1 j=0 t=0 p q n l ; m=n . P+p (2.6) See Appendix A for details concerning the derivation of (2.6). The second equation has been written in a slightly more compact form than the corresponding equation in [17]. However, note in particular that its range of validity is actually l<m<n- l rather than l<m<n/2 as stated in [17]. On the surface, (2.5) and (2.6) appear to provide quite an efficient means of computing 2. Approximation By Averaging 15 P(m,n) with order n 3 summations. This is not so however as the evaluation of P0(m,n) is actually quite problematical. The difficulty is in the amount of numerical accuracy required. The data in Table 2.1 illustrates the point. We have tabulated there each of the m terms, s(j) for j=0,...,m-l, in the summation for P0(m,n), calculated using double precision C language routines. The state transition probabilities are fixed at P=.001, p=0.1 and m=20 is taken. Note that for ease of reading, only 7 significant digits are displayed; the actual computations in C are performed to the full 15 digits accuracy. The sums of s(j) over j are compared with the results obtained using the infinite precision computation of Maple®. For n=30, agreement is quite precise. For n=40, the deviation is slight. For n=50, however, the C program generated result is a nonsensical one. The problem is that in general, each term s(j) in the summation for P0(m,n) is large (absolute value much greater than one) but such large numbers must add together to give a small (positive) number (much less than one). In the above case n=50, the largest s(j) are of the order 101 3 while the actual P0(m,n) is of the order 10"4 so a posteriori, it is evident that at least 17 decimal digits accuracy is required for proper computation. Similarly, for the case n=40, we see that only 2 or so digits of accuracy can be expected with C routines. With n=500 say, hundreds of digits of accuracy are typically required. This order of numerical accuracy is supported by specialized computing packages such as Maple® though the drawback is extremely long execution time. Using a SPARC® station 1, the computation for a single P0(m,n) when n « 5 0 0 can take up to one hour. Here we develop a method to circumvent the problems that are encoun-tered by the previous methods. For most channel conditions which the Gilbert 2. Approximation By Averaging 16 P=0 .001 , p=0.1 j s(j) n=30 n=40 n=50 0 1.924937e+04 5.373783e+07 1.270440e+10 1 -1 .293079e+05 -5 .322397e+08 - 1 . 4 9 4 8 8 7 e + l l 2 3.840047e+05 2.436174e+09 8 . 2 2 5 3 2 3 e + l l 3 -6 .627478e+05 -6 .835843e+09 -2 .811022e+12 4 7.347007e+05 1.315219e+10 6.684053e+12 5 -5 .453716e+05 -1 .838674e+10 -1 .173890e+13 6 2.737843e+05 1.930997e+10 1.577618e+13 7 -9 .148801e+04 -1 .553482e+10 -1 .658501e+13 8 1.940020e+04 9.675176e+09 1.382414e+13 9 -2 .346285e+03 -4 .683027e+09 -9 .204888e+12 10 1.221875e+02 1.758559e+09 4.909311e+12 11 -0 .000000e+00 -5 .086999e+08 -2 .093842e+12 12 0.000000e+00 1.119024e+08 7 . 0 9 9 2 2 9 e + l l 13 -0 .000000e+00 -1 .835178e+07 - 1 . 8 9 2 8 7 5 e + l l 14 0.000000e+00 2.179734e+06 3.902036e+10 15 -0 .000000e+00 -1 .797807e+05 -6 .061153e+09 16 0.000000e+00 9.669907e+03 6.820174e+08 17 -0 .000000e+00 -3.070471e+02 -5 .213714e+07 18 0.000000e+00 4.822956e+00 2.405282e+06 19 -0 .000000e+00 - 2 . 5 1 1 7 9 4 e - 0 2 -5 .013373e+04 E s(j) 3 .934083e-04 5 .258999e-04 - 2 . 9 1 0 6 6 2 e - 0 3 Actual P(20,n) 3 .934082e-04 5 .302741e-04 6 .672299e -04 Table 2.1 Computation of P0(20,n) using (2.6) and double precision C language routines. 2. Approximation By Averaging 17 channel is designed to reasonably simulate, P^O .01 and P/p^0.1. That is, P is nominally a small parameter so we can express (2.6) as a power series in P. Then P(m,n) can be computed from (2.5) using an approximation for P0(m,n) that is its truncated power series, the number of terms taken corresponding to the accuracy desired. The point is that successive terms in the expansion get progressively smaller so that severe numerical precision is not required. In order to derive a compact expression for the order r contribution to Pn(m,n), it is simpler to start from an alternate form of (2.6). Since P0(m,n)(P,p) = prob(system in state B m times out of n) = prob(system in state G n-m times out of n) = prob(system in state B n-m times out of n with state transition probabilities reversed) = P0(n-m,n)(p,P) , (2.7) the second line in (2.6) can be written n-m-l 2 P0(m,n) = P G £ £ cnT_1 cm-^  °t Q^V^ P - q ) 1 i+t j=0 t=0 2 n-m-1 = P c q m £ cK-lf £ c-F1 (-1)* t=0 j=0 j+t n-m-j-1 • L E o>£C^ Y1 (-P) a +v (2.8) a=0 0=0 Denote the coefficient of the P r term in (2.8) by Pn(m,n;r), that is, write n-m+1 P0(m,n) = p G I P0(m,n;r) P r , (2.9) r=0 then 2. Approximation By Averaging 18 n-m-l P0(m,n;r) = (-l)rqm-r £ c2(-l)1 £ c ^ f 1 c m ^ (-1)1 £ cj+t c n-m-j - i q P = (-D rqm- r £ q p t CtC-D4 P = o t = o n-m-l Vj ( . ^ c n ^ l c n - t ^ t c n - m - j - l I j=0 (2.10) The last expression can be simplified using the results in Appendix B. For t=2, the factor in braces is n-m-l £ ( - l ^ C n - m - l ^ ^ c n - m - j - l j=0 — 1 (n-2-,j)! (i+2)(i+l) (m-2-j)!(n-m) (r-p)![j+2-(r-P)]! p!(n-m-j-l-(3)! j = o ^ n-m-l , . - (n-m)R!(r-P)! 1 i y [n-(n-m)-2-j]![j-(r-P-2)]![(n-m)-P-j-l]! j = o . {[j-(r-P-2)]|j-(r-P-l)] + 2(r-P)[j-(r-P-2)] + (r-P)(r-P-l)} = (n-m)PKr-P)! (U(n,n-m,2,r-p\P) + 2(r-P)U(n,n-m,2,r-P-l,P) + (r-p)(r-p-l)U(n,n-m,2,r-p-2,p)} . (2.11) where the function U(-,•,-,-,•) is as defined in (B.l). Similarly, for t=l and t=0, we have n-m-l Vj ( . l y c n - m - l c n - ^ ^ l c n - m - j - l j=0 = (n-m)p!(r-P)! {U(n+l,n-m,2,r-P,P) + (r-P)U(n+l,n-m,2,r-p-l,P)} . (2.12) and 2. Approximation By Averaging 19 g \ . 1 ) i c n - m - l c ^ C r J p C n . r n - j - l j = 0 = (n-m)p!(r-P)! {U(n+2,n-m,2,r-p,p)} (2.13) respectively. Hence, summing over t in (2.10) gives P o ( m ' n ' r ) = (- 1 ) r^ r p? 0^(n-m)p!(r-P)! . {[U(n+2,n-m,2,r-p,p) - 2U(n+l,n-m,2,r-p,P) + U(n,n-m,2,r-p,P)] + 2(r-p)[-U(n+l,n-m,2,r-P-l,P) + U(n,n-m,2,r-p-l,P)] + (r-P)(r-P-l)U(n,n-m,2,r-P-2,P)} . (2.14) The fourth argument of the U(v>V>') function must be nonnegative in order to use the results of Appendix B so the summand in the last equation has to be treated in three parts. Using (B.5), (B.10), (B.ll), (2.14) becomes after some algebra P0(m,n;r) = ( - l ) r q m r 7 ^ ) rA2 „ (n-m-p+l)(n-m-p) * ] pV"q ^ ( n + P> n - m ' 2 > r - 2 > -r(n-m-r+2)  + ^ ( m + r . i X r - l ) ! V(n+r,n-m,2,r-l) L,n-m,2,r) j + ( " q ) r m(mtr)r! V(n+r+l,n-m,2,r) \ (2.15) for r>2 where the function V(v>v) is defined in (B.2). For r=l, the same ex-pression holds except that the sum should be ignored. For r=0, only the last term inside the braces is present. Thus 2. Approximation By Averaging 20 Po(m,n;0) = 0 P0(m,n;l) = q 1 ^ 1 ^ {(n-m+l)V(n+l,n-m,2,0) + ^V(n+2,n-m,2,l)} P0(m,n;r) = ( - I W T ^ ) £ „ (n-m-B+l)(n-m-B) £ ( ^ ) P ^RMRi t LV(n+P,n-m,2,r-2) • ; r > 2 . (r-B)!B! (2.16) Evaluating V(-,•,-,-) by means of (B.9), the three equations can in fact be rewritten as one. For consistency, we also write the first and third equations of (2.6) as series expansions. The final result is P0(0,n;r) = (-l) rcn;1 1 v» (n-m-B+l)(n-m-B) Po(m,n;r) = q™ — c ™ L m ^ P0(n,n;r) = q*"1 p 1 5r>0 c^cmr+P(-q)P ; l < m < n - l (2.17) where 5r/0 equals 1 if r=0 and 0 otherwise. To summarize, what we have accomplished thus far is to give a means of calculating P0(m,n) through the use of (2.9) and (2.17). However, the summation in r can be extended only as far as necessary depending on the precision required (and the available computational resources). That is, we can truncate (2.9) and take the partial sums as successively better approximations: P0(m,n) = p G t P0(m,n;r)Pr = P0(R)(m,n) (2.18) r=0 The subsequent computation of P(m,n), the quantity which we are ultimately interested in, is straightforward as (2.5) involves a relatively small (average n/2) 2. Approximation By Averaging 21 number of positive terms so that great numerical accuracy is not required. The number of summations required when determining P(m,n) this way is of order n 2 r 2 . Clearly, if the P0(m,n) are known to within a certain percentage accuracy, P(m,n) computed from (2.5) can be no worse. The efficiency of our technique is summarized in Table 2.2. The entries in the table represent the number of terms in the expansion required in order to compute P0(m,n) to within 0.1% of the exact value. As expected, the method works best for small P since P is the expansion parameter. Also, the number of terms required increases with p and n." This dependence can be explained as follows. Keeping only terms up to order R in P corresponds approximately to discarding processes with more than R G-»B state transitions. A flip-flopping of states involves the factor pP so when p is large, processes with relatively large numbers of G-»B transitions can less likely be ignored. Similary, if n is large, there is more opportunity for the system to make a G-*B transition. Some remarks concerning how the table is compiled are in order. The ratios IP0(m,n;r+l)Pr + 1 l/lP0 ( r )(m,n)I are computed for all m for successive values of r. When for a particular value of r=R, the ratio becomes less than 0.1% for all m, the computation is stopped and R is entered into the table. We have checked extensively using the recursion techique (2.3) for n< 511 and selectively using the exact Maple® calculation (2.6) for n>511 that indeed {I P0(m,n)-P0(R)(m,n) l/P0(m,n)} £0 .1%. This justifies our choice of 0.1% both as a reasonable level of accuracy and as a "convergence" criterion. However, one caveat on the use of our method is pointed to by parenthesized entries in Table 2.2. For those entries, greater than 15 digits precision is required to compute the individual terms P0(m,n;r) for r close to R so that Maple® must be used. 2. Approximation By Averaging 22 n P=.00001 p=.0001 p=.001 p=.01 P=.l D=.3 7 1 1 1 1 1 15 1 1 1 1 1 31 , 1 1 1 1 1 63 1 1 1 1 2 127 2 2 2 2 2 255 2 2 2 2 2 511 2 2 2 2 3 1023 2 2 2 3 4 2047 2 2 2 4 6 4095 2 2 3 6 10 n P=.0001 p=.001 p=.01 P=.l p=.3 7 1 1 1 1 15 2 2 2 2 31 2 2 2 2 63 2 2 2 2 127 2 2 2 3 255 2 2 3 4 511 3 3 4 5 1023 3 3 5 8 2047 4 4 8 13 4095 5 6 (13) (21) n P=.001 P=.01 p=.01 p=.l p=.3 P=.l p=.3 7 2 2 2 2 2 15 2 2 2 3 3 31 2 2 2 4 4 63 3 3 3 6 6 127 3 3 4 8 8 255 4 4 7 13 13 511 5 6 10 (22) (22) 1023 7 10 17 (40) (41) 2047 ( ID (17) (30) >50 >50 4095 (19) (29) >50 >50 >50 Table 2.2 Efficiency of series expansion technique for computing P0(m,n). 2. Approximation By Averaging 23 From (2.17), we can see that once again the problem is the appearance of large numbers with alternating signs in the sum. The dominant factor in the summand is cm+p which for large m will give rise to r+1 large numbers unless r is quite small. For example, with n=2047, m=1023, p=.01, and r=l l , the summand ranges from 3.2xl0 2 8 to 1.4xl031 in absolute value. Note however that the problem may not arise for certain combinations of parameters. In particular, if P is small enough, PrP0(m,n;r) may in fact be negligibly small for large r. Also, if the accuracy requirement is made somewhat less stringent, then fewer terms will be needed in the expansion so that the numerical problem can be delayed from occurring. Hence, the resummation of (2.6) leading to our present formalism has accomplished what we expected in that the terms in the expansion (2.9) do indeed become small rapidly for reasonable channel parameters although the numerical accuracy problem has not entirely disappeared. We should also mention that since pG=( 1+P/p)"1, one might think that a more consistent technique would also take into account the expansion of this factor in (2.8) and (2.9). In fact, we find that doing this results in no improvement, neither in accuracy nor efficiency. Rather, it just serves to complicate matters as the sum in (2.18) becomes a double sum, and the expansion corresponding to (2.9) becomes infinite. 2. Approximation By Averaging 24 2.3 Comparison with Other Methods of Computing P(m,n) Cuperman gives the following approximations based on (2.6) together with their claimed ranges of validity [17], P0(m,n) & Pqm-1Qn-2m[(n-m-l)p+2] ; m P « l , P « p , q P(m,n) & P(l-h)m t h imq i l[(n-i-l)p+2] ; n P « l , P « p , q . (2.19) i=m Both expressions are meant to be used for l<m<n/2. The ad hoc manner in which (2.19) is derived does not allow for consistent interpretation or improvement. Table 2.3 gives a comparison of approximations for P0(m,n) using the first equation in (2.19) and our series expansion method with (2.9) and (2.17). For n=256, P=.0001, and p=.3, we take R=4 to give us the approximation P0(4 )(m,n) which is within 0.1% of the exact P0(m,n). As can be seen from the table, Cuperman's approximation deviates by as much as 28%. Even for m=50 so that mP=.005, the disagreement is still 19%. Generally, the second equation in (2.19) is even more inaccurate. With channel parameters P=.0001, p=.l, h=.7, even with a block length as small as n=16 (so nP=.0016), Cuperman's approximation is off by as much as 55% as Table 2.4 shows. There, P (2)(m,n) is the series expansion approximation calculated from (2.5) using Po(2)(m,n) and is again within 0.1% of the exact value of P(m,n). Note that only two terms in the expansion for P0(m,n) are required to reach that level of accuracy. The computation of P(m,n) by the recursion relations of (2.3) and (2.4) is in general reasonably efficient for n £ 1000. We have already mentioned the 2. Approximation By Averaging 11=256, P=.0001, p=.3 m P0(m,n) (Cuperman) Po(4)(m,n) % difference 1 7 .62e -03 7 .62e -03 0 .0 5 1 .80e-03 1.84e-03 2 . 1 10 2 . 9 8 e - 0 4 3 .12e -04 4 . 6 15 4 . 9 1 e - 0 5 5 .27e -05 6 .9 20 8 .09e -06 8 .89e -06 9 . 0 25 1 .33e-06 1 .50e-06 1 1 . 0 30 2 . 1 9 e - 0 7 2 . 5 2 e - 0 7 12 .8 35 3 .61e -08 4 . 2 2 e - 0 8 1 4 . 5 40 5 .94e -09 7 . 0 8 e - 0 9 16 . 0 45 9 . 7 7 e - 1 0 1 .18e-09 1 7 . 5 50 1 .61e-10 1 .98e-10 18 .8 55 2 . 6 4 e - l l 3 . 3 0 e - l l 2 0 . 0 60 4 . 3 3 e - 1 2 5 .49e -12 2 1 . 1 65 7 .11e -13 9 .13e -13 2 2 . 1 70 1 .17e-13 1.51e-13 23 .1 75 1.91e-14 2 .51e -14 2 3 . 9 80 3 . 1 3 e - 1 5 4 . 1 5 e - 1 5 2 4 . 6 85 5 .11e -16 6 .85e -16 2 5 . 3 90 8 .36e -17 1 .13e-16 2 5 . 9 95 1 .37e-17 1 .86e-17 2 6 . 4 100 2 . 2 3 e - 1 8 3 .05e -18 2 6 . 8 105 3 .63e -19 4 . 9 9 e - 1 9 2 7 . 2 110 5 .92e -20 8 .16e -20 2 7 . 5 115 9 . 6 3 e - 2 1 1 .33e-20 2 7 . 7 120 1 .56e-21 2 . 1 7 e - 2 1 2 7 . 9 125 2 . 5 4 e - 2 2 3 .52e -22 2 7 . 9 128 8 .52e -23 1 .18e-22 2 7 . 9 Table 2.3 Comparison with Cuperman's expression for Po(m, 2. Approximation By Averaging n=16, P=.0001, p=.l h=.7 m P(m,n) (Cuperman) P<2)(m,n) % difference 1 6 .69e -04 6 .72e -04 0 .4 2 4 . 4 0 e - 0 4 4 .52e -04 2 . 6 3 2 .80e -04 3 .05e -04 7 . 9 4 1 .68e-04 2 .02e -04 1 6 . 8 5 9 . 2 0 e - 0 5 1.27e-04 2 7 . 5 6 4 . 4 6 e - 0 5 7 . 2 1 e - 0 5 3 8 . 1 7 1 .88e-05 3 .56e -05 47 . 3 8 6 .68e -06 1 .48e-05 5 4 . 8 Table 2.4 Comparison with Cuperman's expression for P(m,n). 2. Approximation By Averaging 27 excessive memory requirement, however. Note that reducing the precision in order to conserve memory in this method is risky because of the large number of terms involved. Consider Table 2.5 which is a fairly representative example comparing the times required to compute all P(m,n) for n=1000 with P=0.0001, p=0.1, h=0.5 using the recursion and series expansion methods on a SPARC® station 2. When the series expansion computations can be handled by C language programs, the improvement in computational time is substantial. Even when Maple programming must be used, times to compute the P0(m,n) are comparable to that of the recursion technique in computing P(m,n) (see remark below), especially when we take only r=5, which will nonetheless produce results to within 0.1% accuracy. An additional practical advantage in favor of the series expansion that is not reflected in the results in Table 2.5 is that once P0(m,n) has been computed, only a little further computational effort is required to obtain P(m,n) using (2.5) for as many values of h as desired. With the recursion method, the entire computation must be repeated for each set of parmeters P, p, h. In summary, by using the series expansion in P, we can approximate P(m,n) in a consistent and efficient manner where accuracy may be traded off for computational effort if necessary. For most channel parameters, only a few terms in the expansion are required in order to obtain good accuracy; typically such computations are at least an order of magnitude faster than when using the recursion method to obtain the exact result. For large values of n and P, efficiency suffers due to the requirement for greater numerical accuracy but at least excessive amounts of memory are not required. 2. Approximation By Averaging n=1000 P=.0001, p=.l, h=.5 Method Time to compute all P(m,n) recursion 8 min. series expansion (C program) r=5 10 sec. r=10 15 sec. r=20 33 sec. series expansion (Maple program) Time to compute all P0(m,n) r=5 3 min. r=10 5 min. r=20 45 min. Table 2.5 Computational times of P(m,n) for the recursion series expansion methods. 2. Approximation By Averaging 29 2.4 Some Observations on Short BCH Codes For codes with k<20, it is feasible to compute their P u exhaustively by using the code generator matrices to generate all codewords and summing over their probabilities as given by (1.3). By comparing with the results from applying (2.1) or (2.2), we can study the validity of the approximation (2.1), and at the same time, the relative performance of a given code and its equivalent codes over the Gilbert channel. We will consider in particular the (primitive, cyclic) BCH codes [1-3] with n<31 and k<16 as they are quite common in usage and have well known weight distributions. In addition, BCH codes are not unique in the sense described below; this property will also be interesting when considering P u in the Gilbert channel. A BCH code of length 2 m - l corresponds to a certain primitive polynomial of degree m. A primitive polynomial p(x) of degree m is an irreducible polynomial (one not divisible by any polynomial of degree less than m but greater than 0) such that the smallest integer n for which it divides x n+l is n=2m-l. Since there may be more than one primitive polynomial of degree m, BCH codes with identical n and k are in general not identical. They do however have the same weight distributions and so on the BSC have the same error detection performance. Hence, it should be interesting to see whether a similar situation holds in the Gilbert channel. The problem is similar in nature to that of the relative performance of equivalent codes that we have already mentioned. However, in this case, there is one reduction that we can make immediately because of the time-reversal symmetry of the error bits produced in the Gilbert channel. 2. Approximation By Averaging 30 It turns out that primitive polynomials come in pairs. If p(x) is a primitive polynomial, then so too must be its reciprocal polynomial p*(x) = xmp(l/x) . (2.20) (It is possible that p*(x)=p(x).) It follows that the generator polynomials defining the BCH codes of given (n,k) come in such reciprocal pairs [2]. Later in Chapter 3, we show that reciprocal generators give rise to cyclic codes with the same P u on the Gilbert Channel. Hence, we need to be concerned with only half of the possible number of BCH codes. We consider 7 BCH codes of different (n,k). The channel parameters taken are p={0.01,0.1,0.3}, h={0.5,0.7,0.9} and P={10-6,10-5-5,...,10-25,10-2}. Results are summarized below. 1. Harnming (7,4) code A Hamming code is a special case of a BCH code with k=2m-m-l. For m=3, there is only the primitive polynomial x3+x+l and its reciprocal. By the discussion above, it suffices to study just one of them in the Gilbert channel. The code appears well behaved with P u monotonically increasing with P, decreasing with p and h. The actual P u is always less than the approximate Pu, which we henceforth denote by E[PJ, implying that the Hamming code is superior in error detection performance than the average of its equivalent codes. The percentage difference is only mildly sensitive to and decreases with P, but increases with p and decreases with h. For p=.3, h=.5, the difference is 18%. Generally, P u rises linearly with P for P £ 10'3. 2. Approximation By Averaging 31 2. Hamming (15,11) code For m=4, there is also just one primitive polynomial x4+x+l and its reciprocal. The qualitative behavior is the same as for the Hamming (7,4) code. The maximum difference between the exact and approximate P u reaches a maximum of 51%. 3. B C H (15,7) code Again, P u is monotonic with increasing P and decreasing p and h. Generally, the exact P u lies below EtPJ with a maximum difference of 140% when p=.3, h=.5. At h=.9, p=.01,.l the difference is actually negative (i.e. the approximate P u is larger than the exact one). However, the maximum difference is only -.5%. 4. B C H (15,5) code The qualitative properties are very similar to those of the Hamming codes. The maximum difference between the exact and approximate P u reaches 90% at P=.3, h=.5. 5. B C H (31,16) code For m=5, there are three primitive polynomials: (i) x5+x2+l, (ii) x5+x4+x3+x2+l, and (iii) x5+x4+x2+x+l and their respective reciprocals. We will call the codes corresponding to these polynomials by the same numbers. Qualitatively, the P u of the three different codes behave the same as that of the Hamming codes. The largest deviation between the exact Pu's and the approximate EtPJ occurs at p=.3, h=.5 and is about 950%. At h=.5, the differences between the three codes appear indistinguishable. The differences begin at emerge at low Pu, i.e., small P, large p, h, however; then code (i) has consistently lower P u than code (ii) which in turn has consistently lower P u than code (iii). At P=10"6, p=.3, h=.9, 2. Approximation By Averaging 32 Pu(code (i))=4.0 x IO"15, Pu(code(ii))=6.1 x IO"15, and Pu(code(iii))=9.6 x IO"15. The differences with the approximate P u are 176%, 47%, and 18% espectively. 6. B C H (31,1D code The P u of these three codes behave similarly to that of the B C H (15,7) code. Maximum difference between exact and approximate P u is 540% at p=.3, h=.5. At h=.9, p=.01,.l the difference is negative with a maximum difference of only -1% for code (i), -.1% for code (ii), and -.2% for code (iii). There is a slight difference in P u of 5% to 10% in regions of low P u. 7. B C H (31,6) code The three codes here have identical P u. See Chapter 4, Part 3. Qualitatively, they are similar to the Hamming codes. The maximum difference between the exact P u and E[PJ is 400% at p=.3, h=.5. We can make some general comments based on the observations above. All the BCH codes appear to have proper behaviour in that P u increases monotonically with P and decreases monotonically with p and h. This is probably not surprising as the Hamming codes, double error correcting BCH codes like the B C H (15,7) code, the BCH (15,5) and triple error correcting codes with odd m such as the B C H (31,16) code are known to be proper [4,11]. The percentage difference between exact and approximate P u is only mildly dependent on P and usually decreases with P. This is related to the linearity of P u with respect to P as the exact and approximate P u have almost parallel trajectories. The reason for this is that for small enough P, the pervasive PB factor present in all probabilities (as given by (1.3) and (1.4)) contains the dominant linear contribution to P u from P. The maximum difference between 2. Approximation By Averaging 33 exact and approximate P u always occurs at the largest p and the smallest h values we considered, namely at p=.3 and h=.5. The exact P u is practically upper bounded by EtPJ. When the exact P u is observed to be greater, the difference is slight and always at low P u. The differences between codes with the same (n,k) become most noticeable also at low P u . It should be mentioned that a similar study of some 35 codes with lengths less than 25 was originally reported by Elliott [14]. However, no names of codes or details were given other than the comment that in general P u and E[PJ agreed within an order of magnitude except for extreme cases. 3. Some Exact Results 3.1 No Coding For later use, we record here the probability of obtaining the zero error vector. Since the event not 0n is the union of the mutually exclusive events lf01,.,0^1 [14], n-l prob(0n) = 1 - E probC^i) i=0 n-l = 1 - I W(i) i=0 = i . p B ( i . h ) i ; « ^ i=0 pB(l-h) [ 1-Jn 1-Ln1 = i - J X ~ (j+p-Q)-Tg--(L+p-Q)-n; • (3.D If no coding is used, any nonzero error vector results in an undetected error so the probability of undetected error is P U = £ ^ M ( J + p . Q ) l ^ . ( L + p . Q ) k ^ ( 3 2 ) 3.2 Single Parity Check Encoding Let G(m,n) and B(m,n) be the probability of m errors in a block of length n given that the channel is initially in state G and state B respectively. The cor-34 3. Some Exact Results 35 responding unconditioned probability is just the counting distribution of Chap-ter 2 and can be decomposed as P(m,n) = pGG(m,n) + pBB(m,n) . (3.3) The conditioned probabilities may be found recursively from [14] G(0,1) = 1, G(1,1) = 0, B(0,l) = h , B(l,l) = h.' G(m,n) = QG(m,n-l) + PB(m,n-l) B(m,n) = hqB(m,n-l) + h'qB(m-l,n-l) + hpG(m,n-l) + h'pG(m-l,n-l) (3.4) where of course G(m,n) and B(m,n) are zero if m<0 or m>n. Observe that the last equation gives a relationship between the probabilities of odd and even numbers of errors. By exploiting this fact, we will derive an analytic expression for the probability of undetected error for a single parity check code. A codeword in a single parity check code of length n consists of n-l mes-sage bits and a parity check bit which is 1 or 0 depending on whether the mes-sage contains an odd or even number of 2's. Hence a single parity check code consists of all possible even weight vectors and no odd weight vector. From the standpoint of probability of undetected error, it does not matter where the parity bit is inserted (or whether the original message bits are rearranged). Ultimately, what we are interested in is P u = I P(m,n) . (3.5) meven m^O Note that what we considered above is the even parity version of the single parity check code. For the odd parity version of the code, the parity check bit is 3. Some Exact Results 36 instead the complement of that for the even parity case so the codewords are all the odd weight vectors. The odd parity version is not a linear block code since it does not contain a zero vector. However, since an odd/even weight vector added to an odd weight vector results in an even/odd weight vector, an undetected error results if and only if the error vector has even weight. That is, the odd parity and even parity codes have in fact the same performance as given by (3.5). While we may think of using (2.6) in (3.5), it is immediately seen that the resulting summations are not easily tractable. The idea is to use (3.3) and (3.4) to effectively divide the summation into two parts and then solve for each by more elementary means. To this end, define G„(n) = E G(m,n) modd Ge(n) = E G(m,n) m even B0(n) = E B(m,n) m odd Be(n) = E B(m,n) (3.6) m even where necessarily G0(n) + Ge(n) = 1 B0(n) + Be(n) = 1 . (3.7) Summing over the recursion relations (3.4) gives G0(n) = E QG(m,n-l)+PB(m,n-l) modd = QG 0(n-l) + PB0(n-l) 3. Some Exact Results 37 B0(n) = E hqB(m,n-l) + h'qB(m-l,n-l) + hpG(m,n-l) + h'pG(m-l,n-l) m odd = h qB 0(n-l) + h'qBe(n-l) + hpG0(n-l) + h'pGe(n-l) = (h-h')qB0(n-l) + (h-h')pG0(n-l) + h' . (3.8) The last equations do not quite allow us to solve for B0(n) and G0(n) because of the different arguments on the right hand side. It is convenient at this point to introduce the formalism of generating functions. For a random variable over the nonnegative integers X(n), its generat-ing function is given by the formal sum OO x(t) = E X(n)tn . (3.9) n=0 For sufficiently well behaved X(n), x(t) will be well defined for some range of t. The individual probability X(n) is then the coefficient of t n in a power series expansion of x(t). Applying the definition (3.9), if g0(t) and b0(t) are the generating func-tions for G0(n) and B0(n) respectively, then go(t) = go(t)-G o ( l)t-G o(0) = E G0(n)tn n=2 OO 00 = Q E G 0(n-l)tn + P E B0(n-l)tn n=2 n=2 = Qtg0(t) + Ptbo(t) 3. Some Exact Results 38 bo(t) = b o (t)-B o (l)t-B o (0) + h ,t oo = £ B0(n)tn + h*t n=2 = (h-h')q E B0(n-l)tn + (h-h')p E G0(n-l)tn + h' E tn n=2 n=2 n = l = (h-h')qtbo(t) + (h-h')ptgo(t) + • (3.10) Gathering similar terms gives ( l - Q t ) g o ( t ) - Ptbo(t) = 0 -(h-h')ptg0(t) + [l-(h-h')qt]b0(t) = ^ | . (3.11) The solution of the simultaneous equations is g 0 ( t ) = h'Pt2 i.[(h-h,)q+Q]t-(h-h,)(l-q-Q)t2  h'Pt2 "(l-tXl-JitXl-Lxt) b 0 ( t ) = ^ g „ ( t ) (3.12) where J i = J i ( Q , g , h ) = J(Q,q,h-h') = J(Q,q,2h-l) L i = Lx(Q,g,h) = L(Q,q,h-h') = L(Q,q,2h-l) . (3.13) Now if P0(n) is the unconditioned probability of obtaining an odd number of errors in a block of length n, its generating function is just Po(t) = PGgo(t) + pnbo(t) = p^(p-Q+Y)go(t) • (3.14) 3. Some Exact Results 39 Expanding g„(t) by partial fractions and the result by power series gives p0(t) = h'pBt2 (p-Q+|) f J i 2 1 W 1 * [(Jx-LxXJi-1) l-Jit + ax-JiKLx-l) 1-L; i t T ( J r l X L i - l ) 1-tJ = h'pBt+ ^ | { l + p ( ^ ^ [ ( J i + P - Q ) J i n ( L i - l H L 1 + p - Q ) L i « 1 ( J 1 - l ) ] t n n=2 (3.15) The last expression allows us to recover P„(n) as the coefficient of the tn term in the power series expansion, that is Po(l) = h'pB P 0 ( n ) = | j l + p ( 3 ^ [ ( ^ ; n>2 . (3.16) The probability of undetected error for a single parity check code is therefore P u = 1 - P0(n) - prob(On) = - \ {l+ p ( j P x ! L l ) [(Ji+p-Q)Jin(L1-l)-(L1+p-Q)L1n(J1-l)]} + ^ ( W + U ) ^ - ( L + P - Q ) ^ 1 . (3-17) As a check of this expression, we can consider for instance its large n limit. First, note that the quantity under the radical sign in the definition for J and L (1.5) can be written [(Q-hp)-h]2+4h(l-h)p>0 and a straight forward application of Descarte's Rule of Signs shows that for p>Q, 1>J>0, 0>L>-1 and for p<Q, 1 > J ,L > 0. That is, IJI, ILI, IJXI, I L x I < 1. Hence all the exponentiated terms in (3.17) vanish for large n. The remaining terms combine to give Pu=l/2, expectedly the same result as for the BSC. 3. Some Exact Results 40 Figures 3.1 to 3.12 are plots showing P u as a function of P and p for fixed values of h and n. We see graphically that P u is monotonically increasing with P and decreasing with p and h. Hence for the single parity check code in the Gilbert channel, as channel conditions worsen, the probability of undetected error increases. This type of behaviour corresponds to that of a proper code in the BSC, that is a code for which P u increases monotonically up to (l/2)n"k with the bit error rate [4]. The P u also increases quite sharply with n for all channel parameter values, similar to the BSC case, as the single parity check code is quite limited in its error detection capabilities since it detects no even number of errors. 3.3 Probability of a Vector and its Cyclic Shifts According to (1.3) and (1.4), the probability of an error pattern e is prob(e) = pBu(a) • wt(e>l MvCbi) i=l u(c) (3.18) where a is the number of 0's before the first 1, bi is the number of 0's between the ith and (i+l)th 1, and c is the number of 0's following the last 1 in e. Denote a right cyclic shift of s places by 7ts. A negative s will correspond to a left cyclic shift. Then the probability of e and its next c right cyclic shifts and a left cyclic shifts is £ prob(7tse) = p B wt(e)-l n v ( b i ) i=l bo £ u(s)u(b0-s) (3.19) s=0 Figure 3.1 P u for single parity check code, n=8, h=.5 Figure 3.2 P u for single parity check code, n=8, h=.7 Figure 3.3 P u for single parity check code, n=8, h=.9 3. Some Exact Results 44 Figure 3.4 P u for single parity check code, n=16, h=.5 Figure 3.5 P u for single parity check code, n=16, h=.7 3. Some Exact Results 46 Figure 3.6 P u for single parity check code, n=16, h=.9 Exact Results Figure 3.7 purorsingl 3. Some Exact Results 48 Figure 3.8 P u for single parity check code, n=32, h=.7 Figure 3.9 Pu for single parity check code, n=32, h=.9 3. Some Exact Results 50 Pu n=64, h=.5 Figure 3,0 P . i - P * * « * " " ^ Figure 3.11 P u for single parity check code, n=64, h=.7 3. Some Exact Results 52 Figure 3.12 P u for single parity check code, n=64, h=.9 3. Some Exact Results 53 where for convenience we have put b0=a+c. The sum d u2(d) = E u(s)u(d-s) (3.20) s=0 is easily handled using generating function methods. From [12], the generating function for u(-) is l +(p-Q)tj_J___L_j u w _ J-L J 1-Jt 1-Ltj Since u 2 is the self convolution of u, the generating function for u 2 is (3.21) U2(t) = U(t)2 l+2(p-Q)t+(p-Q)2t2 (J-L)2 J 2 2JL (1-Jt)2" J-L 1-Jt 1-Lt L 1 L 2 1 1- tJ + (1-Lt)2| • ( 3 , 2 2 ) By expanding U2(t) as a power series in t, it follows that u2(k) = {(k+l)(Jk+2-rLk+2)+2k(p-Q)(Jk+1+Lk+1)+(k-l)(p-Q)2(Jk+Lk)} 2 JL - TJT^i {(Jk+1+Lk+1)+2(p-Q)(Jk+Lk)+(p-Q)2(Jk-1+Lk-1)} . (3.23) Putting b=(b0 v..,bw t(e)-i), the total probability of e and all n of its cyclic shifts is therefore n-l wt(e>l X prob(7tse) = £ u2((7ttb)o)v((7rtb)1)...v((7ctb)wt(e)-i) • (3.24) s=0 t=0 The last equation is invariant with respect to e-*Jtse as it must be since it is the probability of a cyclic shift class of vectors of which e is but one member. As an example, consider the dual Hamming codes. (A dual code of a lin-ear block code C(n,k) is a linear block code C*(n,n-k) consisting of all vectors 3. Some Exact Results 54 which have vanishing inner product with the codewords of C(n,k) [1-3].) The dual Hamming (2m-l,m) code consists only of the zero vector and a single cyclic class of vectors of weight 2m'1 generated by the reciprocal of the parity polynomial h(x)=(xn+l)/g(x) where g(x) is the generating polynomial of the Hamming (2m-l,2m-m-l) code. For m=3 and g(x)=l+x+x3, h*(x)=l+x2+x3+x4 so the probability of undetected error is the probability of the cyclic class generated by the codeword pattern (101110 0), that is P u= u2(2)v(l)v(0)v(0) + u2(l)v(0)v(0)v(2) + u2(0)v(0)v(2)v(l) + u2(0)v(2)v(l)v(0) = v(0){u2(2)v(l)v(0)+u2(l)v(0)v(2)+2u2(0)v(2)v(l)} . (3.25) Of course, in this simple case, one could just as well have summed the 7 codeword patterns explicitly and obtained the same result. As another example, consider the BCH (31,6) code encountered in Chapter 2. The distinct generator polynomials (not counting reciprocals) are known to be (i) . l+x+x2+x5+x9+xu+x13+x14+x15+x16+x18+x19+x21+x24+x25 , (ii) . l+x4+x5+x6+x8+x10+x13+x15+x16+x17+x18+x21+x22+x24+x25 , (iii) . l+x+x2+x3+x6+x7+x9+xn+x14+x15+x16+x17+x21+x22+x25 . A BCH code always contains the all i's vector since the generator polynomial never contains the (x+1) factor. The weights of the codewords corresponding to the generator polynomials are all 15, so the B C H (31,6) code must just consist of two cyclic classes of 31 codewords each, with respective weights 15 and 16 in addition to the all 0's and all l's codewords. Writing out the components of the b vector as defined after (3.23) corresponding to the above code polynomials, we get 3. Some Exact Results 55 (i) . (5,0,0,2,3,1,1,0,0,0,1,0,1,2,0) (ii) . (5,3,0,0,1,1,2,1,0,0,0,2,0,1,0) (iii) . (5,1,0,0,1,0,1,1,2,0,0,0,3,0,2) which clearly give rise to the same probabilities for the respective cyclic classes according to (3.24). The situation is the same with the weight 16 codewords which are just the complements of the weight 15 codewords. Hence without doing much work, we have shown that all the B C H (31,6) codes have the same P u in the Gilbert channel. Similarly, an explicit expression for P u can be written using the rule (3.23) for other simple codes which have but a few readily identifiable cyclic classes of codewords such as the Hamming (7,4) code, etc. However, that is not the point of the exercise. Rather, what we want to point out to conclude this section is the following. In a general channel, all 2k codewords can conceivably have different probabilities. If there is no way to analytically relate the codewords, the exact determination of P u (for large k) is surely an intractable problem. For the BSC, the codewords are divided into weight classes, each characterized by a certain probability. For the Gilbert channel, we know this tactic will not be immediately applicable. However, due to the underlying symmetry of the Gilbert channel, we see from the above that n vectors which are cyclic permutations of each other can be grouped together and have a total probability dependent on only a single pattern. The cyclic codes appear to fit naturally into this scheme. As there are no other obvious channel and code symmetries that can be identified this way, the cyclic codes probably are the best candidates for exact analysis. The problems are to first identify or enumerate the various cyclic classes and then 3. Some Exact Results 5 6 provide summation rules between them. Of the latter, we would expect to apply techniques similar to that which has been given here for summing over collections of u2(-) and v(-)-3. 4 Reciprocal Generator Polynomials Here we derive an interesting fact regarding the cyclic code generated by the reciprocal of the generator polynomial of a given cyclic code. If g(x) is a generator polynomial of a cyclic code, it necessarily divides x n+l. It follows that the reciprocal polynomial g*(x) must also divide x n+l and therefore it also generates a cyclic code of the same length [1-3]. Denote the degree of g(x) (and g*(x)) by q. A code polynomial in the cyclic code generated by g(x) may be written c(x)=g(x)s(x) where s(x) is of degree at most n-q-1. Write s(x)=(xa+...+l)xb (where a+b<n-q-l) so that c(x)=g(x)(xa+...+l)xb. Now in the code generated by g*(x), there is a codeword given by d(x)=g*(x)(xa+...+l)*xb. Suppose first for simplicity that b=0, then d(x)=c*(x). In coordinate form, the first q+a+1 bits of the vectors c and d will be mirror image patterns of each other with the 1st and (q+a+l)th bits l's. The remaining n-q-a-1 bits are 0's. From (1.3) and (1.4), it is clear that mirror image patterns have the same probability. The 0's following the last 1 just contribute a factor u(n-q-a-l) to the net probability of the codeword in both cases. Hence prob(c)=prob(d). Now if b ^  0, then the vectors c and d will both have b 0's at one end (contributing a factor w(b) to the net probability) and n-q-a-b-1 0's at the other end (contributing a factor u(n-q-a-b-l) to the net probability) with mirror image patterns of length 3. Some Exact Results 57 q+a+1 delimited by i's in the middle. Therefore, again prob(c)=prob(d). Clearly, the c(x) and d(x) above define a one to one correspondence between codewords in the cyclic codes generated by g(x) and g*(x) respectively. Hence cyclic codes generated by reciprocal generator polynomial pairs have identical P u on the Gilbert channel. 4. Monte Carlo Simulation 4.1 The Monte Carlo Method When analytical methods are unavailable, a common technique used to study the performance of communication systems is computer simulation. Usually, the formalism is phrased in terms of that for bit-error rate estimation [19] so first we will rewrite the terminology so that it is appropriate for our case. The probability of undetected error of a code with block length n can be written P u = I H(v)prob(v) (4.1) v € V n where j 1 v e C , v * 0 . . . . H ( v ) = 0 otherwise ( 4 - 2 ) is the error detector. In other words, P u is the expectation value of the error detector. A natural estimator of the expectation is just the sample mean P\, = ^ l H ( v i ) (4.3) i=l where v; is the ith of a total of N blocks sent through the channel. Hence a suit-able basis for estimating the error rate is by the observation of errors; this defines the Monte Carlo method. As N-*oo, the estimate $ u converges to P u . For finite N, the reliability of the estimator is quantified in terms of a confidence level a 5 8 4. Monte Carlo Simulation 59 and its associated confidence limits, y+(F\i) and y.(P\0 defined through the rela-tion prob(y+<Pu<y.)= 1-a . (4.4) P\j is binomially distributed. By applying the normal approximation to the binomial distribution, it can be shown that 1+ d „ 2 2NPU . 1±A i+-4NP\, da 2 (4.5) where d a satisfies da dt e-t2/2 = 1- a (4.6) -da For our purposes, we always take a=.01, then da=2.575829, and the corre-sponding confidence limits define a 99% confidence interval. 4.2 Simulation Study of CRC-16 Codes A CRC code is defined by a qth degree generator polynomial g(x)=gqxq+gq.iXq"1+...+go. Denote the polynomial representing the message by s(x)=Sk-ixk"1+Sk.2Xk"2+...-rSo. A codeword is formed by appending to the message bits certain parity bits corresponding to the so called CRC polynomial r(x)=remainder[xqs(x)/g(x)]. The polynomial representation of the codeword is then c(x)=s(x)xq+r(x)=Sk.ixq+k-1+sk.2x<i+k-2+.. .+s0xq+rq.1xq-1+rq.2xq-2+.. .+r0. CRC codes are often used in practical error control as the encoding (modulo-2 4. Monte Carlo Simulation 60 polynomial division) is particularly simple to implement by feedback shift regis-ter circuitry. If the generator polynomial is formed from the product of (x+1) and a primitive polynomial of degree q-1 (and the block length n is less than 21-1), it can be shown that the resulting CRC code detects all double errors, all odd numbers of errors, and all bursts (including end-around) of length q or less [21]. There are two commonly used standard CRC codes with q=16, both have the required form discussed above. Their generator polynomials are X i 6 + X i 5 + X 2 + L for CRC-ANSI , x16+x12+x5+l, forCRC-CCITT . The performance of these CRC codes have been studied in the BSC by Witzke and Leung [9]. Both codes are not proper for all values of k. The probability of undetected error as a function of the error rate e reaches a maximum before decreasing to the asymptotic value of (1/2)"* at e=l/2. The CRC-CCITT code performs better (in terms of probability of undetected error) than the CRC-ANSI code for all values of e. In view of the connection with burst detection and the previous study, we have chosen the CRC-ANSI and CRC-CCITT codes as the subject for Monte Carlo simulation study. While we are mainly interested in questions of practical concern such as whether the superiority of CRC-CCITT over CRC-ANSI persists over a Gilbert channel and how the results compare with those of the BSC, the study will also give some general idea of the size of problem that can be tackled in a practical implementation of the Monte Carlo technique on the Gilbert channel. We take k=25,50. The simulation program is written in the C language, using the library function rand() for random number generation. Random bit-vectors 4. Monte Carlo Simulation 61 v=(i> 0 ,Ui, . . . , i>n-i) are generated according to the distribution implied by Figure 1.1. That is, before the first bit is generated, the state of the channel is set ran-domly to the G or B state according to the ratio PG/PB- At each successive bit, if the channel is in the G state, a 0 is generated and a transition to the B state occurs randomly with probability P. If the channel is in the B state, a 0 or 1 bit is generated according to the ratio h/h' and a transition to the G state occurs ran-domly with probability p. The procedure ends after n bits corresponding to an entire vector are generated. The probability of any vector so produced will be in accordance with (1.3) and (1.4). The vector is checked to determine whether it is a codeword using shift register arithmetic. If so, the error counter is incremented as indicated by (4.2). This procedure of vector generation and checking is repeated N times. The estimated probability of error and the 99% confidence levels are then computed using (4.3) and (4.5). Running the simulation on the SPARC® station 1, we have found that for N=107 and k=25, a single run (corresponding to one set of parameters P,p,h) requires on the average of approximately 2.5 hours. For k=50, the average time required is approximately 4.0 hours. It is clear from the nature of this simulation that the computational time increases linearly with n. The parameter values which we have considered are P={10-4,10-3-5,10-3,10-2-5,10-2>10-L5,10-1,10-5,l}, p={0.01,0.1,0.3}, h={0.5,0.7,0.9}; therefore, 81 "data points" $ n in all. The results of our simulation are plotted in Figures 4.1 to 4.9. For clarity, we do not draw lines through the simulation data points. The error bars around these points bound the 99% confidence intervals. The lines on the plots are the P u of the CRC codes for the BSC at the corresponding 6eff. We have chosen to display the plots with P along the x-axis and fixed values of p and h. This is Monte Carlo Simulation Figure 4.1 S imula t ion of CRC-16 codes, p = .01, h=.5 Monte Carlo Simulation Figure 4.2 Simulation of CRC-16 codes, p = .01, h = .7 Monte Carlo Simulation Figure 4.3 Simulation of CRC-16 codes, p = .01, h=.9 4. Monte Carlo Simulation 65 Figure 4.4 Simulation of CRC-16 codes, p = . l , h = .5 Monte Carlo Simulation Figure 4.5 Simulation of CRC-16 codes, p=.l, h=.7 Monte Carlo Simulation Figure 4.6 Simulation of CRC-16 codes, p=.l, h=.9 4. Monte Carlo Simulation 6 8 Figure 4.7 Simulation of CRC-16 codes, p=.3, h=.5 4. Monte Carlo Simulation 69 Figure 4.8 Simulation of CRC-16 codes, p = .3, h = .7 4. Monte Carlo Simulation 70 Figure 4.9 Simulation of CRC-16 codes, p = .3, h=.9 4. Monte Carlo Simulation 71 reasonable since P is the quantity most like the bit error rate in the BSC, that is, it is usually the small quantity. It is evident from the plots (and even before hand) that the minimum P u we can hope to simulate reasonably with our sample size is ^10/N=10"6. This corresponds to the range 1 0 " 4 £ P £ 1 0 " 2 depending on the values of p and h. In spite of the limitations, a number of features in the plots are evident. Firstly, CRC-CCITT appears to maintain its superiority over CRC-ANSI over all parameter values. The difference between them is small for p=0.01, h=.5, but that difference increases to more than an order of magnitude difference in P u with increasing p and h. This is reasonable since for small p, the tendency is for the channel to stay in the B state if it gets there, and if h=.5, which is the random vector limit where all bit strings are equally likely, the two codes have the same probability of undetected error 1/216. Secondly, the qualitative behaviour of both codes for k=25 is similar to k=50. Quantatively, the k=50 codes reach the plateau at 1/216 for smaller values of P which is just as expected. Thirdly, the general behaviour of both CRC-ANSI and CRC-CCITT are most like their BSC effective counterparts when p=.3 where a maximum at approximately 10"4 is also evident. This is expected since for large p, the mean time that the system stays in the B state is short so the correlation between errors is small. Finally, note that there are ranges of parameter values where the Gilbert channel P u are quite different (order of magnitude or more) from the BSC effective values. For instance, when p=.01 and h=.5 and 10"4<P<10"2, P u is consistently overestimated by the BSC effective values; and perhaps more significantly, the opposite is true for the region p=.01, h=.9, and P< IO 3 . These observations point out that while some of the properties of CRC codes on the BSC may persist over the Gilbert channel, there are certain parameter value 4. Monte Carlo Simulation 72 ranges where interpolation is not possible, presumably because of fundamental differences between the channels. 5. Conclusions We have examined in this thesis the general problem of computing the probability of undetected error for binary linear block codes on the Gilbert channel. Little on the subject is known to date so we have tried different approaches. The accuracy of the code evaluation technique using the P(m,n) distribu-tion is difficult to predict as we have seen even in limited observations on short BCH codes. Certainly, precise numerical results are not possible. The fact that equivalent codes, and in the case of cyclic codes, codes which have different generator polynomials but same weight distribution, give rise to different P u is an additional albeit interesting complication. Perhaps the role of this approxi-mation technique should be that of preliminary or cursory investigation. In that case, the efficient computation of P(m,n) which we have developed can be uti-lized most advantageously. Even though exact analysis in the Gilbert channel is quite difficult in gen-eral, it is probably one of the few channel models other than the BSC where there is a realistic chance to obtain analytical results. We have made some exploratory attempts in this direction and have obtained a closed form expression for the single parity check code. The identification of the channel's underlying symme-tries also allowed us to derive two results on cyclic codes which can be used to facilitate the computation of P u . The results of our simulation of CRC codes give us impetus for further study. For example, it is clear that the P u for a code on a Gilbert Channel can be quite different from that on a BSC. In addition, we have observed some proper-73 5. Conclusions 74 ties of CRC codes which are themselves of practical and theoretical interest, esp-ecially in comparison with the corresponding BSC results. On the subject of the simulation itself, the standard Monte Carlo technique which we used is the most rudimentary, and therefore not very efficient computationally. Other techniques for simulation can make use of partial analytical information and involve more specific details than standard Monte Carlo [19]. Perhaps one of these alternate methods can be adapted for use on the Gilbert channel. Finally, though we have not considered it here, it should be worthwhile to examine the related problem of error correction [2]. On the Gilbert channel, neighbouring codewords which are the same Hamming distance away from a given vector can have different a posteriori probabilities. Therefore, it would be interesting to examine for instance, the effectivness of the nearest neighbour decoding rule (which is optimal on the BSC) or the modification of this rule to improve performance. References [I] F. J. MacWilliams and N. J. A. Sloane, The Theory of Error-Correcting Codes, North-Holland, Amsterdam, 1977 [2] S. Lin and D. J. Costello, Error Control Coding: Fundamentals and Applications, Prentice-Hall, New Jersey, 1983 [3] V. Pless, Introduction to the Theory of Error-Correcting Codes, Wiley, New York, 1989 [4] C. Leung, E. R. Barnes, and D. U. Friedman, "On some properties of the undetected error probability of linear codes," I.E.E.E. Trans. Info. Theory, vol. 22, pp. 235-237, Mar. 1979 [5] J. K. Wolf, A. M. Michelson, and A. H. Levesque, "On the probability of undetected error probability of linear codes," I.E.E.E. Trans. Commun., vol. 30, pp. 317-324, Feb. 1982 [6] T. Kasami, T. Klove, and S. Lin, "Linear block codes for error detection," I.E.E.E. Trans. Info. Theory, vol. 29, pp. 131-136, Jan. 1983 [7] C. Leung, "Evaluation of the undetected error probability of single-parity check product codes," IEEE Trans. Commun., vol. 31, pp. 250-253, Feb. 1983 [8] T. Fujiwara, T. Kasami, A. Kitai, and S. Lin, "On the undetected error probability for shortened Hamming codes," IEEE Trans. Commun., vol. 33, pp. 570-574, June 1985 [9] K. A. Witzke and C. Leung, "A comparison of some error detecting CRC code standards," I.E.E.E. Trans. Commun., vol. 33, pp. 996-998, Sept. 1985 [10] K. A. Witzke, "Examination of the undetected error probability of linear block codes," MASc Thesis, Univ. of British Columbia, Aug. 1984 [II] C. T. Ong, "On the undetected error probability of linear codes," MASc Thesis, Univ. of British Columbia, Mar. 1990 [12] E. N. Gilbert, "Capacity of a burst-noise channel," Bell Sys. Tech. J., vol. 39, pp. 1253-1265, Sept. 1960 75 References 76 [13] A. A. Alexander, R. M. Gryb, and D. W. Nast, "Capability of the telephone network for data transmission," Bell Sys. Tech. J., vol. 39, pp. 431-476, May 1960 [14] E. O. Elliott, "Estimates of error rates for codes on burst-noise channels," Bell Sys. Tech. J., vol. 42, pp. 1977-1997, Sept. 1963 [15] W. Feller, An Introduction to Probability Theory and Its Applications, vol. 1, Wiley, New York, 1968 [16] E. O. Elliott, "A Model of the switched telephone network for data communications," Bell Sys. Tech. }., vol. 44, pp. 89-109, Jan. 1965 [17] V. Cuperman, "An upper bound for the error probability on the Gilbert channel," I.E.E.E. Trans. Commun. Tech., vol. 17, pp. 532-535, Oct. 1969 [18] J. B. Cain and R. S. Simpson, "The distribution of burst lengths on a Gilbert channel," I.E.E.E. Trans. Info. Theory, vol. 15, pp. 624-627, Sept. 1969 [19] M . C. Jeruchim, "Techniques for estimating the bit error rate in the simulation of digital communication systems," IEEE J. Select. Areas Commun., vol. 2, pp. 153-170, Jan. 1984 [20] L. N. Kanal and A. R. K. Sastry, "Models for channels with memory and their applications to error control," Proc. IEEE, vol. 66, pp. 724-744, July, 1978 [21] D. Bertsekas and R. Gallager, Data Networks, Prentice-Hall, New Jersey, 1987 Appendix A The derivation of (2.6) is presented here for more than just the sake of completeness as there is apparent confusion regarding the range of m for which the equation is valid in [17]. The generating function (see (3.9)) of the P0(m,n) distribution relative to n for m> 1 is given by [16,17] £ P0(m,n) z n = p^T z m [q-^q-P)]-1 (l-zQ)--"1 [l-z(q-P)]2 . (A.l) Expanding the binomials and then gathering factors of z gives p oo m-1 2 I P0(m,n) z n = ^ — I I I c ^ c ^ c 2 qm* 1Q i(P-q) l + t z i + J + t + m . (A.2) n=m r + P i=0 j=0 t=0 The coefficient of z n on the right hand side is just P0(m,n), hence P0(m,n) = p ^ I c ^ ^ c ^ c 2 q-HQ^-HP-q)*^ . (A3) 0<t<2 j+t<n-m There are three cases to consider depending on the value of n-m. For n-m=0, (A.3) reduces to P0(n,n) = p ^ q n - 1 (A.4) For n-m=l, (A.3) reads P0(m,n) = p^" I C ^ C ^ C ? q m-j- lQn.m.j - t ( p. q ) i + t ( A > g ) ^ O^j^m-l Ost<2 j+t<l 77 Appendix A 78 but n-m-j-t=l-j-t<0 if j+t > 1 while n-j-t=m+l-j-t>0 so the cn.n^;|.t factor in the summand is zero for j+t > 1 hence (A. 5) is the same as m-l 2 P0(m,n) = p — E E CJV+KC? qm-J-lQn-m-j-t (p.q ) j +t ( A g ) r + P j=0 t=0 For n-m>2, (A.3) becomes p 2 rnin(m-l,n-m-t) P0(m,n) = p T T E E C ^ 1 ^ . t C ? qm-J-lQn-m-j-t ( p_ q ) i +t ( A ? ) r + P t=0 j=0 and there are three subcases to consider. If m-l<n-m-2, the upper limit of the summation over j is m-l so (A.7) becomes identical to the expression in (A.6). If m-l >n-m, however p 2 n-m-t P0(m,n) = E Eq c ^ X ^ t C 2 q^MQ^J-KP-q)"*1 p 2 I m-t m-l I = pTT E E - E C ^ X ^ j - t C 2 q m-J-lQn-m-j-t ( p . q ) j + t ( A g ) ^ v t=0 L j=0 j=n-m-t+lJ The second sum inside the braces (which should be ignored for m-l=n-m, t=0) is zero since n-j-t-m<n-(n-m-t+l)-t-m=-l and n-j-t>n-(m-l)-t=(n-m)+l-t>3-t> 1 means that c^j.,. is zero in the summand. Hence (A.8) is also the same as (A.6). Finally, if m-l=n-m-l, from (A.7) p 2 min(m-l,m-t) >o(m,n) = p^ p" EQ EO tf^C^C? qm-J-lQn-m-j-t (p.q ) j +t ' c n j " l c n n # t c 2 q^J^Q^-KP-q) 1 p 1 m-l 2 m-2 ^ I E + E E L T ^ 11=0 j=0 t=2 j=0 P 2 m-l = p— E E tf^c^c? q-^Q—J-t(P-q)i T^ 1.1=0 j=0 +t r.m-1 -1 (A.9) Appendix A 79 so that (A.6) is again obtained since c^^O. Consider now m=0. Po(0,n) is just the probability of being in state G for all n transmission intervals. For this to occur, the system must be initially in state G (with probability pG) and remain there (with probability Q) for the next n-l intervals hence Summarizing for the various cases considered above, we have the result given in (2.6), namely Po(0,n) = R (A. 10) P+p R Q ; m=0 P+p P m-l 2 P0(m,n) = < P+p I I c^-1cnPnJ;jitc2qm-J-1Qn-m-J-t(P-q)i+t ; l < m < n - l j=0 t=0 (A.11) Appendix B Here we evaluate the sums required for the computation of P0(m,n;r) in Chapter 2. We assume in the below that a>b+l and b> 1 so the following are certainly well defined: b-l U(a,b,c,d,e) = (-l)i(a-i-2)! j=0 (a-b-c-j)!(j-d)!(b-e-l-j)! V( W( b-l a,b,c,d)= Y I (-l)Ka-j-2)! j = o b-l a-b-c-j)!(j-d)!(b-j-D! (-l)i(a-i-2)! a,b,cj- ^ ( a . b . c . j ) ! j ! ( b . j . 1 ) ! • j=0 (B.l) (B.2) (B.3) By shifting the range of the summation variable, we can obtain relationships between the above quantities. Thus U(a,b,c,d,e) b-l (-DKa-j-2)! j = o (a-b-c-j)!(j-d)![b-l-(j+e)]! b+e-l so that U(a,b,c,d,e) = (-l)i[(a+e)-,i-2]! [(a+e)-b-c-j]![j-(d+e)]!(b-l-j)! j=e {b-l e-1 I - I j=0 j=0 (-l)i[(a+e),i-2] [(a+e)-b-c-j]!|j-(d+e)]!(b-l-j)! -l)eU(a+e,b,c,d+e,0) if d>0 -l) e U(a+e,b,c,d+e,0) ; d>0, e>0 -l) e V(a+e,b,c,d+e). if e > l (B.4) (B.5) 80 Appendix B 81 Also vr K . i w itf ^  (-l)»[(a-dH-2]! vta,b,c,al-t-i; ^ [(a-d)-(b-d)-(c+d)-j]!j![(b-d)-j-l]! j=-d — < 0 ; d>0,b-d<0 ("1)d2^ [(a-d)-(b-d)-(c+d)-j]!j![(b-d)-j-l]! ; d ^ 0 ' * 5 ' ^ 1 j=0 (-l)i[(a-d)-i-2]! [0 ; d>0,b-d<0 [(-l)dW(a-d,b-d,c+d) ; d>0,b-d>l (B.6) From standard references, for instance on pp. 62 of [15], W(a,b,c)=7^f](a-b-j)(b-l+j) ; a>2,b>l ,c>2 (B.7) We can apply (B.7) to the last line in (B.6) if a-d>2, b-d> 1, c+d>2. The first condition is redundant since we always take a>b+l. Hence for b-d> 1, c+d>2, d>0, i c+d-l V(a,b,c,d) = (-l)d (c+d-1), n (a-b-j)(b-d-l+j) • (B.8) However, note that in fact (B.8) holds even for b-d<0 if c > l since then b-d-l+l=b-d<0 and b-d-l+(c+d-l)=b+c-2>0 so one of the factors (b-d-l+j) must be zero. That is, we have I c+d-l V(a,b,c,d)= (-l)d (c+d.1), fi (a-b-j)(b-d-l+j) ; a>b+l, b > l , c> l , d>0, c+d>2 = (-l)d(c+d-l)!cc^.1 1ccb+1:1 . (B.9) Appendix B 82 We require two further identities regarding V(a,b,c,d). Using the explicit form on the right hand side of (B.9), it is simple to verify that V(a+l,b,c,d) = a . b a c " . d + l V ( a > b > c > d ) ( B 1 0 ) and \T, v J i \ -(a-b-c-d)Qb-d-l) V(a,b,c,d+U = ^ V(a,b,c,d) . (B.ll) Appendix C The following pages contain the source code listings for the C and Maple® programs used in the numerical computations and simulation. A brief descrip-tion of the purpose of each program follows. The parameters to be entered at runtime and any include files required are specified by comments at the beginning of each listing. 1. pOmn_c.c - C program that computes P0(m,n) from the closed form summation expression (2.6). 2. pmn_r.c - C program that computes P(m,n) from the recursion relations (2.3) and (2.4). 3. pmn_s.c - C program that computes P0(r)(m,n) and P^dn.n) from (2.5), (2.17), and (2.18). 4. pu.c - C program that computes Pu(C(n,k)) by exhaustive computation using explicit codeword generation and (1.3), E[PJ using (2.2) and the series expansion method of computing P(m,n), and their difference in percentage. 5. mc_crc.c - C program that generates a Monte Carlo simulation for CRC codes according to the procedure in Chapter 4. 6. pmn_c.map - Maple program that computes P0(m,n) from the closed form summation expression (2.6) and then P(m,n) from (2.5). 7. pmn_s.map - Maple program that computes P0(r)(m,n) and P^Xnijn) from (2.5), (2.17), and (2.18). 83 pOmn_c.c / * p0mn_c.c * / /* Computes each term of P0(m,n) and the total using a closed form summation expression given by Cuperman. */ / * parameters prompted for at runtime: n = block length m = number of errors in block P = G to B transition probability p = B to G transition probability h = B state bit error rate */ •include <stdio.h> Jinclude <math.h> / * factorial function */ double fact(n) int n; < int i ; double x=1.0; for (i=l; i<=n; i++) x = x * i ; return x; ) / * inverse factorial function */ double invfact(n) int n; ( i f (n<0) return 0.0; else return 1.0/fact(n); ) / * function A(j) = (n-m-j)p+mQ */ double a(m,n,pg,pb,j) int m, n, j ; double pg, pb; ( return (n-m-j)*pb + m*(l-pg); ) /* jth term in summation expression for P0(m,n) * / double pO(m,n,pg,pb,j) int m, n, j ; double pg, pb; ( return (pg/ (pg+pb)) *fact (n- j-2) *invfact (j) *invf act (m-j-1) *invfact (n-m-j) *pow(l-pb,(double)(m-j-1))*pow(pg+pb-1,(double)(j)) *pow(l-pg,(double)(n-m-j-2)) * (a(m, n,pg,pb, j) *a (m,n,pg,pb, j) - pb*a (m,n,pg,pb, j) - m*(1-pg-pb)*(1-pg)) / m; ) void main() { int m, n, j ; double pg, pb, x, y; printf("n = "); scanf("%d", in); printf Cm = "); scanf("*d", tm) ; printf("P = "); scanf("%lf", Spg); printf("p = "); scanf("%lf", Spb); pr int fP j p0(m,n, j)\n"); printf (" \n"); x = 0.0; for (j=0; j<=m-l; j++) ( y = p0(m,n,pg,pb, j) ; x = x + y; printf("%5d %13e\n", j , y); ) printf("total %13e\n", x); return; pmnr.c /* pmn_r.c */ /* Conputes P(m,n) for a l l m exactly using recurrence relations. */ •include <stdio.h> •include <math.h> •define N 1000 /* Define the value of N to be the block length. */ void main() ( i n t i , m, n; double x; double pg, pb, h; double u[N+l), v[N+l), p[N+l); double r[N+l][N+ll; p r i n t f ( " P = " ) ; scanf("»lf", Spg); p r i n t f ( " p = " ) ; s c a n f ( " % l f " , spb); p r i n t f ( " h = " ) ; s c a n f ( " % l f " , Sh); u[0] = 1.0; u[l] = h + (l-h)*pb; for (i=2; i<=N; i++) { u [ i ] = <l-pg+h*(l-pb))*u[i-l] + h*(pg+pb-1)*u[i-2]; ) v[0] = ( l - h ) * ( l - p b ) ; v [ l ] = (l-h)*(pb*pg+h*(l-pb)*(l-pb)); for (i=2; i<=N; i++) ( v [ i ] = (1-pg+h*(1-pb))*v[i-l] + h*(pg+pb-1)*v[i-2]; ) for (n=l; n<=N; n++) ( r [ l ] [n) = u[n - U ; for (m=2; nK=n; m++) ( x = 0.0; for (i=l; i<=n-nrt-l; i++) x = x + v [ i - l ] * r [ m - l ] [ n - i ) ; r[m)[n] = x; ) ) for (m=l; nK=N; m++) ( x = 0.0; for (i=l; i<=N-m+l; i++) x = x + (pg/(pg+pb))* (l-h)*u[i-l]*r[m)[N-i+ 1 ] ; ptm) = x; ) %d\n' printf("N p r i n t f ( " "> p r i n t f ( " for ( n F = l ; m<=N, pri n t f ( " % 4 . d N); P(m,N)\n"); \n"); m++) %e\n", m, p[m]); return; (D Q. O CO Cn pmn_s.c / * pmn_s.c * / / * C a l c u l a t e s PO(m, n) approximately by series expansion i n P and then the P(m,n) d i s t r i b u t i o n for a l l m i n the G i l b e r t channel. / * parameters prompted for at runtime: P = G to B t r a n s i t i o n p r o b a b i l i t y p = B to G t r a n s i t i o n p r o b a b i l i t y h = B s tate b i t e r r o r rate n = b l o c k length r = number of terms taken i n expansion * / •include <stdio.h> Iinclude <math.h> •define N 4095 / * Set N to at least n * / •define D 4.0 / * Set D to f i x overflow p o s s i b i l i t y i n p(m,n,h,vector) / * (-1)A1 for integer 1. * / int parity(1) int 1; ( if (<l /2)*2-l=0) return 1; else return - 1 ; ) /* Binomial coefficients. */ double binom(m, n) int m, n; { int i ; double x=1.0; i f (n>=0 SS m>=0) I i f (n<m) return 0.0; else ( for (i=0; i<=m-l; i++) x = x * (double)(n-i)/(double)(m-i); return x; ) ) i f (n>=0 SS m<0) return 0.0; i f (n<0 ii m>=0) ( for (i=0; i<=m-l; i++) x = x * (double)(m-n-l-i)/(double)(m-i); return parity(m)*x; ) i f (n<0 SS m<0) ( if (n<m) return 0.0; e l s e < for (i=0; i<=n-m-l; i++) x = x * (double) ( -m- l - i ) / (double) (n-m-i ) ; return parity(n-m)*x; ) ) ) / * The r t h term i n the expansion of P0(m,n)/pG i n P. * / double pOr(m, n ,pg,pb,r ) i n t m, n, r ; double pg, pb; ( i n t b , c; double x=0.0; i f (r=0) { i f (m=0) return 1.0; e l s e return 0.0; ) else i f (r==l) ( i f <m=0) return (1-n) *pg; i f (m>0 SS nxn) return 2*pg*pow(l-pb,(double)(m-1))*(1+((n-m-1)/2.0)*pb); i f (III—n) return p g * p o » ( l - p b , ( d o u b l e ) ( n - l ) ) / p b ; ) else ( i f (m=0) ( return parity(r)*pow(pg, (double) ( r ) ) *binom(r ,n- l ) ; ) i f (m>0 16 m<n) ( for (b=0; tK=r; b++) x = x + (n-m-b+l)*(n-m-b) *binom(b,r)*binom(r-l ,m+b-l)*pow(pb-l , (double)(b)) ; return pow(pg, (double)(r))*pow(1-pb,(double)(m-r)) *binom(r- l , n-m)*x/(r*(n-m)); ) i f (III—n) return 0.0; ) ) / * P0(m,n) accurate to order r i n P. * / double p0(m,n,pg,pb,r) i n t m, n; double pg, pb; ( i n t 1; double x=0.0; for (1=0; K = r ; 1++) x = x + pOr(m,n,pg,pb,l) ; return (pb/(pg+pb))*x; ) I* Sum of (l-h)*m*h*(i-m)*C(m,i)*vector[il over i=m..n. */ double p(m,n,h,vector) i n t m, n; double h; double vector[ ]; < i n t i ; double x=0.0, y=1.0; for (i=m; i<=n; i++) ( x = x + y * v e c t o r [ i ] ; y = (y*h*(i+l)) / (i+l-m) ; ) return pow( (pow (x,l/D) *pow(l-h,m/D)) ,D); v o i d main() ( i n t m, n, r; double pg, pb, h; double pO_o[N+l]; p r i n t f ( " P = " ) ; s c a n f ( " % l f " , Spg); p r i n t f ( " p = " ) ; s c a n f ( " * l f " , *pb); p r i n t f ( " h = " ) ; s c a n f ( " * l f " , i h ) ; p r i n t f ( " n = " ) ; scanf("%d", Sn); p r i n t f ( " r = " ) ; scanf("%d", t r ) ; p r i n t f (" m P0(m, n)\n"); p r i n t f (" \n"); for (m=0; m<=n; m++) pO_o[m) = pO(m,n,pg,pb,r) ; for (m=0; m<=n; nrt-+) ( p r i n t f ( " % 5 d m); printf ( " % . 4 e \ n " , pO_o[m]); ) p r i n t f ( " \ n " ) ; p r i n t f (" m P(m,n)\n"); p r i n t f (" \n"); for (m=0; m<=n; n»++) ( p r i n t f ( " % 5 d ", m); printf("%.4e\n", p(m,n,h,p0_o)) ; ) return; ) ' Q. O CO pu.c /* pu.c */ /* Calculates occurrence probability of a l l nonzero codewords of a code by exhaustive computation, the average of a l l equivalents to that code using the P(m,n) (approximate) di s t r i b u t i o n , and the i r percentage difference i n the G i l b e r t channel. */ /* include f i l e required: "codegen" - see below for format /* parameters prompted for at runtime: P = G to B t r a n s i t i o n probability p = B to G t r a n s i t i o n probability h = B state b i t error rate r ~ number of terms taken i n expansion */ /* Sample "codegen" Code Specification F i l e •define N 7 •define K 4 •define M 16 const char codename[] = "Hamming (7,4), c y c l i c , l+z+z*3"; const i n t gen[K] [N] = ( 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1 ); const i n t wt[N+l) = ( 1, 0, 0, 7, 7, 0, 0, 1 ); */ •include <stdio.h> •include <math.n> •include "codegen" double ul[N+l], vl[N+l], wl[N+l); /* Binomial c o e f f i c i e n t s . */ double binom(m,n) i n t m, n; ( i n t i ; double x=1.0; i f (n>=0 SS m>=0) ( i f (n<m) return 0.0; e l s e ( for (i=0; i<=m-l; i++) x = x * (double)(n-i)/(double)(m-i) ; return x; ) ) i f (n>=0 SS m<0) return 0.0; i f (n<0 SS m>=0) < for (i=0; i<=m-l; i++) x = x * (double)(m-n-l-i)/(double)(m-i); return parity(m)*x; ) i f (n<0 SS m<0) ( i f <n<m) return 0.0; else ( for (i=0; i<=n-m-l; i++) x = x * (double)(-m-l-i)/(double)(n-m-return parity(n-m)*x; ) ) ) /* Binary addition routine. */ in t x o r ( i , j ) i n t i , j ; ( i f (i = j) return 0; else return 1; ) /* ( - l ) ' l for integer 1. */ int parity(1) i n t 1; ( i f ((l/2)*2-l=0) return 1; else return -1; ) /* F i r s t root of x*2-(Q+h*q)*x+h*(Q-p)=0. */ double jp(pg,pb,h) double pg, pb, h; ( double b, c; b = -(1-pg+h*(1-pb)); c = h*(1-pg-pb); return (-b+sqrt(b*b-4*c))/2; ) /* Second root of xA2-(Q+h*q)*x+h*(Q-p)=0. */ double jm(pg,pb,h) double pg, pb, h; ( double b, c; b = -(1-pg+h*(1-pb)); c = h*(1-pg-pb); return (-b-sqrt(b*b-4*c))/2; I /* P r o b a b i l i t y of k O'a following a 1. * / double u(k,pg,pb,h) i n t k; double pg, pb, h; ( return ( (jp(pg,pb,h)+pg+pb-l)*pow(jp(pg,pb,h),(double)(k)) -(jm(pg,pb,h)+pg+pb-l)*pow(jm(pg,pb,h), (double) (k)) ) / (JP (pg,pb,h) - jm(pg,pb,h)) ; ) /* P r o b a b i l i t y of k O'a between l ' s . */ double v(k,pg,pb,h) i n t k; double pg, pb, h; ( return ( ((1-pb)*jp (pg,pb,h)+pg+pb-l)*pow(jp(pg,pb,h),(double)(k)) - ((1-pb)* jm(pg,pb,h)+pg+pb-1) *pow(jm(pg,pb,h),(double)(k)) *(1-h)/(jp(pg,pb,h)-jm(pg,pb,h)); ) /* P r o b a b i l i t y of k O'a preceding a 1. * / double w(k,pg,pb,h) i n t k; double pg, pb, h; ( return ( (jp(pg,pb,h)+pg+pb-l)*pow(jp(pg,pb,h), (double)(k)) -(jm(pg,pb,h)+pg+pb-l)*pow(jm(pg,pb,h), (double) (k)) ) * (Pg/ (pg+pb)) * (1-h) / (jp (pg, pb, h) - jm (pg, pb, h)) ; /* P r o b a b i l i t y of b i t pattern s. */ double prob(s) i n t a [ ) ; ( i n t i ; i n t c=0; char firat_one='Y'; double x=1.0; for (i=0; i<N; i++) ( i f (a[i] = 0) C++; el s e ( i f (firat_one = 'Y') ( x = x*wl[c); fi r s t _ o n e = 'N'; ) e l s e { x = x * v l [ c ] ; ) c = 0; ) ) x = x * u l [ c ] ; return x; /* The rth term i n the expanaion of P0(m,n)/pG i n P. */ double pOr(m,n,pg,pb,r) i n t m, n, r; double pg, pb; ( int b, c; double x=0.0; i f (r=0) ( i f (m=0) return 1.0; else return 0.0; } else i f (r=l> I i f (ro=0) return (l-n)*pg; i f (m>0 SS nxn) return 2*pg*pow(l-pb, (double)(m-1))*(1+((n-m-1)12.0)*pb); i f (III—n) return pg*pow(l-pb, (double)(n-l))/pb; ) else { i f (m=0) { return parity(r)*pow(pg, (double)(r))*binom(r,n-l); ) i f (m>0 tS nKn) ( for (b=0; b<=r; b++) x = x + (n-m-b+1)*(n-m-b) *binom(b,r)*binom(r-l,m+b-l)*pow(pb-l,(double) (b)); return pow(pg, (double) (r)) *pow(l-pb, (double) (m-r)) *binom(r-l,n-m)*x/(r*(n-m)); I i f (m—n) return 0.0; ) ) / * P0(m,n) accurate to order r i n P. */ double pO(m,n,pg,pb,r) in t m, n; double pg, pb; ( i n t 1; double x=0.0; for (1=0; K=r; 1++) x = x + p0r(m,n,pg,pb,l); return (pb/(pg+pb))*x; ) /* Sum of (l-h) Am*h A(i-m)*C(m,i)*vector[i) over i=m..n. */ double p(m,n,h, vector) i n t m, n; double h; double vector[] ; ( i n t i ; double x=0.0; f o r (i=m; i<=n; i++) x = x + binom(m,i)*pow(h, (double) (i-m))*vector[i]; return (((<x*pow(l-h,m/4.0))*pow(l-h,m/4.0))*pow(l-h,m/4.0))*pow(l-h,m/4.0)) v o i d main() ( char carry; i n t b, c, i , m, r; i n t k=0, n=0; i n t mes[K], cw[N]; double pg, pb, h; double total=0.0 , totala=0.0; double pO_o[N+l]; p r i n t f ( " % s \ n " , codename); p r i n t f C P = " ) ; s c a n f ( " % l f " , i p g ) ; p r i n t f ( " p = " ) ; s c a n f ( " % l f " , Spb); p r i n t f ( " h = " ) ; s c a n f ( " % l f " , *h); p r i n t f ( " r = " ) ; scanf("%d", t r ) ; for (m=0; m<=N; ro++) ( ul[m] = u(m,pg,pb,h); vl[m] = v (m,pg,pb,h) ; wl[m] = (pg/(pg+pb))*(l-h)*ul[m]; ) for (k=0; k<K; k++) mes[kl = 0; for (b=l; b<M; b++) ( carry = 'X'; for (k=0; k<K t t carry='Y'; k++) ( i f ((mes[k]=xor(mes[k],l)) = 1) carry = 'N'; ) for (n=0; n<N; n++) ( c = 0; for (k=0; k<K; k++) c = xor(c,mes[k]*gen[k][n]); cw[n] = c; ) t o t a l = t o t a l + prob(cw); ) for (m=0; m<=N; m++) p0_o(m] = pO (m,N,pg,pb,r) ; for ( i=l; i<=N; i++) t o t a l a = t o t a l a + (wt[i)/binom(i,N))*p(i,N,h ,p0_o); p r i n t f ( " P u = %.4e ", t o t a l ) ; printf("Pua = %.4e ", total a ) ; p r i n t f ( " % % d i f f . = %.4e\n", 1 0 0 *(totala-total)/total); return; TJ (D x" O o o ( i f (randl) > h*MAXINT) mc_crc.c /* mc_crc.c */ ( z [ i ] = 1; /* Computes the approximate probability of undetected error for CRC codes fl a g = 1; i n the G i l b e r t channel by Monte Carlo simulation and the 99% confidence ) l i m i t s . */ else z t i ] = 0; /* include f i l e required: i f (randO > (1-pb)*MAXMT) "codepar" - see below for format state = 'G'; /* parameters prompted for at runtime: ) P - G to B t r a n s i t i o n probability ) p = B to G t r a n s i t i o n probability 1 h = B state b i t error rate N = number of samples */ /* Binary addition routine. */ i n t xor(i,j) int i , j ; /* Sample "codepar" Code Specification F i l e ( •define N 21 i f ( i = j) •define P 16 return 0; else const char codename[] = "CRC-16, l+z*2+z A15+z A16"; return 1; const i n t g[P+l) = (1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1); */ ) /* Check to see i f vector z[N] i s a non-zero codeword. */ i n t check() •include <stdio.h> ( •include <math.h> in t i , j , q; •include <values.h> int sr[P); •include "codepar" i f (flag = 0) return 0; const seed[4] = (1, 2, 3, 4) ; /» random number seed values */ for (i=0; i<P; i++) const double d99 = 2.575829; s r [ P - l - i ) = z [ i ] ; for (i=P; i<N; i++) i n t z[N]; ( i n t f l a g ; q = s r [ P - l ] ; for (j=P-l; j>0; j — ) s r [ j ) = x o r ( s r [ j - l ] , g [ j ] * q ) ; /* Generate a random vector z[N), set f l a g to 1 i f i t ' s nonzero. */ sr[0]=xor(z[i],g[0]*q); ) void generate (pg,pb,h) double pg,pb,h; i for (i=0; i<P; i++) t char state; i f (sr[i] != 0) i n t i ; return 0; f l a g = 0; } return 1; i f (rand() > (pb/(pg+pb))*MAXINT) 1 state = 'B'; e l s e state = ' G' ; void main() for (i=0; i<N; i++) { < int sample, fseed, a, t ; i f (state == 'G') int m=0; ( double pg, pb, h, p; z[i]=0; double ym99, yp99; i f (randO > (1-pg)*MAXINT) state = 'B'; p r i n t f ( " % s \ n " , codename); ) p r i n t f ( " k = %d\n", N-16); e l s e printf("P = " ) ; scanf("%1£", Ipg); p r i n t f ( " p = " ) ; s c a n f ( " % l f " , Spb); p r i n t f ( " h = " ) ; s c a n f ( " % l f " , Sh); p r i n t f ( " N = " ) ; scanf("%d", Ssample); for (t=0; t<4; t++) ( srand(seed[t]); for (s=0; s<sample/4; s++) { generate(pg,pb,h); i f (check () = 1) m = m + 1; ) 1 p = ((double)(m)) / sample; ym99 = p * (1+((d99*d99)/(2*p*sample))*(1-sqrt((4*p*sample)/(d99*d99)+l))); yp99 = p * (l+((d99*d99)/(2*p»sample))*(l+sqrt((4»p*sample)/(d99*d99)+l))); pr i n t f ( " P u = %e\n", p); printf("lower error for 99%* confidence = %e\n", p-ym99); printf("upper error for 99%* confidence = %e\n", yp99-p); return; pmnc.map # pmn_c.map t C a l c u l a t e s PO(m,n) u s i n g t h e e x a c t c l o s e d fo rm e x p r e s s i o n g i v e n by # Cuperman and t h e n P ( m , n ) . t p r o c e d u r e s : # s - computes a t e r m i n t h e sum f o r PO (m, n) # s a i l - computes a l l t e r m s i n t h e sum f o r P0(m,n) i pOc - computes P0(m,n) # p O c a l l - computes P0(m,n) f o r a l l m # p c - computes P(m, n) # p c a l l - computes P(m,n) f o r a l l m a := p r o c ( m , n , p g , p b , j ) (n-m-j )*pb+m*(1-pg) e n d ; s := p r o c ( m , n , p g , p b , j ) ( p g / ( p g + p b ) ) * ( ( 1 - p b ) * (m- j - 1 ) ) * ( ( 1 - p g ) A ( n - m - j - 2 ) ) * ( ( p g + p b - 1 ) A j ) * ( a ( m , n , p g , p b , j ) A 2 - p b * a ( m , n , p g , p b , j ) - m * ( 1 - p g - p b ) * ( 1 - p g ) ) * ( n - j - 2 ) ! / ( m * j ! * ( m - j - l ) ! * ( n - m - j ) !) e n d ; s a i l := p r o c ( m , n , p g , p b ) f o r j f rom 0 t o m-1 do p r i n t ( j , s ( m , n , p g , p b , j ) ) ; o d ; e n d ; pOc := p r o c ( m , n , p g , p b ) i f m=0 t h e n ( p b / ( p g + p b ) ) * ( ( 1 - p g ) A ( n - l ) ) e l i f n x n t h e n s u m ( ' s ( m , n , p g , p b , j ) ' , ' j ' = 0 . .min (m-1 ,n -m) ) e l s e ( p g / ( p g + p b ) ) * ( ( l - p b ) A ( m - 1 ) ) f i e n d ; f i l l a r r a y := p r o c ( n , p g , p b ) p O c a r r : = a r r a y ( 0 . . n ) ; f o r m f rom 0 t o n do p O c a r r [m] :=pOc (m,n ,pg ,pb) ; o d ; e n d ; p O c a l l := p r o c ( n , p g , p b ) f i l l a r r a y ( n , p g , p b ) ; f o r m f r o m 0 t o n do p r i n t (m, p O c a r r [m]) ; o d ; e n d ; p c := p r o c ( m , n , p g , p b , h ) f i l l a r r a y ( n , p g , p b ) ; y := 0 . 0 ; f o r i f rom m t o n do y := y + b i n o m i a l ( i , m ) * ( h A ( i - m ) ) * p 0 c a r r [ i ] ; o d ; y * ( ( l - h ) ' m ) ; e n d ; p e l := p r o c ( m , n , p g , p b , h ) y := 0 . 0 ; p c a l l f o r i f rom m t o n do y := y + b i n o m i a l ( i , m ) * ( h A ( i - m ) ) * p 0 c a r r [ i ] ; o d ; y * ( ( 1 - h ) A m ) ; end ; := p r o c ( n , p g , p b , h ) f i l l a r r a y ( n , p g , p b ) ; f o r m from 0 t o n do p r i n t ( m , p c l ( m , n , p g , p b , h ) ) ; o d ; end; "O (D r> Q. X" O o CO od; pmn_s.map end; # pmn s.map p := proc(m,n,pg,pb,h,r) • Compute P0(m, n) by series expansion in P and then P(m,n). fillarray(n,pg,pb,r); y := 0.0; 1 procedures: for i from m to n do • pOr - computes rth term of series expansion for P0(m,n) y := y + binomial(i,m)*(h*(i-m))*p0arr[i); * pOrall - computes a l l terms of series expansion for P0(m,n) up to r od; • pO - computes P0(m,n) up to order r accuracy y*((l-h)-m); 1 pOall - computes P0(m,n) up to order r accuracy for a l l m end; •  - computes P(m,n) to order r accu acy t pal l - computes P(m, n) to order r accuracy for a l l m pi := proc(m,n,pg,pb,h, r) y := 0.0; for i from m to n do pOr : = proc(m,n,pg,pb,r) y := y + binomial(i,m)*(h*(i-m))*p0arr[i]; i f r=0 then od; i f m=0 then 1 y*((l-h)*m); else 0 end; f i e l i f r=l then pall := proc(n,pg,pb,h,r) i f m=0 then (l-n)*pg fillarray(n,pg,pb,r); e l i f m=n then (pg/pb)*(1-pb)*(n-l) for m from 0 to n do else 2*pg*(1+(n-m-l)*pb/2)*(1-pb)A(m-l) print (m, p i (m,n,pg,pb,h, r)) ; f i od; else end; i f m=0 then ((-pg)*r)"binomial(n-l,r) e l i f m=n then 0 else ((pg"r)*((l-pb)A(m-r))"binomial(n-m, r-1)/(r*(n-m))) *sum(' (n-m-b+1)*(n-m-b)'binomial(r,b)'binomial(m+b-1,r-1) *(pb-l)'b','b'=0..r) f i f i end; pOrall := proc(m,n,pg,pb,r) for j from 0 to r do print (pOr (m,n,pg,pb, j)) ; od; end; pO := proc(m,n,pg,pb,r) x := 0.0; for 1 from 0 to r do x := x + pOr(m,n,pg,pb,l) ; od; x*(pb/(pg+pb)) ; end; f i l larray := proc(n,pg,pb,r) pOarr:=array(0..n); for m from 0 to n do pOarr [m] :=p0 (m, n,pg, pb, r); od; end; pOall := proc(n,pg,pb,r) fillarray(n.pg.pb, r); for m from 0 to n do print (m, pOarr [m]) ; 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items