library( binfunest)
set.seed( 31394)

Bit Error Rate Functions

The package include a number of functions that report the theoretical bit error rate (BER) of many common modulation formats (Proakis, John G. 1995), with the signal-to-noise ratio (SNR) specified as the energy in a single bit divided by the noise power in a one Hertz bandwidth (E_b/N₀). See `help(“Theoretical”) for a complete list.

plot( QPSKdB, 3, 20, log="y", ylim=c( 1e-10, 0.5), panel.first = grid(),
      main="Modulation Performance Curves", xlab="SNR in dB",
      ylab="BER")
curve( QAMdB.8.star, 3, 20, col="blue", add=TRUE)
curve( PSQPSKdB, 3, 20, col="red", add=TRUE)
curve( QAMdB.16, 3, 20, col="green", add=TRUE)
legend( "bottomleft", legend=c( "QPSKdB", "QAMdB.8.star", "PSQPSKdB",
                                "QAMdB.16"),
        lty=c( 1, 1, 1, 1),
        col=c("black", "blue", "red", "green"))

Offset and a Back-to-Back SNR

In high performance communications systems it is common for both the transmitter (Tx) and receiver (Rx) to add noise to the signal. Unless you have a perfect example of one or the other, it is impossible to measure the source: both the Tx and Rx will contribute. This has the effect of having a minimum bit error rate (BER) even when no noise is added to the signal, called the Back-to-Back (B2B) BER, and is often expressed in Decibels as the signal-to-noise ratio (SNR) which would the B2B BER in a perfect Rx.

The other common effect is to add an offset to the SNR. That is, the system may always perform as though the SNR were a few Decibels less.

We can take a BER function and add the offset and B2B BER to it as follows. The linear SNR is γ = E_b/N₀, but in the B2B case the input noise N₀ is zero, yet a finite error rate is still observed. This is modeled as γ = E_b/(N₀ + β) where β is 1/B2B.

$$ P_{B2B}(\gamma) = P\left[ O\left( \frac{1}{1/\gamma + 1/B}\right)\right]$$

where O is the linear offset (i.e., undB( offset)), B is the linear B2B SNR, and P_B2B(γ) is the probability of error at a linear SNR of γ.

Creating BER Functions with B2B & Offset

If we have a function P_dB(s_dB) of the SNR in Decibels, we can express the B2B and offset in decibels too.

P_e(s_dB, B_dB, O_dB) = BER(−dB(undB(−s_dB) + undB(−B_dB) − O_dB)

The B2BQ package provides a function factory to convert any such function P_dB into a function with B2B and Offset in Decibels B2BConvert.

QPSKdB.B2B <- B2BConvert( QPSKdB)
O1 <- 3
B1 <- 16
s <- 0:20
plot( s, y=QPSKdB( s), typ="l", log="y", ylim=c( 1e-10, 0.5), 
      panel.first = grid(), main="QPSK Performance Limits", xlab="SNR in dB",
      ylab="BER")
lines( s, y=QPSKdB.B2B( s, Inf, O1), col="blue", lwd=2)
lines( s, y=QPSKdB.B2B( s, B1, O1), col='red')
abline( h=QPSKdB.B2B( B1, +Inf, O1), col='black', lty=2)
legend( "bottomleft", legend=c( "Theory", "Offset", "Offset + B2B", "B2B"),
        lty=c( 1, 1, 1, 2), lwd=c(1, 2, 1, 1), 
        col=c( 'black', 'blue', 'red', 'black'))

Dotted line shows B2B limiting BER.

Using `nls`

Let’s generate some example data and use nls to fit the curve to the data.

N <- 1000000L
(r <- rbinom( length( s), N, QPSKdB.B2B( s, B1, O1)))
#>  [1] 161193 134359 107870  83728  61830  43368  28246  17266   9817   5083
#> [11]   2195   1024    393    134     44     11      7      1      0      0
#> [21]      0
bplot1 <- function( s, r, N, O, B, ylim=c( 1e-10, 0.5)) {
  plot( s, y=QPSKdB.B2B( s, B, O), col='red', log='y', type='l', ylim=ylim,
        main="QPSK Performance Limits", xlab="SNR in dB",
      ylab="BER", panel.first = grid())
  points( s, r/N)
  lines( s, y=QPSKdB( s))
  lines( s, y=QPSKdB.B2B( s, Inf, O), col="blue", lwd=2)
  lines( s, y=QPSKdB.B2B( rep( B1, length( s)), Inf, O1), col='black', lty=2)
  legend( "bottomleft",
        legend=c( "Data", "Theory", "Offset", "Offset + B2B", "B2B"),
        lty=c( NA, 1, 1, 1, 2), lwd=c(NA, 1, 2, 1, 1), pch=c( 1, NA, NA, NA, NA),
        col=c( 'black', 'black', 'blue', 'red', 'black'))
}
bplot1( s, r, N, O1, B1)

Samples above 17 dB resulted in zero errors (in this example) and so are not plotted.

nls creates a non-linear least squares fit to estimate the parameters.

df1 <- data.frame( Errors=r, Trials=rep( N, length( s)), SNR=s)
nls.fit1 <- nls( Errors/Trials ~ QPSKdB.B2B( SNR, b, o), data=df1, 
                 start=list( b=10, o=0))
summary( nls.fit1)
#> 
#> Formula: Errors/Trials ~ QPSKdB.B2B(SNR, b, o)
#> 
#> Parameters:
#>    Estimate Std. Error t value Pr(>|t|)    
#> b 15.948539   0.052466     304   <2e-16 ***
#> o  2.992966   0.002874    1041   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 8.646e-05 on 19 degrees of freedom
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 6.26e-08
bplot1(s, r, N, O1, B1)
c <- coef( nls.fit1)
lines( s, y=QPSKdB.B2B( s, c['b'], c['o']), col='green')

Non-linear least squares fit (green line) of the line to the data.

Those are remarkably close estimates. However, we have fit the function to the error rate, not the error counts. There are zeros in the data, and this doesn’t match the function at high SNRs very well. The following will perform the fit without the zero points.

nls.fit2 <- nls( Errors/Trials ~  QPSKdB.B2B( SNR, b, o), 
                 data=subset( df1, Errors > 0), start=list( b=10, o=0))
summary( nls.fit2)
#> 
#> Formula: Errors/Trials ~ QPSKdB.B2B(SNR, b, o)
#> 
#> Parameters:
#>    Estimate Std. Error t value Pr(>|t|)    
#> b 15.948539   0.057174   278.9   <2e-16 ***
#> o  2.992966   0.003132   955.5   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 9.422e-05 on 16 degrees of freedom
#> 
#> Number of iterations to convergence: 7 
#> Achieved convergence tolerance: 9.639e-08
bplot1(s, r, N, O1, B1)
c <- coef( nls.fit2)
lines( s, y=QPSKdB.B2B( s, c['b'], c['o']), col='green')

Non-linear fit (green line) excluding the zeros.

That was unsatisfying; the answer is almost identical to the version including the zeros. Plotting the data on a linear scale indicates why.

 plot( s, y=QPSKdB.B2B( s, B1, O1), col='red', type='l', panel.first = grid(),
       main="QPSK Performance Limits", xlab="SNR in dB", ylab="BER")
  points( s, r/N, panel.first = grid())
#> Warning in plot.xy(xy.coords(x, y), type = type, ...): "panel.first" is not a
#> graphical parameter
  lines( s, y=QPSKdB( s))
  lines( s, y=QPSKdB.B2B( s, Inf, O1), col="blue", lwd=2)
  lines( s, y=QPSKdB.B2B( s, c['b'], c['o']), col='green')
  legend( "bottomleft",
        legend=c( "Data", "Theory", "3 dB Offset", "Offset + B2B", "Fit"),
        lty=c( NA, 1, 1, 1, 1), lwd=c(NA, 1, 2, 1, 1), pch=c( 1, NA, NA, NA, NA),
        col=c( 'black', 'black', 'blue', 'red', 'green'))

Linear plot of performance showing how close to zero the estimates lie.

The error rate is so close to zero past 10 dB that a least squares error minimization is not affected by any fit errors. A better technique might be to fit the inverse “Q” function of both sides of the equation. The Q function is Markum’s Q function and is given by

$$ Q(s) = \frac{1}{2} \mathrm{Erfc} \left( s / \sqrt{2}\right) $$

In R, this is simply pnorm( x, lower.tail=FALSE), although the package provides the function Q_(x) for convenience. The inverse “Q” function returns the argument of Q_(x) which results in the supplied error rate, conventionally returning the number in Decibels (thus the idea of back-to-back Q in Decibels). This package provides the function defined as Q_Inv(x) = 2 * dB( -qnorm( x)). Since Q_Inv of zero is infinite, we have to exclude zero results again.

nls.fit3 <- nls( Q_Inv( Errors/Trials) ~  Q_Inv( QPSKdB.B2B( SNR, b, o)), 
                 data=subset( df1, Errors > 0), start=list( b=10, o=0))
summary( nls.fit3)
#> 
#> Formula: Q_Inv(Errors/Trials) ~ Q_Inv(QPSKdB.B2B(SNR, b, o))
#> 
#> Parameters:
#>   Estimate Std. Error t value Pr(>|t|)    
#> b 15.92088    0.09440   168.7   <2e-16 ***
#> o  2.99072    0.02542   117.6   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.06776 on 16 degrees of freedom
#> 
#> Number of iterations to convergence: 4 
#> Achieved convergence tolerance: 4.165e-08
bplot1(s, r, N, O1, B1)
c <- coef( nls.fit3)
lines( s, y=QPSKdB.B2B( s, c['b'], c['o']), col='green')

Results of non-linear fit (green line) to Q function. Zeros are still filtered out.

Maximum Likelihood Estimation

In order to use all of the data, and to derive statistical significance from the estimates, we must turn to maximum likelihood estimation. This will model the parameters that most likely created the observed data (including the zeros).

Log Likelihood of observations

The likelihood is the probability that the specific observation took place. This is the product of the probability of each observation, and we will maximize the sum of log of the probabilities. The dbinom function returns exactly this:

dbinom( r, N, QPSKdB.B2B( s, B1, O1), log=TRUE)
#>  [1] -6.94159646 -6.75102233 -7.28795647 -6.56415934 -6.45336248 -6.29411569
#>  [7] -6.52783185 -5.90179937 -5.88238628 -6.01136681 -8.43372668 -5.37619530
#> [13] -4.34938113 -3.40025874 -2.83464924 -2.31452040 -2.84437356 -1.01822454
#> [19] -0.37789410 -0.12626314 -0.04583595

It is easy to see that the high-SNR observations (the last ones: the SNR ranges from zero to 20 with the samples numbered 1 to 21) have the highest likelihood. This is because when the SNR is high, zero errors is about the only likely event, while at a lower SNR, anything the the range of $N p \pm 2\sqrt{N p}$ is pretty likely (where N is the number of samples 1,000,000, and p is the probability of an error; still assuming p ≪ 1 here), so each individual result has a small likelihood. For example, at an SNR of 1 dB (2nd item) the BER is 0.13 and the expected value is 134,364. While this is the peak likelihood, the probability is spread around plus or minus about 733 values.

Using `optim`

We can use optim to find maximum likelihood estimates of the offset and B2B. We have to cast the parameters as a vector and create a function llsb to return the log-likelihood of the observations. optim will naturally find minimums, so the negative of the sum is used. Sorting the log likelihood before the sum will ensure that precision isn’t lost, although such loss of precision is “unlikely.”

llsb <- function( par) 
  -sum( sort( dbinom( r, N, QPSKdB.B2B( s, par[1], par[2]), log=TRUE),
              decreasing=TRUE))
mle1 <- optim( par=c( b2b=20, off=0), llsb, hessian=TRUE)
mle1$par
#>       b2b       off 
#> 15.937961  2.992359
(mle1s <- sqrt( diag( solve( mle1$hessian))))
#>         b2b         off 
#> 0.060770870 0.006420547

Those are pretty close estimates (although not as close as the prior fits). The diagonal of the inverse of the Hessian matrix is an estimate of the variance of the parameters, so their square root is the standard deviation.

Below is a simple function that will take a function created by B2BConvert and estimate the two parameters. (A few extra parameters are added for possible later use.)

mlef <- function( N, s, r, f, ..., startpar=c( b2b=20, off=0),
                    method="Nelder-Mead", lower, upper, control) {
  stopifnot( length( N) == length( s) && length( s) == length( r))
  ll <- function( par) {
    -sum( dbinom( r, N, f( s, par[1], par[2], ...), log=TRUE))
  }
  res <- optim( startpar, ll, hessian=TRUE)
}

mle2 <- mlef( rep( N, length( s)), s, r, QPSKdB.B2B)
mle2$par
#>       b2b       off 
#> 15.937961  2.992359
sqrt( diag( solve( mle2$hessian)))
#>         b2b         off 
#> 0.060770870 0.006420547

The output is in $par, the first being the B2B SNR, and the second being the offset. Note the names come from the startpar vector.

Here is a plot with the original data and the estimated line. Note the estimated line lies right on top of the the “Offset + 15 dB B2B”.

bplot1(s, r, N, O1, B1)
lines( s, QPSKdB.B2B( s, mle2$par[1], mle2$par[2]), col='green', lty=4)
legend( "bottomleft",
        legend=c( "Data", "Theory", "3 dB Offset", "Offset + 15dB B2B",
                  "Estimated"),
        lty=c( NA, 1, 1, 1, 4), pch=c( 1, NA, NA, NA, NA),
        col=c( 'black', 'black', 'blue', 'red', 'green'))

Estimated as green dashed line (dashing allows the red line to peek through).

Using `mle`

The stats4::mle function can do a better job, and it does not need the parameters cast into a list.

require( stats4)
#> Loading required package: stats4
llsb2 <- function( b2b, off) 
   -sum( dbinom( r, N, QPSKdB.B2B( s, b2b, off), log=TRUE))
mle3 <- mle( llsb2, start=c( b2b=20, off=0), nobs=length(s))
summary( mle3)
#> Maximum likelihood estimation
#> 
#> Call:
#> mle(minuslogl = llsb2, start = c(b2b = 20, off = 0), nobs = length(s))
#> 
#> Coefficients:
#>      Estimate  Std. Error
#> b2b 15.939294 0.060817557
#> off  2.992328 0.006422466
#> 
#> -2 log L: 190.0448

stats4 provides other functions to work with the result:

logLik( mle3)
#> 'log Lik.' -95.02241 (df=2)
vcov( mle3)
#>              b2b          off
#> b2b 0.0036987752 3.170244e-04
#> off 0.0003170244 4.124807e-05
sqrt(diag(vcov( mle3)))
#>         b2b         off 
#> 0.060817557 0.006422466
plot(profile( mle3), absVal = FALSE)

confint( mle3)
#> Profiling...
#>         2.5 %    97.5 %
#> b2b 15.821521 16.059970
#> off  2.979737  3.004913

The covariance indicates that there is a relationship between the two parameters, although it is hard to determine how strong. We’ll look into this further below.

Using `mleB2B`

The same result can be reached easier with this package’s mleB2B function, which facilitates the data being hosted in a data frame or list as well as vectors (the dataframe form is shown here):

df <- data.frame( Errors=r, SNR=s)
est1 <- mleB2B( data=df, Errors="Errors", N=N, f=QPSKdB.B2B, 
                fparms=list( x="SNR"), start=c(b2b=20, offset=0))
summary( est1)
#> Maximum likelihood estimation
#> 
#> Call:
#> stats4::mle(minuslogl = function (b2b = 20, offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QPSKdB(-dB((undB(-x) + undB(-B2B))) - offset))(x = 0:20, b2b, 
#>     offset), log = TRUE))), start = c(b2b = 20, offset = 0), 
#>     method = "Nelder-Mead", nobs = 21L)
#> 
#> Coefficients:
#>         Estimate  Std. Error
#> b2b    15.937961 0.060770870
#> offset  2.992359 0.006420547
#> 
#> -2 log L: 190.0468

Correlation between B2B and Offset

The correlation matrix can be computed as follows:

c <- vcov( est1)
v <- diag(1 / sqrt( diag( c)))
v %*% c %*% v
#>           [,1]      [,2]
#> [1,] 1.0000000 0.8115133
#> [2,] 0.8115133 1.0000000
rm( c, v)

The package provides a function, as we’ll find it useful.

cor.mle( est1)
#>           [,1]      [,2]
#> [1,] 1.0000000 0.8115133
#> [2,] 0.8115133 1.0000000

A correlation of 1 is rather large. This means that the two random variables of the estimated B2B and offset are not independent. This can be demonstrated:

L <- 100
tmp <- function(K){
  r <- rbinom( length( s), N, QPSKdB.B2B( s, B1, O1))
  coef( mleB2B( Errors=r, N=N, f=QPSKdB.B2B, 
                fparms=list( x=s), start=c(b2b=20, offset=0)))
}
samps <- vapply(1:L, tmp, FUN.VALUE=c("B2B"=0, "Offset"=0)) 
plot( t(samps))

100 Samples of parameter estimates.

Very High B2B Q

It is common for many simulated systems or lower data-rate systems to reach almost theoretical performance, perhaps with an offset, but no, or almost no, detectable B2B. That is, the system can be run back-to-back indefinitely and generate no errors. In this case, estimating a B2B might corrupt a simple offset estimate. The following has no visible B2B Q.

O2 <- 3
B2 <- 24
N <- 1000000
(r2 <- rbinom( length( s), N, QPSKdB.B2B( s, B2, O2)))
#>  [1] 159138 131179 104411  79601  57103  38680  23711  13334   6439   2786
#> [11]    957    268     53      8      0      0      0      0      0      0
#> [21]      0
mle4 <- mleB2B(  Errors=r2, N=N, f=QPSKdB.B2B,
                 fparms=list( x=s), start=c(b2b=20, offset=0))
summary(mle4)
#> Maximum likelihood estimation
#> 
#> Call:
#> stats4::mle(minuslogl = function (b2b = 20, offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QPSKdB(-dB((undB(-x) + undB(-B2B))) - offset))(x = 0:20, b2b, 
#>     offset), log = TRUE))), start = c(b2b = 20, offset = 0), 
#>     method = "Nelder-Mead", nobs = 21L)
#> 
#> Coefficients:
#>         Estimate  Std. Error
#> b2b    23.977009 0.429759639
#> offset  3.003006 0.006813417
#> 
#> -2 log L: 154.5195
mle4coef <- coef(mle4)
plot(  s, r2/N, log='y', panel.first = grid(), ylim=c(1e-14, 0.5))
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 7 y values <= 0 omitted from
#> logarithmic plot
lines( s, y=QPSKdB( s), col='black')
lines( s, y=QPSKdB.B2B( s, Inf, O1), col="blue", lwd=2)
lines( s, y=QPSKdB.B2B( s, B2, O2), col="red")
lines( s, y=QPSKdB.B2B( s, mle4coef[1],  mle4coef[2]), col="green")
legend( "bottomleft",
        legend=c( "Data", "Theory", "Offset", "Offset + B2B",  
                  "Estimated"),
        lty=c( NA, 1, 1, 1), lwd=c(NA, 1, 2, 1, 1), 
        col=c( 'black', 'black', 'blue', 'red', 'green'),
        pch=c( 1, NA, NA, NA, NA))

24 dB B2B Q and 3 dB Offset estimates.

The plot shows the data running down the Offset line, but it is quite unconvincing that the B2B estimate is correct. Typically the B2B parameter will have a very large variance in these instances (this does not), and multiple runs will demonstrate this. Also, the correlation between the two parameters is now 0.82. Ignoring our guilty knowledge of how the data was generated, we might want to test the hypothesis that the B2B Q is really infinite (this is a possibility). We can estimate the offset only using the fixed option of mle or simply setting the B2B parameter. Note that when optimizing over a single variable, a few other changes are required.

(mle5 <- mleB2B( Errors=r2, N=N, f=QPSKdB.B2B, method="Brent", 
                 fparms=list( x=s, B2B=+Inf), start=c( offset=0), 
                 lower=-6, upper=10))
#> 
#> Call:
#> stats4::mle(minuslogl = function (offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QPSKdB(-dB((undB(-x) + undB(-B2B))) - offset))(x = 0:20, B2B = Inf, 
#>     offset), log = TRUE))), start = c(offset = 0), method = "Brent", 
#>     nobs = 21L, lower = -6, upper = 10)
#> 
#> Coefficients:
#>  offset 
#> 3.05971

AIC(mle4, mle5)
#>      df      AIC
#> mle4  2 158.5195
#> mle5  1 262.0605

The substantially lower AIC for mle4 indicates that mle5 is not a statistically better estimate. Again ignoring our guilty knowledge that the B2B Q is really 24 dB, we might now generate some more data. A 100 Gbps system that we are willing to test for an hour will generate 3.6E+14 bits, and thus could be used to test a system (resulting in a BER of 3.7E-15, and an expected value of 1.3 errors).

Adding Observations

The following will simulate an additional run at 19 dB with many more samples, and add a single error a to the current observations. We need to make a vector of the number of trials now as they are not all the same. The plot below has an extended range.

N2 <- 100e9*3600
Nvec <- c( rep(N, length(s)), N2)
s2 <- c(s, 19)
(r4 <- c( r2, 1))
#>  [1] 159138 131179 104411  79601  57103  38680  23711  13334   6439   2786
#> [11]    957    268     53      8      0      0      0      0      0      0
#> [21]      0      1
mle6 <- mleB2B(  Errors=r4, N=Nvec, f=QPSKdB.B2B,
                 fparms=list( x=s2), start=c(b2b=20, offset=0))
summary( mle6)
#> Maximum likelihood estimation
#> 
#> Call:
#> stats4::mle(minuslogl = function (b2b = 20, offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QPSKdB(-dB((undB(-x) + undB(-B2B))) - offset))(x = c(0, 1, 2, 
#> 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
#> 20, 19), b2b, offset), log = TRUE))), start = c(b2b = 20, offset = 0
#> ), method = "Nelder-Mead", nobs = 22L)
#> 
#> Coefficients:
#>         Estimate  Std. Error
#> b2b    24.050326 0.345298503
#> offset  3.003859 0.005947601
#> 
#> -2 log L: 156.6186
mle6coef <- coef(mle6)
plot(  s2, r4/Nvec, log='y', panel.first = grid(), ylim=c(1e-15, 0.5))
#> Warning in xy.coords(x, y, xlabel, ylabel, log): 7 y values <= 0 omitted from
#> logarithmic plot
lines( s2, y=QPSKdB( s2), col='black')
lines( s2, y=QPSKdB.B2B( s2, Inf, O1), col="blue", lwd=2)
lines( s2, y=QPSKdB.B2B( s2, B2, O2), col="red")
lines( s2, y=QPSKdB.B2B( s2, mle6coef[1],  mle6coef[2]), col="green")
legend( "bottomleft",
        legend=c( "Data", "Theory", "Offset", "Offset + B2B",  
                  "Estimated"),
        lty=c( NA, 1, 1, 1, 1), lwd=c(NA, 1, 2, 1, 1), 
        col=c( 'black', 'black', 'blue', 'red', 'green'),
        pch=c( 1, NA, NA, NA, NA))

$24 dB B2B Q and 3 dB Offset estimates with additional data point.$

24 dB B2B Q and 3 dB Offset estimates with additional data point.

Note that the B2B estimate increased a bit, and the std error decreased a bit. The example below will show the same process with an infinite B2B, showing that when it truly is infinite, the AIC of an infinite model is lower.

(r3 <- rbinom( length( s), N, QPSKdB( s - O2)))
#>  [1] 158099 130613 103985  78763  56182  37374  22964  12536   6116   2396
#> [11]    737    170     25      5      1      0      0      0      0      0
#> [21]      0
mle7 <- mleB2B( N=N, Errors=r3, f=QPSKdB.B2B, fparms=list( x=s),
                start=c(b2b=20, offset=0))
summary( mle7)
#> Maximum likelihood estimation
#> 
#> Call:
#> stats4::mle(minuslogl = function (b2b = 20, offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QPSKdB(-dB((undB(-x) + undB(-B2B))) - offset))(x = 0:20, b2b, 
#>     offset), log = TRUE))), start = c(b2b = 20, offset = 0), 
#>     method = "Nelder-Mead", nobs = 21L)
#> 
#> Coefficients:
#>         Estimate   Std. Error
#> b2b    55.012504 93.667567141
#> offset  3.000775  0.004051487
#> 
#> -2 log L: 162.4777
cor.mle(mle7)
#>           [,1]      [,2]
#> [1,] 1.0000000 0.2350425
#> [2,] 0.2350425 1.0000000
mle8 <- mleB2B( N=N, Errors=r3, f=QPSKdB.B2B, fparms=list( x=s, B2B=+Inf),
                start=list(offset=0), method="Brent", lower=-6, upper=10)
summary( mle8)
#> Maximum likelihood estimation
#> 
#> Call:
#> stats4::mle(minuslogl = function (offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QPSKdB(-dB((undB(-x) + undB(-B2B))) - offset))(x = 0:20, B2B = Inf, 
#>     offset), log = TRUE))), start = list(offset = 0), method = "Brent", 
#>     nobs = 21L, lower = -6, upper = 10)
#> 
#> Coefficients:
#>      Estimate  Std. Error
#> [1,] 3.000947 0.003938021
#> 
#> -2 log L: 162.4718


AIC( mle7, mle8)
#>      df      AIC
#> mle7  2 166.4777
#> mle8  1 164.4718
exp( (AIC( mle7) - AIC( mle6)) / 2)
#> [1] 18.71879

Real Data

We have some simulations of various constellations and their bit maps used with perfect detectors (i.e., the constellation points were generated, noise was added, and then the points were detected; there was no actual modulation/demodulation). The BERDFc data included with the package has a count of symbols and bits per symbol, so we need to make a data vector of the count of bits. The code below examines an 8-ary constellation which does not have a perfect gray code (i.e., there are adjacent symbols with more than one bit difference). This causes some bit errors to come in pairs, and breaks the binomial assumption.

Q8Sbits <- BERDFc[BERDFc$Name=="8QAM Star",] # Select one constellation
Q8Sbits$Nbits <- Q8Sbits$Bps * Q8Sbits$N # Add  a Nbits column
QAMdB.8.star.B2B <- B2BConvert( QAMdB.8.star) 
mle10 <- mleB2B( data=Q8Sbits, N="Nbits", Errors="BER", 
                 f=QAMdB.8.star.B2B , start=c( B2B=30, offset=0), 
                 fparms=list( x="SNR"))
summary( mle10)
#> Maximum likelihood estimation
#> 
#> Call:
#> stats4::mle(minuslogl = function (B2B = 30, offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QAMdB.8.star(-dB((undB(-x) + undB(-B2B))) - offset))(x = c(3, 
#> 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14), B2B, offset), log = TRUE))), 
#>     start = c(B2B = 30, offset = 0), method = "Nelder-Mead", 
#>     nobs = 12L)
#> 
#> Coefficients:
#>         Estimate  Std. Error
#> B2B    20.045714 0.128527021
#> offset -0.240504 0.004910512
#> 
#> -2 log L: 891.1396
cor.mle( mle10)
#>           [,1]      [,2]
#> [1,] 1.0000000 0.9012394
#> [2,] 0.9012394 1.0000000

This looks like the B2B may be finite, and this can be verified by such a fit.

mle11 <- mleB2B( data=Q8Sbits, N="Nbits", Errors="BER", 
                 f=QAMdB.8.star.B2B , method="Brent", 
                 fparms=list( x="SNR", B2B=+Inf), start=c( offset=0), 
                 lower=-6, upper=10)
summary( mle11)
#> Maximum likelihood estimation
#> 
#> Call:
#> stats4::mle(minuslogl = function (offset = 0) 
#> -sum(sort(stats::dbinom(Errors, N, (function (x, B2B = Inf, offset = 0) 
#> QAMdB.8.star(-dB((undB(-x) + undB(-B2B))) - offset))(x = c(3, 
#> 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14), B2B = Inf, offset), log = TRUE))), 
#>     start = c(offset = 0), method = "Brent", nobs = 12L, lower = -6, 
#>     upper = 10)
#> 
#> Coefficients:
#>         Estimate  Std. Error
#> [1,] -0.08920179 0.002145889
#> 
#> -2 log L: 2078.837
AIC( mle10, mle11)
#>       df       AIC
#> mle10  2  895.1396
#> mle11  1 2080.8373

That result indicates that the B2B is actually needed. We can see what is happening here by looking at the data.

Q8Sbits
#>         Name SNR Bps NoisePower        N    SER    BER    Nbits
#> 25 8QAM Star   3   3 0.79032060  1600000 230974 306269  4800000
#> 26 8QAM Star   4   3 0.62789043  1600000 158180 210328  4800000
#> 27 8QAM Star   5   3 0.49863525  1600000  98785 131498  4800000
#> 28 8QAM Star   6   3 0.39581559  1600000  54622  72968  4800000
#> 29 8QAM Star   7   3 0.31499330  1600000  26713  35513  4800000
#> 30 8QAM Star   8   3 0.25009840  1600000  10919  14550  4800000
#> 31 8QAM Star   9   3 0.19861226  1600000   3445   4604  4800000
#> 32 8QAM Star  10   3 0.15778619  1600000    896   1181  4800000
#> 33 8QAM Star  11   3 0.12526456  1600000    137    177  4800000
#> 34 8QAM Star  12   3 0.09954161  1600000     17     23  4800000
#> 35 8QAM Star  13   3 0.07893708  1600000      2      2  4800000
#> 36 8QAM Star  14   3 0.06276492 17600000      1      2 52800000

The 14 dB SNR observation has only one symbol error, but reports two bit errors. This is because the 8QAM Star constellation has a few adjacent symbols that have two bit changes, and that one symbol error doubles the estimated error rate. Therefore, this estimator indicates that the BER is not binomialy distributed.

Proakis, John G. 1995. Digital Communication. Third Edition. New York, NY: McGraw-Hill.

Estimating B2B Q and Offsets