computing odd Ap using bessel function

quantization_noise.py

"""
Quantization noise of a sinusoid: harmonic coefficients A_p, noise power,
and the impact of the bit count n.

Model: uniform quantizer, step Delta, input x(t) = A*sin(w t).
The quantization error expands into ODD harmonics only:

    e(t) = sum_{p odd} A_p sin(p w t)

For a mid-tread quantizer  y = Delta*round(x/Delta)  the error e = y - x has

    A_p = (2/pi) * sum_{m=1..inf} ((-1)^m / m) * J_p(2*pi*m*A/Delta)      (p odd)
    A_p = 0                                                              (p even)

The lecture-slide form drops the (-1)^m (a different quantizer phase). It gives
different individual line values but the SAME envelope, the same roll-off cliff
at p ~= 2*pi*A, and the same total power -> set convention='slide' to reproduce it.
The slide also writes e_n as the quantizer OUTPUT, adding a signal term
delta_{p1}*A at p=1; use Ap_full() for that.

For an n-bit full-scale sinusoid, A ~= 2^(n-1) LSB (with Delta = 1).

Reproduced + verified here:
  * even harmonics vanish (Jacobi-Anger has no J_0 -> no DC, mean = 0)
  * mean-square error -> Delta^2/12
  * SQNR = 6.02 n + 1.76 dB
  * roll-off cliff at p ~= 2*pi*A = pi*2^n            (doubles every bit)
  * each spur ~ -3 dB/bit in LSB, ~ -9 dB/bit relative to the signal
  * the analytic A_p matches a direct quantizer to 4+ digits

Requires: numpy, scipy, matplotlib.
"""

import numpy as np
from scipy.special import jv
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator


# --------------------------------------------------------------------------
# 1. Harmonic coefficient A_p
# --------------------------------------------------------------------------
def Ap_dist(p, A, Delta=1.0, convention='mid_tread', M=40000):
    """
    Distortion part of the p-th harmonic amplitude (excludes the signal term).
    `p` may be a scalar or 1-D array. `convention`:
        'mid_tread' -> (-1)^m sign, matches y = Delta*round(x/Delta)   [default]
        'slide'     -> all-positive sign, the lecture-slide expression
    The m-series tail decays ~ m^(-3/2); M~4e4 gives ~3-4 digit accuracy.
    """
    p = np.atleast_1d(np.asarray(p, dtype=float))
    m = np.arange(1, M + 1)
    sign = (-1.0) ** m if convention == 'mid_tread' else np.ones_like(m, float)
    coef = sign * 2.0 / (m * np.pi)                 # (M,)
    arg = 2 * np.pi * A * m / Delta                 # (M,)
    out = (jv(p[:, None], arg[None, :]) * coef[None, :]).sum(axis=1)
    out[np.mod(p, 2) == 0] = 0.0                     # even harmonics are exactly 0
    return out


def Ap_full(p, A, Delta=1.0, **kw):
    """Quantizer-output coefficient: distortion part + signal term A at p=1."""
    p = np.atleast_1d(np.asarray(p, dtype=float))
    out = Ap_dist(p, A, Delta, **kw)
    out[p == 1] += A
    return out


# --------------------------------------------------------------------------
# 2. Validation: project a real quantizer's error onto each harmonic
#    (single finely-sampled period -> robust against aliasing)
# --------------------------------------------------------------------------
def verify_against_quantizer(A, Delta=1.0, harmonics=(1, 3, 5, 7, 21, 41), N=4_000_000):
    th = 2 * np.pi * np.arange(N) / N
    x = A * np.sin(th)
    e = Delta * np.round(x / Delta) - x                       # mid-tread error
    print(f"  A={A} LSB :   direct projection   vs   analytic A_p")
    for p in harmonics:
        b = (2.0 / N) * np.sum(e * np.sin(p * th))            # Fourier sine coeff
        print(f"    p={p:>2}:  {b:+.5e}     {Ap_dist(p, A, Delta)[0]:+.5e}")


# --------------------------------------------------------------------------
# 3. Moments:  mean -> 0,  mean-square -> Delta^2/12
# --------------------------------------------------------------------------
def verify_moments(A, Delta=1.0, N=1 << 20):
    th = 2 * np.pi * np.arange(N) / N
    e = Delta * np.round(A * np.sin(th) / Delta) - A * np.sin(th)
    print(f"  A={A:>5} LSB :  mean(e) = {e.mean(): .2e}  (->0) | "
          f"mean(e^2) = {np.mean(e**2):.5f}  (-> Delta^2/12 = {Delta**2/12:.5f})")


