Example on redigitizing floating point values to scaled and offset integers

digitize_example.py

#!/usr/bin/env python3
import numpy as np
import matplotlib.pyplot as plt

def digitize(x, smin, smax, nbits, signed=False, zeromean=False):
    # Stats
    xmean = x.mean()
    xstd = x.std()

    if zeromean:
        xmean = 0.0
    
    xmin = 0
    xmax = (1 << nbits) - 1
    xoffset = xmean + smin * xstd
    xscale = (smax - smin) * xstd / (1 << nbits)

    # Convert to ints
    x_int = np.clip(np.floor((x - xoffset) / xscale), xmin, xmax)

    if signed:
        x_int += -(1 << (nbits - 1))
        xoffset += (1 << (nbits - 1)) * xscale

    # Offset by half a step
    xoffset += xscale / 2
        
    return x_int, xscale, xoffset

if __name__ == "__main__":    
    # Time series length (~1 s block at 5.12 us sampling)
    n = 48 * 4096
    data_sigma = 1.0

    # Generate zero mean data or not?
    generate_zero_mean_data = True

    # Digitization settings
    smin, smax = -3, 3
    nbits = 8
    signed = False
    zeromean = True
    
    # Generate time axis
    t = np.arange(n) * 5.12e-6

    # Generate data
    np.random.seed(4)
    if generate_zero_mean_data:
        x = np.random.normal(0.0, data_sigma, n)
    else:
        xr = np.random.normal(0.0, data_sigma, n)
        xi = np.random.normal(0.0, data_sigma, n)
        yr = np.random.normal(0.0, data_sigma, n)
        yi = np.random.normal(0.0, data_sigma, n)
        x = xr**2 + xi**2 + yr**2 + yi**2
        
    # Digitize
    x_int, xscale, xoffset = digitize(x, smin, smax, nbits, signed, zeromean)

    # Reconstruct
    x_float = x_int * xscale + xoffset
    
    # Difference
    xdiff = x - x_float

    print(f"Scale value:  {xscale:+.5f}\nOffset value: {xoffset:+.5f}")
    print(f"Input data,         mean: {x.mean():+.5f}, standard deviation: {x.std():+.5f}")
    print(f"Reconstructed data, mean: {x_float.mean():+.5f}, standard deviation: {x_float.std():+.5f}")
    print(f"Difference,         mean: {xdiff.mean():+.5f}, standard deviation: {xdiff.std():+.5f}")

    if signed:
        bmin = -(1 << (nbits - 1))
        bmax = (1 << (nbits - 1)) - 1
        xmin = xoffset - (1 << nbits) * xscale / 2 - xscale / 2
        xmax = xoffset + (1 << nbits) * xscale / 2 - xscale / 2
    else:
        bmin = 0
        bmax = (1 << nbits) - 1
        xmin = xoffset - xscale / 2
        xmax = xoffset + (1 << nbits) * xscale - xscale / 2

    # Count clipped values
    c = (x < xmin) | (x >= xmax)
    print(f"{np.sum(c):d} samples clipped ({100 * np.sum(c) / n:.3f} percent)")

    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 8))

    ax1.hist(x, bins=200, range=[x.min(), x.max()], label="Input", histtype="step", density=True)
    ax1.hist(x_float, bins=257, range=[xmin, xmax], label="Reconstructed", histtype="step", density=True)
    ax1.legend()
    ax1.set_xlabel("Floating point values")
    ax1.set_ylabel("Fraction")
    ax1.axvline(xmin, color="C2")
    ax1.axvline(xmax, color="C2")
    
    ax2.hist(x_int, bins=2**8, range=[bmin, bmax], histtype="step", density=True)
    ax2.set_xlabel("Integer values")
    ax2.set_ylabel("Fraction")

    ax3.plot(t, xdiff)
    ax3.plot(t[c], xdiff[c], marker=".", linestyle="None", label="Clipped values")
    ax3.legend()
    ax3.set_xlabel("Time (s)")
    ax3.set_ylabel("Difference")
    ax3.axhline(xscale / 2, color="C2")
    ax3.axhline(-xscale / 2, color="C2")
    
    plt.tight_layout()
    plt.savefig("digitize.png", bbox_inches="tight")