ADC quantization efficiency calculations for radio astronomy

quantization_efficiency.py

"""
# quantization_efficiency.py

This script implements the formulas for quantization efficiency
as used in radio astronomy. The formulas appear in:

    Thompson, A. R., D. T. Emerson, and F. R. Schwab (2007), 
    "Convenient formulas for quantization efficiency",
    Radio Sci., 42, RS3022, doi:10.1029/2006RS003585.
    
    A. R. Thompson, J. M.  Moran, G. W. Swenson, 
    "Interferometry and Synthesis in Radio Astronomy",
    3rd Edition, Springer Open
    https://link.springer.com/book/10.1007%2F978-3-319-44431-4
    
The second reference (TMS), is frely available online under a CC-BY license.

The two main functions here are:
compute_adc_nq - compute the quantization efficiency for 
                 an ADC with a given input RMS
compute_nq     - Compute quantization efficiency for Q levels 
                 with level spacing eps = 1/rms
"""

import numpy as np
from scipy.special import erf
from numpy import exp, sum, sqrt

def nq_even(N, eps):
    """ Compute quantization efficiency for even number of levels Q
    
    Args:
        N (int): Even level number, N = Q/2, where Q is number of levels
        eps (float): epsilon, level spacing (1/rms)

    Notes:
        Computes TMS Equation 8.79, pg 334
    """
    m  = np.arange(1, N)
    
    S0 = sum(exp(-m**2 * eps**2 / 2.0))
    num = (2.0 / np.pi) * (0.5 + S0)**2
    
    S1 = 2.0*sum(m * erf(m * eps / sqrt(2.0)))
    den = (N-0.5)**2 - S1
    
    return num / den

def nq_odd(N, eps):
    """ Compute quantization efficiency for odd number of levels Q
    
    Args:
        N (int): Odd level number, N = (Q-1)/2, where Q is number of levels
        eps (float): epsilon, level spacing in sigma units (1/rms)

    Notes:
        Computes TMS Equation 8.83, pg 335
    """
    m = np.arange(1, N+1)
    
    S0 = sum(exp(-(m-0.5)**2 * eps**2 / 2.0))
    num = (2.0 / np.pi) * S0**2
    
    S1 = 2.0*sum((m-0.5) * erf((m-0.5) * eps / sqrt(2.0)))
    den = N**2 - S1
    
    return num / den

def compute_nq(Q, eps):
    """ Compute quantization efficiency for Q levels with level spacing eps
    
    Args:
        Q (int or np.array): number of levels. Eg Q=256 for an 8-bit ADC
        eps (float or np.array): epsilon, level spacing in sigma units (1/rms)
        
    Notes:
        * Computes nq via TMS Eqns 8.79 and 8.83 as required
        * Functions are vectorized so numpy arrays can be used as inputs
    """
    N = int(Q/2)
    if Q % 2:
        compute_nq_vec = np.vectorize(nq_odd)
    else:
        compute_nq_vec = np.vectorize(nq_even)
        
    out = compute_nq_vec(N, eps)

    if out.size == 1:
        return float(out)
    else:
        return out
        #return nq_odd(N, eps)
        #return nq_even(N, eps)

def compute_adc_nq(n_bit, rms):
    """ Compute quantization efficiency for ADC
    
    Args:
        n_bit: Number of ADC bits
        rms: root mean square input level in ADC counts
    
    Notes:
        This is essentially a wrapper of compute_nq. Applies the 
        quantization equations from TMS.
    """
    eps = 1.0/rms
    Q = 2**(n_bit)
    return compute_nq(Q, eps)


def test_compute_nq():
    """ Unit test for compute_nq function
    
    Notes:
        Test data is from Table 8.2 in TMS.
    """
    Q_N_eps_nq = [
        [3, 1, 1.224, 0.80983],
        [4, 2, 0.995, 0.88115],
        [8, 4, 0.586, 0.96256],
        [9, 4, 0.534, 0.96930],
        [16, 8, 0.335, 0.98846],
        [32, 16, 0.188, 0.99651],
        [256, 128, 0.0312, 0.99991]
    ]
    
    for Q, N, eps, nq in Q_N_eps_nq:
        out = compute_nq(Q, eps)
        assert np.isclose(nq, out)

if __name__ == "__main__":  
    # Run unit test
    test_compute_nq()
    
    # Print out a list of ADC bit depth, RMS and efficiency
    print("\n ADC quantization efficiency table: \n")
    print(" N_bits  |   RMS    | n_Q  ")
    print("-----------------------------")
    for n_bits in (2, 3, 4, 5, 6, 7, 8):
        for rms in (4, 8, 16, 32, 64):
            if rms > 2**n_bits:
                print("%8s | %8d | - " % (n_bits, rms))
            else:
                eff = compute_adc_nq(n_bits, rms)
                print("%8s | %8d | %0.5f" % (n_bits, rms, eff))
        print("-----------------------------")
    
    # plot quantization efficiency as a function of epsilon
    #     This recreates Fig. 8.11 in TMS (pg 337):
    #     Quantization efficiency as a function of the threshold spacing, epsilon, 
    #     in units equal to the rms amplitude, sigma.
    print("\n Plotting quantization efficiency curves...")
    import pylab as plt
    
    eps = np.linspace(0.001, 2, 101)
          	  
    for Q in (128, 32, 16, 8, 4):
    	eff = compute_nq(Q, eps)
    	plt.plot(eps, eff, label='$Q=%i$' % Q )
    
    plt.ylim(0.6, 1.01)
    plt.legend(loc=4, ncol=2)
    plt.ylabel("$\\eta_Q$")
    plt.xlabel("$\\epsilon$")
    plt.title("Quant eff $\\eta_Q$ as fn. of threshold spacing $\\epsilon$ in units of rms ampl. $\\sigma$.")
    plt.show()
    
    print("Done.")

Output when run is a plot of quantization efficiency:

And a list of efficiencies for ADC input RMS:

 ADC quantization efficiency table:

 N_bits  |   RMS    | n_Q
-----------------------------
       2 |        4 | 0.74076
       2 |        8 | -
       2 |       16 | -
       2 |       32 | -
       2 |       64 | -
-----------------------------
       3 |        4 | 0.87234
       3 |        8 | 0.77113
       3 |       16 | -
       3 |       32 | -
       3 |       64 | -
-----------------------------
       4 |        4 | 0.98020
       4 |        8 | 0.89020
       4 |       16 | 0.78261
       4 |       32 | -
       4 |       64 | -
-----------------------------
       5 |        4 | 0.99481
       5 |        8 | 0.98646
       5 |       16 | 0.89725
       5 |       32 | 0.78757
       5 |       64 | -
-----------------------------
       6 |        4 | 0.99482
       6 |        8 | 0.99869
       6 |       16 | 0.98846
       6 |       32 | 0.90033
       6 |       64 | 0.78987
-----------------------------
       7 |        4 | 0.99482
       7 |        8 | 0.99870
       7 |       16 | 0.99967
       7 |       32 | 0.98918
       7 |       64 | 0.90176
-----------------------------
       8 |        4 | 0.99482
       8 |        8 | 0.99870
       8 |       16 | 0.99967
       8 |       32 | 0.99991
       8 |       64 | 0.98947
-----------------------------