# --------------------------------------------------------------------------
# 4. Per-bit scaling table (rms spur in LSB and relative to signal)
# --------------------------------------------------------------------------
def scaling_table(bits=range(3, 9)):
    print(f"{'n':>2} {'A':>5} {'cliff~2piA':>11} {'rms(LSB)':>12} "
          f"{'rms/A(dBc)':>11} {'dB/bit LSB':>11} {'dB/bit dBc':>11}")
    p_lsb = p_dbc = None
    for nb in bits:
        A = 2 ** (nb - 1)
        odd = np.arange(1, int(2 * np.pi * A), 2)
        rms = np.sqrt(np.mean(np.abs(Ap_dist(odd, A, M=12000)) ** 2))
        dbc = 20 * np.log10(rms / A)
        d1 = (20*np.log10(rms) - p_lsb) if p_lsb is not None else np.nan
        d2 = (dbc - p_dbc) if p_dbc is not None else np.nan
        print(f"{nb:>2} {A:>5} {int(2*np.pi*A):>11} {rms:>12.4e} "
              f"{dbc:>11.2f} {d1:>11.2f} {d2:>11.2f}")
        p_lsb, p_dbc = 20*np.log10(rms), dbc


# --------------------------------------------------------------------------
# 5. Plots
# --------------------------------------------------------------------------
def _style():
    plt.rcParams.update({'font.size': 11, 'font.family': 'DejaVu Sans',
                         'axes.spines.top': False, 'axes.spines.right': False})


def plot_Ap_vs_p(amps=(2, 8, 32), pmax=260, convention='mid_tread', fname='Ap_vs_p.png'):
    _style()
    colors = ['#2563eb', '#d97706', '#16a34a']
    odd = np.arange(1, pmax + 1, 2)
    fig, ax = plt.subplots(figsize=(10.5, 6.2), dpi=150)
    for A, c in zip(amps, colors):
        nb = int(round(np.log2(A))) + 1
        v = np.abs(Ap_dist(odd, A, convention=convention))
        ax.semilogy(odd, v, '-o', color=c, lw=1.0, ms=4.0, alpha=0.8,
                    label=f'$A={A}$ ($n\\approx{nb}$), cliff $2\\pi A\\approx{2*np.pi*A:.0f}$')
        ax.axvline(2 * np.pi * A, color=c, ls=':', lw=1.2, alpha=0.8)
    ax.set_xlim(0, pmax); ax.set_ylim(2e-4, 0.8)
    ax.xaxis.set_major_locator(MultipleLocator(20))
    ax.set_xlabel('harmonic index  p'); ax.set_ylabel(r'$|A_p|$  (LSB, $\Delta=1$)')
    ax.set_title(f'$|A_p|$ vs $p$  (odd $p$ only; even $p\\equiv0$) — {convention}')
    ax.legend(loc='upper right', frameon=False, fontsize=9)
    ax.grid(True, which='both', axis='y', alpha=0.18)
    fig.tight_layout(); fig.savefig(fname, dpi=150, bbox_inches='tight')
    print(f"  wrote {fname}")


def plot_n_impact(fname='n_impact.png'):
    _style()
    fig, (axL, axR) = plt.subplots(1, 2, figsize=(13.5, 5.6), dpi=150)
    n = np.arange(1, 17)
    axL.plot(n, 6.02 * n + 1.76, '-o', color='#2563eb', ms=5)
    axL.set_xlabel('number of bits  n'); axL.set_ylabel('SQNR (dB)')
    axL.set_title('SQNR $=6.02\\,n+1.76$ dB  (~6 dB / bit)')
    axL.grid(alpha=0.2); axL.xaxis.set_major_locator(MultipleLocator(2))
    for nb, c in {4: '#d97706', 6: '#16a34a', 8: '#9333ea'}.items():
        A = 2 ** (nb - 1)
        odd = np.arange(1, 1101, 2)
        dbc = 20 * np.log10(np.clip(np.abs(Ap_dist(odd, A, M=12000)) / A, 1e-12, None))
        axR.plot(odd, dbc, '-', color=c, lw=1.0, alpha=0.85,
                 label=f'$n={nb}$ ($A={A}$), cliff $\\approx{2*np.pi*A:.0f}$')
        axR.axvline(2 * np.pi * A, color=c, ls=':', lw=1.2, alpha=0.7)
    axR.set_xlim(0, 1100); axR.set_ylim(-120, 0)
    axR.set_xlabel('harmonic index  p (odd)'); axR.set_ylabel('spur rel. to signal (dBc)')
    axR.set_title('per-harmonic $A_p$: floor $\\approx-9$ dB/bit, cliff $\\propto 2^{\\,n}$')
    axR.grid(alpha=0.2); axR.legend(loc='lower left', fontsize=8.5, frameon=False)
    fig.tight_layout(); fig.savefig(fname, dpi=150, bbox_inches='tight')
    print(f"  wrote {fname}")


# --------------------------------------------------------------------------
if __name__ == '__main__':
    print("Analytic A_p vs a real round() quantizer (mid-tread):")
    verify_against_quantizer(8)
    print("\nMoments:")
    for A in (8, 32):
        verify_moments(A)
    print("\nPer-bit scaling:")
    scaling_table()
    print("\nFigures:")
    plot_Ap_vs_p()        # add convention='slide' to reproduce the lecture form
    plot_n_impact()
    print("\nDone.")