"""Biquad filter design and processing.

Design functions implement the Audio EQ Cookbook formulas
(Bristow-Johnson, 2005) for peaking and shelving EQ filters.
Processing uses Direct Form II Transposed for numerical stability.
"""

import numpy as np


def design_peaking_eq(fc, gain_db, Q, fs):
    """Design a peaking (bell) EQ biquad section.

    Parameters
    ----------
    fc : float
        Center frequency in Hz.
    gain_db : float
        Gain at center frequency in dB.
    Q : float
        Quality factor (higher = narrower bandwidth).
    fs : float
        Sample rate in Hz.

    Returns
    -------
    sos : ndarray, shape (1, 6)
        Second-order section ``[b0, b1, b2, 1, a1, a2]``.
    """
    A = 10 ** (gain_db / 40)
    w0 = 2 * np.pi * fc / fs
    alpha = np.sin(w0) / (2 * Q)

    b0 = 1 + alpha * A
    b1 = -2 * np.cos(w0)
    b2 = 1 - alpha * A
    a0 = 1 + alpha / A
    a1 = -2 * np.cos(w0)
    a2 = 1 - alpha / A

    return np.array([[b0 / a0, b1 / a0, b2 / a0, 1.0, a1 / a0, a2 / a0]])


def design_low_shelf(fc, gain_db, Q, fs):
    """Design a low shelf biquad section.

    Parameters
    ----------
    fc : float
        Corner frequency in Hz.
    gain_db : float
        Shelf gain in dB (positive = boost, negative = cut).
    Q : float
        Quality factor controlling the transition slope.
    fs : float
        Sample rate in Hz.

    Returns
    -------
    sos : ndarray, shape (1, 6)
        Second-order section ``[b0, b1, b2, 1, a1, a2]``.
    """
    A = 10 ** (gain_db / 40)
    w0 = 2 * np.pi * fc / fs
    alpha = np.sin(w0) / (2 * Q)
    two_sqrt_A_alpha = 2 * np.sqrt(A) * alpha

    b0 = A * ((A + 1) - (A - 1) * np.cos(w0) + two_sqrt_A_alpha)
    b1 = 2 * A * ((A - 1) - (A + 1) * np.cos(w0))
    b2 = A * ((A + 1) - (A - 1) * np.cos(w0) - two_sqrt_A_alpha)
    a0 = (A + 1) + (A - 1) * np.cos(w0) + two_sqrt_A_alpha
    a1 = -2 * ((A - 1) + (A + 1) * np.cos(w0))
    a2 = (A + 1) + (A - 1) * np.cos(w0) - two_sqrt_A_alpha

    return np.array([[b0 / a0, b1 / a0, b2 / a0, 1.0, a1 / a0, a2 / a0]])


def design_high_shelf(fc, gain_db, Q, fs):
    """Design a high shelf biquad section.

    Parameters
    ----------
    fc : float
        Corner frequency in Hz.
    gain_db : float
        Shelf gain in dB (positive = boost, negative = cut).
    Q : float
        Quality factor controlling the transition slope.
    fs : float
        Sample rate in Hz.

    Returns
    -------
    sos : ndarray, shape (1, 6)
        Second-order section ``[b0, b1, b2, 1, a1, a2]``.
    """
    A = 10 ** (gain_db / 40)
    w0 = 2 * np.pi * fc / fs
    alpha = np.sin(w0) / (2 * Q)
    two_sqrt_A_alpha = 2 * np.sqrt(A) * alpha

    b0 = A * ((A + 1) + (A - 1) * np.cos(w0) + two_sqrt_A_alpha)
    b1 = -2 * A * ((A - 1) + (A + 1) * np.cos(w0))
    b2 = A * ((A + 1) + (A - 1) * np.cos(w0) - two_sqrt_A_alpha)
    a0 = (A + 1) - (A - 1) * np.cos(w0) + two_sqrt_A_alpha
    a1 = 2 * ((A - 1) - (A + 1) * np.cos(w0))
    a2 = (A + 1) - (A - 1) * np.cos(w0) - two_sqrt_A_alpha

    return np.array([[b0 / a0, b1 / a0, b2 / a0, 1.0, a1 / a0, a2 / a0]])


def design_parametric_eq(bands, fs):
    """Design a multi-band parametric equalizer as cascaded biquads.

    Parameters
    ----------
    bands : list of tuples
        Each tuple is ``(fc, gain_db, Q, type)`` where *type* is one of
        ``'peak'``, ``'lowshelf'``, or ``'highshelf'``.
    fs : float
        Sample rate in Hz.

    Returns
    -------
    sos : ndarray, shape (N, 6)
        Cascaded SOS matrix with one row per band.

    Raises
    ------
    ValueError
        If an unknown band type is specified.
    """
    designers = {
        "peak": design_peaking_eq,
        "lowshelf": design_low_shelf,
        "highshelf": design_high_shelf,
    }
    sections = []
    for fc, gain_db, Q, btype in bands:
        if btype not in designers:
            raise ValueError(
                f"Unknown band type '{btype}'. "
                f"Expected one of {list(designers.keys())}."
            )
        sections.append(designers[btype](fc, gain_db, Q, fs))
    return np.vstack(sections)


def quantize_sos(sos, bits):
    """Quantize SOS coefficients to fixed-point representation.

    Simulates the effect of representing biquad coefficients in
    Q(bits-1) fixed-point format. Coefficients are rounded to the
    nearest representable value and clipped to [-1, max_val] where
    ``max_val = (2^(bits-1) - 1) / 2^(bits-1)``.

    The ``a[0] = 1`` column (index 3) is left unmodified.

    Parameters
    ----------
    sos : ndarray, shape (N, 6)
        SOS coefficient matrix in ``[b0, b1, b2, 1, a1, a2]`` format.
    bits : int
        Word length (e.g. 16 for Q15, 32 for Q31).

    Returns
    -------
    sos_q : ndarray, shape (N, 6)
        Quantized SOS matrix.
    """
    scale = 2 ** (bits - 1)
    max_val = (scale - 1) / scale
    sos_q = sos.copy()
    for col in [0, 1, 2, 4, 5]:
        sos_q[:, col] = np.round(sos[:, col] * scale) / scale
        sos_q[:, col] = np.clip(sos_q[:, col], -1.0, max_val)
    return sos_q


def biquad_process(sos, x):
    """Process a signal through cascaded biquads (DF-II Transposed).

    Pure Python/NumPy implementation useful for understanding the
    algorithm and for porting to C or fixed-point platforms.

    Implements per section::

        y[n] = b0 * x[n] + s1
        s1   = b1 * x[n] - a1 * y[n] + s2
        s2   = b2 * x[n] - a2 * y[n]

    Parameters
    ----------
    sos : ndarray, shape (N, 6)
        SOS coefficient matrix in ``[b0, b1, b2, 1, a1, a2]`` format.
    x : ndarray, shape (M,)
        Input signal.

    Returns
    -------
    y : ndarray, shape (M,)
        Filtered output signal.
    """
    x = np.asarray(x, dtype=np.float64)
    y = x.copy()

    for section in sos:
        b0, b1, b2 = section[0], section[1], section[2]
        a1, a2 = section[4], section[5]
        s1 = 0.0
        s2 = 0.0
        for n in range(len(y)):
            xn = y[n]
            yn = b0 * xn + s1
            s1 = b1 * xn - a1 * yn + s2
            s2 = b2 * xn - a2 * yn
            y[n] = yn

    return y