Signal Processing

General-purpose DSP utilities: filtering, windowing, peak finding, dithering, Teager energy, zero-crossings, and resampling.

Filters and windows

v_windows

V_WINDOWS - Generate a standard windowing function.

v_windows

v_windows(
    wtype, n=256, mode=None, p=None, ov=None
) -> ndarray

Generate a standard windowing function.

Parameters:

Name	Type	Description	Default
wtype	str or int	Window type name or code.	required
n	int	Number of output points. Default 256.	256
mode	str	Options string controlling scaling and sampling.	None
p	float or array_like	Parameter(s) for parameterized windows.	None
ov	int	Overlap for convolution with rectangle ('o' option).	None

Returns:

Name	Type	Description
w	ndarray	Window values (1D array of length n).

Source code in pyvoicebox/v_windows.py

def v_windows(wtype, n=256, mode=None, p=None, ov=None) -> np.ndarray:
    """Generate a standard windowing function.

    Parameters
    ----------
    wtype : str or int
        Window type name or code.
    n : int, optional
        Number of output points. Default 256.
    mode : str, optional
        Options string controlling scaling and sampling.
    p : float or array_like, optional
        Parameter(s) for parameterized windows.
    ov : int, optional
        Overlap for convolution with rectangle ('o' option).

    Returns
    -------
    w : ndarray
        Window values (1D array of length n).
    """
    # Resolve window type
    if isinstance(wtype, str):
        wtype_int = _WNAM.get(wtype.lower())
        if wtype_int is None:
            raise ValueError(f'Unknown window type: {wtype}')
    else:
        wtype_int = int(wtype)

    n = int(np.floor(n))

    if mode is None:
        mode = 'uw'

    # Parse mode flags
    mm = np.zeros(len(mode), dtype=int)
    ll = 'hc lrbns'
    for i, c in enumerate(ll):
        if c != ' ':
            for j, mc in enumerate(mode):
                if mc == c:
                    mm[j] = i - 2  # offset so h=-2, c=-1, l=1, r=2, b=3, n=4, s=5

    max_mm = max(mm) if len(mm) > 0 else 0
    min_mm = min(mm) if len(mm) > 0 else 0

    # k index into kk table (0-based)
    k = 3 * max(max_mm, 0) + max(-min_mm, 0)
    if k < 3:
        if wtype_int in _ZERO_END_WINDOWS:
            k += 12  # default to 'n' option

    kk_row = _KK[k]

    # Handle 'o' option (convolve with rectangle)
    do_conv = 'o' in mode
    if do_conv:
        if ov is None:
            ov = n // 2
        n = n - ov + 1
    else:
        ov = 0

    fn = int(np.floor(n))
    kp = kk_row[2] * n + kk_row[3]
    ks = kk_row[0] * fn + kk_row[1]
    v = ((np.arange(0, 2 * fn, 2) + ks) / kp)

    # Generate window
    if wtype_int == 1:  # rectangle
        w = np.ones_like(v)
    elif wtype_int == 2:  # hanning
        w = 0.5 + 0.5 * np.cos(np.pi * v)
    elif wtype_int == 3:  # hamming
        w = 0.54 + 0.46 * np.cos(np.pi * v)
    elif wtype_int == 4:  # harris3
        w = 0.42323 + 0.49755 * np.cos(np.pi * v) + 0.07922 * np.cos(2 * np.pi * v)
    elif wtype_int == 5:  # harris4
        w = (0.35875 + 0.48829 * np.cos(np.pi * v) +
             0.14128 * np.cos(2 * np.pi * v) + 0.01168 * np.cos(3 * np.pi * v))
    elif wtype_int == 6:  # blackman
        w = 0.42 + 0.5 * np.cos(np.pi * v) + 0.08 * np.cos(2 * np.pi * v)
    elif wtype_int == 7:  # vorbis
        w = np.sin(0.25 * np.pi * (1 + np.cos(np.pi * v)))
    elif wtype_int == 8:  # rsqvorbis
        w = 0.571 - 0.429 * np.cos(0.5 * np.pi * (1 + np.cos(np.pi * v)))
    elif wtype_int == 9:  # triangle
        pp = p if p is not None else 1
        if not np.isscalar(pp):
            pp = pp[0]
        w = 1 - np.abs(v) ** pp
    elif wtype_int == 10:  # cos
        pp = p if p is not None else 1
        if not np.isscalar(pp):
            pp = pp[0]
        w = np.cos(0.5 * np.pi * v) ** pp
    elif wtype_int == 11:  # kaiser
        pp = p if p is not None else 8
        if not np.isscalar(pp):
            pp = pp[0]
        w = besseli0(pp * np.sqrt(np.maximum(1 - v ** 2, 0))) / besseli0(pp)
    elif wtype_int == 12:  # gaussian
        pp = p if p is not None else 3
        if not np.isscalar(pp):
            pp = pp[0]
        w = np.exp(-0.5 * pp ** 2 * v ** 2)
    elif wtype_int == 13:  # cauchy
        pp = p if p is not None else 1
        if not np.isscalar(pp):
            pp = pp[0]
        w = (1 + (pp * v) ** 2) ** -1
    elif wtype_int == 14:  # dolph
        raise NotImplementedError('Dolph-Chebyshev window not yet implemented')
    elif wtype_int == 15:  # tukey
        pp = p if p is not None else 0.5
        if not np.isscalar(pp):
            pp = pp[0]
        if pp > 0:
            pp = min(pp, 1.0)
            w = 0.5 + 0.5 * np.cos(np.pi * np.maximum(1 + (np.abs(v) - 1) / pp, 0))
        else:
            w = np.ones_like(v)
    else:
        raise ValueError(f'Unknown window type code: {wtype_int}')

    # Convolve with rectangle
    if do_conv and ov > 0:
        w = np.cumsum(w)
        orig_n = n
        w_ext = np.zeros(orig_n + ov - 1)
        w_ext[:orig_n] = w
        w_ext[orig_n:orig_n + ov - 1] = w[orig_n - 1] - w[orig_n - ov:orig_n - 1]
        w_ext[ov:orig_n] = w_ext[ov:orig_n] - w[:orig_n - ov]
        w = w_ext
        n = orig_n + ov - 1

    # Scale
    g = 1.0
    for c in mode:
        if '1' <= c <= '9':
            g = 1.0 / (ord(c) - ord('0'))
            break

    if 'd' in mode:
        w = w * (g / np.sum(w))
    elif 'D' in mode or 'a' in mode:
        w = w * (g / np.mean(w))
    elif 'e' in mode.replace('ne', '').replace('se', ''):
        # 'e' for energy, but not 'ne' or 'se' mode letters
        if any(c == 'e' for c in mode if c not in 'ns'):
            w = w * np.sqrt(g / np.sum(w ** 2))
    elif 'E' in mode:
        w = w * np.sqrt(g / np.mean(w ** 2))
    elif 'p' in mode.replace('sp', ''):
        w = w * (g / np.max(w))

    if 'q' in mode:
        w = np.sqrt(w)

    return w

v_windinfo

V_WINDINFO - Window information and figures of merit.

v_windinfo

v_windinfo(w, fs=1) -> dict

Calculate window information and figures of merit.

Parameters:

Name	Type	Description	Default
w	array_like	Window vector.	required
fs	float	Sampling frequency. Default 1.	1

Returns:

Name	Type	Description
info	dict	Dictionary with window properties: - nw: window length - ewgdelay: energy centroid delay from first sample - dcgain: DC gain in dB - enbw: equivalent noise bandwidth - scallop: scalloping loss in dB

Source code in pyvoicebox/v_windinfo.py

def v_windinfo(w, fs=1) -> dict:
    """Calculate window information and figures of merit.

    Parameters
    ----------
    w : array_like
        Window vector.
    fs : float, optional
        Sampling frequency. Default 1.

    Returns
    -------
    info : dict
        Dictionary with window properties:
        - nw: window length
        - ewgdelay: energy centroid delay from first sample
        - dcgain: DC gain in dB
        - enbw: equivalent noise bandwidth
        - scallop: scalloping loss in dB
    """
    w = np.asarray(w, dtype=float).ravel()
    nw = len(w)

    # DC gain
    dc = np.sum(w)
    dcgain = 20 * np.log10(abs(dc)) if dc != 0 else -np.inf

    # Energy centroid delay
    energy = np.sum(w ** 2 * np.arange(nw))
    total_energy = np.sum(w ** 2)
    ewgdelay = energy / total_energy if total_energy > 0 else 0

    # Equivalent noise bandwidth (normalized)
    enbw = nw * total_energy / (dc ** 2) if dc != 0 else np.inf

    # Scalloping loss
    nfft = max(1024, 4 * nw)
    W = np.fft.fft(w, nfft)
    W_half = np.abs(W[:nfft // 2 + 1])
    if W_half[0] > 0:
        # Response at half-bin spacing
        w_shifted = w * np.exp(1j * np.pi * np.arange(nw) / nw)
        W_shift = np.abs(np.sum(w_shifted))
        scallop = 20 * np.log10(W_shift / abs(dc)) if abs(dc) > 0 else -np.inf
    else:
        scallop = 0

    info = {
        'nw': nw,
        'ewgdelay': ewgdelay,
        'dcgain': dcgain,
        'enbw': enbw,
        'scallop': scallop,
    }
    return info

v_filterbank

V_FILTERBANK - Apply a bank of IIR filters to a signal.

v_filterbank

v_filterbank(b, a, x, zi=None) -> tuple[ndarray, ndarray]

Apply a bank of filters to a signal.

Parameters:

Name	Type	Description	Default
b	list of array_like or ndarray	Numerator coefficients. If 2D, each row is a filter.	required
a	list of array_like or ndarray	Denominator coefficients. If 2D, each row is a filter.	required
x	array_like	Input signal.	required
zi	list of array_like	Initial filter states.	None

Returns:

Name	Type	Description
y	ndarray	Output signals, one column per filter.
zf	list of ndarray	Final filter states.

Source code in pyvoicebox/v_filterbank.py

def v_filterbank(b, a, x, zi=None) -> tuple[np.ndarray, np.ndarray]:
    """Apply a bank of filters to a signal.

    Parameters
    ----------
    b : list of array_like or ndarray
        Numerator coefficients. If 2D, each row is a filter.
    a : list of array_like or ndarray
        Denominator coefficients. If 2D, each row is a filter.
    x : array_like
        Input signal.
    zi : list of array_like, optional
        Initial filter states.

    Returns
    -------
    y : ndarray
        Output signals, one column per filter.
    zf : list of ndarray
        Final filter states.
    """
    x = np.asarray(x, dtype=float).ravel()

    if isinstance(b, np.ndarray) and b.ndim == 2:
        nf = b.shape[0]
        b_list = [b[i, :] for i in range(nf)]
        a_list = [a[i, :] for i in range(nf)]
    elif isinstance(b, list):
        nf = len(b)
        b_list = b
        a_list = a
    else:
        nf = 1
        b_list = [np.asarray(b, dtype=float).ravel()]
        a_list = [np.asarray(a, dtype=float).ravel()]

    y = np.zeros((len(x), nf))
    zf = []

    for i in range(nf):
        bi = np.asarray(b_list[i], dtype=float).ravel()
        ai = np.asarray(a_list[i], dtype=float).ravel()
        if zi is not None and i < len(zi):
            yi, zfi = lfilter(bi, ai, x, zi=zi[i])
        else:
            yi = lfilter(bi, ai, x)
            zfi = np.zeros(max(len(bi), len(ai)) - 1)
        y[:, i] = yi
        zf.append(zfi)

    return y, zf

v_maxfilt

V_MAXFILT - Find max of an exponentially weighted sliding window.

v_maxfilt

v_maxfilt(
    x, f=1, n=None, d=None, x0=None
) -> tuple[ndarray, ndarray, ndarray]

Find max of an exponentially weighted sliding window.

Calculates y(p) = max(f^r * x(p-r), r=0:n-1) where x®=-inf for r<0.

Parameters:

Name	Type	Description	Default
x	array_like	Input data (vector or matrix).	required
f	float	Exponential forgetting factor, range [0, 1]. Default 1 (no forgetting).	1
n	int or None	Length of sliding window. Default is inf.	None
d	int or None	Dimension to work along (0-based). Default is first non-singleton.	None
x0	array_like or None	Initial values placed in front of x data.	None

Returns:

Name	Type	Description
y	ndarray	Output array, same shape as x.
k	ndarray	Index array such that y = x[k] (0-based).
y0	ndarray	Last n-1 values for subsequent calls (or last output if n=inf).

Source code in pyvoicebox/v_maxfilt.py

def v_maxfilt(x, f=1, n=None, d=None, x0=None) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Find max of an exponentially weighted sliding window.

    Calculates y(p) = max(f^r * x(p-r), r=0:n-1) where x(r)=-inf for r<0.

    Parameters
    ----------
    x : array_like
        Input data (vector or matrix).
    f : float, optional
        Exponential forgetting factor, range [0, 1]. Default 1 (no forgetting).
    n : int or None, optional
        Length of sliding window. Default is inf.
    d : int or None, optional
        Dimension to work along (0-based). Default is first non-singleton.
    x0 : array_like or None, optional
        Initial values placed in front of x data.

    Returns
    -------
    y : ndarray
        Output array, same shape as x.
    k : ndarray
        Index array such that y = x[k] (0-based).
    y0 : ndarray
        Last n-1 values for subsequent calls (or last output if n=inf).
    """
    x = np.asarray(x, dtype=float)
    s = x.shape

    if d is None:
        d_found = None
        for i, si in enumerate(s):
            if si > 1:
                d_found = i
                break
        d = d_found if d_found is not None else 0

    # Concatenate x0 if provided
    if x0 is not None:
        x0 = np.asarray(x0, dtype=float)
        if x0.size > 0:
            y = np.concatenate([x0, x], axis=d)
            nx0 = x0.shape[d]
        else:
            y = x.copy()
            nx0 = 0
    else:
        y = x.copy()
        nx0 = 0

    # Move working dimension to axis 0
    y = np.moveaxis(y, d, 0)
    s_y = y.shape
    s1 = s_y[0]

    if n is None:
        n0 = np.inf
    else:
        n0 = max(n, 1)

    nn = n0
    nn = int(min(nn, s1)) if not np.isinf(nn) else s1

    # Reshape to 2D: (s1, rest)
    rest_shape = s_y[1:]
    rest_size = int(np.prod(rest_shape)) if len(rest_shape) > 0 else 1
    y_2d = y.reshape(s1, rest_size)
    k_2d = np.tile(np.arange(s1)[:, np.newaxis], (1, rest_size))

    # The MATLAB algorithm: doubling step max-propagation
    if nn > 1:
        j = 1
        j2 = 1
        while j > 0:
            g = f ** j
            # Compare y[j:] with g*y[:-j]
            for col in range(rest_size):
                mask = y_2d[j:, col] <= g * y_2d[:s1 - j, col]
                rows = np.where(mask)[0]
                y_2d[rows + j, col] = g * y_2d[rows, col]
                k_2d[rows + j, col] = k_2d[rows, col]
            j2 = j2 + j
            nj = nn - j2
            if nj <= 0:
                j = 0
            else:
                j = min(j2, nj)

    # Compute y0 output
    if not np.isinf(n0):
        ny0 = min(s1, nn - 1)
    else:
        ny0 = min(s1, 1)

    if ny0 <= 0 or np.isinf(n0):
        # For n=inf, y0 is the last output
        if np.isinf(n0) and ny0 == 1:
            y0_2d = y_2d[-1:, :]
        else:
            y0_2d = np.zeros((max(ny0, 0), rest_size))
    else:
        y0_2d = y_2d[s1 - ny0:, :]

    # Remove prepended x0 data
    if nx0 > 0:
        y_2d = y_2d[nx0:, :]
        k_2d = k_2d[nx0:, :] - nx0
        s1_out = s1 - nx0
    else:
        s1_out = s1

    # Reshape back
    out_shape = (s1_out,) + rest_shape
    y_out = y_2d.reshape(out_shape)
    k_out = k_2d.reshape(out_shape)

    # Move axis back
    y_out = np.moveaxis(y_out, 0, d)
    k_out = np.moveaxis(k_out, 0, d)

    # Reshape y0
    y0_shape = list(s)
    y0_shape[d] = ny0
    y0_2d_reshaped = y0_2d.reshape((ny0,) + rest_shape)
    y0_out = np.moveaxis(y0_2d_reshaped, 0, d)

    return y_out, k_out, y0_out

v_momfilt

V_MOMFILT - Calculate moments of a signal using a sliding window.

v_momfilt

v_momfilt(x, r, w=None, m=None) -> tuple[ndarray, int]

Calculate moments of a signal using a sliding window.

Parameters:

Name	Type	Description	Default
x	array_like	Input signal.	required
r	array_like	List of moments to calculate (+ve relative to mean, -ve relative to zero).	required
w	int or array_like	Window or window length. Default: Hamming of length(x).	None
m	int	Center sample of window (1-based). Default: ceil((1+len(w))/2).	None

Returns:

Name	Type	Description
y	(ndarray, shape(len(x), len(r)))	Moments.
mm	int	Actual value of m used.

Source code in pyvoicebox/v_momfilt.py

def v_momfilt(x, r, w=None, m=None) -> tuple[np.ndarray, int]:
    """Calculate moments of a signal using a sliding window.

    Parameters
    ----------
    x : array_like
        Input signal.
    r : array_like
        List of moments to calculate (+ve relative to mean, -ve relative to zero).
    w : int or array_like, optional
        Window or window length. Default: Hamming of length(x).
    m : int, optional
        Center sample of window (1-based). Default: ceil((1+len(w))/2).

    Returns
    -------
    y : ndarray, shape (len(x), len(r))
        Moments.
    mm : int
        Actual value of m used.
    """
    x = np.asarray(x, dtype=float).ravel()
    r = np.round(np.asarray(r, dtype=float)).astype(int).ravel()

    if w is None:
        w = hamming(len(x), sym=True)
    elif np.isscalar(w):
        w = hamming(int(w), sym=True)
    else:
        w = np.asarray(w, dtype=float).ravel()

    lw = len(w)
    if m is None:
        m = (1 + lw) / 2.0
    m = max(int(round(m)), 1)
    mm = m

    lx = len(x)
    lxx = lx + m - 1
    xx = np.zeros(lxx)
    xx[:lx] = x

    # Compute y0 = filter(w, 1, ones)
    cw = np.cumsum(w)
    sw = cw[-1]
    y0 = np.full(lxx, sw)
    lxw = min(lxx, lw)
    y0[:lxw] = cw[:lxw]
    if m > 1:
        y0[lx:lx + m - 1] -= cw[:m - 1]
    yd = y0[m - 1:]
    yd[np.abs(yd) < np.finfo(float).eps] = 1.0

    nr = len(r)
    y = np.zeros((lx, nr))

    mr = max(np.abs(r)) if len(r) > 0 else 0
    mk = np.zeros(mr, dtype=int)
    neg_r = r[r < 0]
    if len(neg_r) > 0:
        mk[-neg_r - 1] = 1
    maxr = max(r) if len(r) > 0 else 0
    if maxr > 1:
        mk[:maxr] = 1

    ml = np.where(mk > 0)[0]  # 0-based indices of needed moments
    lml = len(ml)

    if lml > 0:
        moments_needed = ml + 1  # actual moment orders (1-based)
        xm = np.zeros((lxx, lml))
        for i, mom in enumerate(moments_needed):
            xm[:, i] = lfilter(w, [1.0], xx ** mom)
        xm = xm[m - 1:, :] / yd[:, np.newaxis]

        # Build mapping from moment order to column index
        mlx = np.zeros(mr, dtype=int)
        for i, mi in enumerate(ml):
            mlx[mi] = i

    # Fill in zero-centered moments (negative r)
    fr = np.where(r < 0)[0]
    if len(fr) > 0:
        for fi in fr:
            y[:, fi] = xm[:, mlx[-r[fi] - 1]]

    # 0th moment
    fr = np.where(r == 0)[0]
    if len(fr) > 0:
        y[:, fr] = 1.0

    # 1st moment about mean
    fr = np.where(r == 1)[0]
    if len(fr) > 0:
        y[:, fr] = 0.0

    # 2nd moment about mean (variance)
    fr = np.where(r == 2)[0]
    if len(fr) > 0:
        yfr = xm[:, mlx[1]] - xm[:, mlx[0]] ** 2
        for fi in fr:
            y[:, fi] = yfr

    # Higher moments about the mean
    if maxr > 2 and lml > 0:
        bc = np.array([1, -2, 1], dtype=float)
        mon = np.array([1, -1], dtype=float)
        am = np.zeros((lx, maxr))
        am[:, 0] = xm[:, mlx[0]]
        for mi in range(1, maxr):
            am[:, mi] = xm[:, mlx[0]] ** (mi + 1)

        for i in range(3, maxr + 1):
            bc = np.convolve(bc, mon)
            fr = np.where(r == i)[0]
            if len(fr) > 0:
                yfr = xm[:, mlx[i - 1]]
                for j in range(1, i - 1):
                    yfr = yfr + xm[:, mlx[i - j - 2]] * am[:, j - 1] * bc[j]
                yfr = yfr + am[:, i - 1] * (bc[i - 1] + bc[i])
                for fi in fr:
                    y[:, fi] = yfr

    return y, mm

v_meansqtf

V_MEANSQTF - Mean square transfer function of a filter.

v_meansqtf

v_meansqtf(b, a=None) -> float

Calculate the mean square transfer function of a filter.

This equals the average output power when the filter is fed with unit variance white noise.

Parameters:

Name	Type	Description	Default
b	array_like	Numerator filter coefficients.	required
a	array_like	Denominator filter coefficients. Default is [1].	None

Returns:

Name	Type	Description
d	float	Mean square transfer function.

Source code in pyvoicebox/v_meansqtf.py

def v_meansqtf(b, a=None) -> float:
    """Calculate the mean square transfer function of a filter.

    This equals the average output power when the filter is fed
    with unit variance white noise.

    Parameters
    ----------
    b : array_like
        Numerator filter coefficients.
    a : array_like, optional
        Denominator filter coefficients. Default is [1].

    Returns
    -------
    d : float
        Mean square transfer function.
    """
    b = np.asarray(b, dtype=float).ravel()
    if a is None or len(np.atleast_1d(a)) == 1:
        return float(b @ b)

    a = np.asarray(a, dtype=float).ravel()
    m = v_lpcar2ra(b.reshape(1, -1))
    m[0, 0] *= 0.5
    rr = v_lpcar2rr(a.reshape(1, -1), len(m[0]) - 1)
    d = 2.0 * (rr @ m.T).item()
    return d

v_resample

V_RESAMPLE - Resample and remove end transients.

v_resample

v_resample(x, p, q, n=10, b=5) -> ndarray

Resample signal and remove end transients.

Parameters:

Name	Type	Description	Default
x	array_like	Input signal.	required
p	int	Upsampling factor.	required
q	int	Downsampling factor.	required
n	int	Filter length. Default 10.	10
b	float	Kaiser window beta. Default 5.	5

Returns:

Name	Type	Description
y	ndarray	Resampled output signal.

Source code in pyvoicebox/v_resample.py

def v_resample(x, p, q, n=10, b=5) -> np.ndarray:
    """Resample signal and remove end transients.

    Parameters
    ----------
    x : array_like
        Input signal.
    p : int
        Upsampling factor.
    q : int
        Downsampling factor.
    n : int, optional
        Filter length. Default 10.
    b : float, optional
        Kaiser window beta. Default 5.

    Returns
    -------
    y : ndarray
        Resampled output signal.
    """
    x = np.asarray(x, dtype=float)
    if p == q:
        return x.copy()

    y = resample_poly(x, p, q)
    m = int(np.ceil(n * max(p / q, 1)))
    if x.ndim == 1:
        y = y[m:-m] if len(y) > 2 * m else y
    else:
        y = y[m:-m, :] if y.shape[0] > 2 * m else y
    return y

v_dlyapsq

V_DLYAPSQ - Solve discrete Lyapunov equation in square root form.

v_dlyapsq

v_dlyapsq(a, b) -> ndarray

Solve the discrete Lyapunov equation AV'VA' - V'V + BB' = 0.

V is upper triangular with real non-negative diagonal entries. Equivalent to chol(dlyap(a, b@b')) but better conditioned.

Parameters:

Name	Type	Description	Default
a	(array_like, shape(n, n))	State matrix.	required
b	(array_like, shape(n, m))	Input matrix.	required

Returns:

Name	Type	Description
v	(ndarray, shape(n, n))	Upper triangular solution.

Source code in pyvoicebox/v_dlyapsq.py

def v_dlyapsq(a, b) -> np.ndarray:
    """Solve the discrete Lyapunov equation AV'VA' - V'V + BB' = 0.

    V is upper triangular with real non-negative diagonal entries.
    Equivalent to chol(dlyap(a, b@b')) but better conditioned.

    Parameters
    ----------
    a : array_like, shape (n, n)
        State matrix.
    b : array_like, shape (n, m)
        Input matrix.

    Returns
    -------
    v : ndarray, shape (n, n)
        Upper triangular solution.
    """
    a = np.asarray(a, dtype=complex if np.iscomplexobj(a) else float)
    b = np.asarray(b, dtype=complex if np.iscomplexobj(b) else float)

    # Schur decomposition of a'
    s, q = schur(a.conj().T)
    s, q = rsf2csf(s, q)

    # QR factorization of b'*q
    r = np.linalg.qr(b.conj().T @ q, mode='r')

    m, n = r.shape if r.ndim == 2 else (1, r.shape[0])
    if r.ndim == 1:
        r = r.reshape(1, -1)

    u = np.zeros((n, n), dtype=complex)

    if m == 1:
        for i in range(n - 1):
            si = s[i, i]
            aa = np.sqrt(1 - si * np.conj(si))
            u[i, i] = r[0, 0] / aa
            rhs = u[i, i] * np.conj(si) * s[i, i + 1:] + aa * r[0, 1:]
            lhs = np.eye(n - i - 1) - np.conj(si) * s[i + 1:, i + 1:]
            u[i, i + 1:] = np.linalg.solve(lhs, rhs)
            r = np.atleast_2d(aa * (u[i, i] * s[i, i + 1:] + u[i, i + 1:] @ s[i + 1:, i + 1:]) - si * r[0, 1:])
        u[n - 1, n - 1] = r[0, 0] / np.sqrt(1 - s[n - 1, n - 1] * np.conj(s[n - 1, n - 1]))
    else:
        # General case with m > 1
        w = np.zeros(m)
        w[m - 1] = 1.0
        em = np.eye(m)
        for i in range(n - m):
            si = s[i, i]
            aa = np.sqrt(1 - si * np.conj(si))
            u[i, i] = r[0, 0] / aa
            rhs = u[i, i] * np.conj(si) * s[i, i + 1:] + aa * r[0, 1:]
            lhs = np.eye(n - i - 1) - np.conj(si) * s[i + 1:, i + 1:]
            u[i, i + 1:] = np.linalg.solve(lhs, rhs)
            vv = aa * (u[i, i] * s[i, i + 1:] + u[i, i + 1:] @ s[i + 1:, i + 1:]) - si * r[0, 1:]
            rr = np.zeros((m, n - i - 1), dtype=complex)
            rr[:m - 1, :] = r[1:, 1:]
            # QR update equivalent
            combined = np.vstack([rr, vv.reshape(1, -1)])
            qq, r_new = np.linalg.qr(combined, mode='reduced')
            r = r_new[:m, :]

        for i in range(max(0, n - m), n - 1):
            si = s[i, i]
            aa = np.sqrt(1 - si * np.conj(si))
            u[i, i] = r[0, 0] / aa
            rhs = u[i, i] * np.conj(si) * s[i, i + 1:] + aa * r[0, 1:]
            lhs = np.eye(n - i - 1) - np.conj(si) * s[i + 1:, i + 1:]
            u[i, i + 1:] = np.linalg.solve(lhs, rhs)
            vv = aa * (u[i, i] * s[i, i + 1:] + u[i, i + 1:] @ s[i + 1:, i + 1:]) - si * r[0, 1:]
            ni = n - i - 1
            rr = np.zeros((ni + 1, ni), dtype=complex)
            rr[:ni, :] = r[1:ni + 1, 1:] if r.shape[0] > 1 else np.zeros((ni, ni), dtype=complex)
            combined = np.vstack([rr[:ni, :], vv.reshape(1, -1)])
            qq, r_new = np.linalg.qr(combined, mode='reduced')
            r = r_new[:ni, :]

        u[n - 1, n - 1] = r.ravel()[0] / np.sqrt(1 - s[n - 1, n - 1] * np.conj(s[n - 1, n - 1]))

    v = np.triu(np.linalg.qr(u @ q.conj().T, mode='r'))
    dv = np.diag(v)
    ix = dv != 0
    signs = np.ones(n, dtype=complex)
    signs[ix] = np.abs(dv[ix]) / dv[ix]
    v = np.diag(signs) @ v

    if np.isrealobj(a) and np.isrealobj(b):
        v = np.real(v)
    return v

Analysis primitives

v_findpeaks

V_FINDPEAKS - Find peaks with optional quadratic interpolation.

v_findpeaks

v_findpeaks(
    y, m="", w=None, x=None
) -> tuple[ndarray, ndarray]

Find peaks with optional quadratic interpolation.

Parameters:

Name	Type	Description	Default
y	array_like	Input signal.	required
m	str	Mode string: 'f' - include first sample if downward initial slope 'l' - include last sample if upward final slope 'm' - return only the maximum peak 'q' - quadratic interpolation 'v' - find valleys instead of peaks	''
w	float	Width tolerance. Peaks closer than w will have the lower one removed.	None
x	array_like	X-axis values for y. Default: 0-based indices (1-based in MATLAB).	None

Returns:

Name	Type	Description
k	ndarray	Positions of peaks (1-based if x not given, fractional if 'q').
v	ndarray	Peak amplitudes.

Source code in pyvoicebox/v_findpeaks.py

def v_findpeaks(y, m='', w=None, x=None) -> tuple[np.ndarray, np.ndarray]:
    """Find peaks with optional quadratic interpolation.

    Parameters
    ----------
    y : array_like
        Input signal.
    m : str, optional
        Mode string:
          'f' - include first sample if downward initial slope
          'l' - include last sample if upward final slope
          'm' - return only the maximum peak
          'q' - quadratic interpolation
          'v' - find valleys instead of peaks
    w : float, optional
        Width tolerance. Peaks closer than w will have the lower one removed.
    x : array_like, optional
        X-axis values for y. Default: 0-based indices (1-based in MATLAB).

    Returns
    -------
    k : ndarray
        Positions of peaks (1-based if x not given, fractional if 'q').
    v : ndarray
        Peak amplitudes.
    """
    y = np.asarray(y, dtype=float).ravel()
    ny = len(y)

    if x is not None:
        x = np.asarray(x, dtype=float).ravel()

    if 'v' in m:
        y = -y.copy()
    else:
        y = y.copy()

    dx = y[1:] - y[:-1]
    r = np.where(dx > 0)[0]  # indices where signal rises
    f_idx = np.where(dx < 0)[0]  # indices where signal falls

    k = np.array([], dtype=float)
    v = np.array([], dtype=float)

    if len(r) > 0 and len(f_idx) > 0:
        # Convert to 1-based indices to match MATLAB logic throughout
        r1 = r + 1  # 1-based rise indices
        f1 = f_idx + 1  # 1-based fall indices

        # Compute rs: time since the last rise
        dr = r1.copy()
        dr[1:] = r1[1:] - r1[:-1]
        rc = np.ones(ny)
        rc[r1] = 1 - dr  # MATLAB: rc(r+1) where r is 1-based
        rc[0] = 0
        rs = np.cumsum(rc).astype(int)

        # Compute fs: time since the last fall
        df = f1.copy()
        df[1:] = f1[1:] - f1[:-1]
        fc = np.ones(ny)
        fc[f1] = 1 - df
        fc[0] = 0
        fs = np.cumsum(fc).astype(int)

        # Compute rq: time to the next rise
        # MATLAB: rp([1; r+1]) = [dr-1; ny-r(end)-1]
        # [1; r+1] in MATLAB (1-based) = [0; r1] in 0-based
        rp = -np.ones(ny)
        rp_indices = np.concatenate([[0], r1])
        rp_values = np.concatenate([dr - 1, [ny - r1[-1] - 1]])
        rp[rp_indices] = rp_values
        rq = np.cumsum(rp).astype(int)

        # Compute fq: time to the next fall
        fp = -np.ones(ny)
        fp_indices = np.concatenate([[0], f1])
        fp_values = np.concatenate([df - 1, [ny - f1[-1] - 1]])
        fp[fp_indices] = fp_values
        fq = np.cumsum(fp).astype(int)

        # Find peaks: (rs<fs) & (fq<rq) & centered in plateau
        k_idx = np.where(
            (rs < fs) & (fq < rq) &
            (np.floor((fq - rs) / 2.0) == 0)
        )[0]
        v = y[k_idx].copy()

        if 'q' in m:
            if x is not None:
                xm = x[k_idx - 1] - x[k_idx]
                xp = x[k_idx + 1] - x[k_idx]
                ym = y[k_idx - 1] - y[k_idx]
                yp = y[k_idx + 1] - y[k_idx]
                d_val = xm * xp * (xm - xp)
                b = 0.5 * (yp * xm ** 2 - ym * xp ** 2)
                a = xm * yp - xp * ym
                j = a > 0  # j=0 on a plateau
                v[j] = y[k_idx[j]] + b[j] ** 2 / (a[j] * d_val[j])
                k = np.empty_like(k_idx, dtype=float)
                k[j] = x[k_idx[j]] + b[j] / a[j]
                k[~j] = 0.5 * (x[k_idx[~j] + fq[k_idx[~j]]] +
                                x[k_idx[~j] - rs[k_idx[~j]]])
            else:
                # k_idx are 0-based; MATLAB uses 1-based indices
                k_1based = k_idx + 1  # convert to 1-based for computation
                b = 0.25 * (y[k_idx + 1] - y[k_idx - 1])
                a = y[k_idx] - 2 * b - y[k_idx - 1]
                j = a > 0
                k = k_1based.astype(float)
                v[j] = y[k_idx[j]] + b[j] ** 2 / a[j]
                k[j] = k_1based[j] + b[j] / a[j]
                k[~j] = k_1based[~j] + (fq[k_idx[~j]] - rs[k_idx[~j]]) / 2.0
        elif x is not None:
            k = x[k_idx].copy()
        else:
            k = (k_idx + 1).astype(float)  # 1-based

    # Add first and last samples if requested
    if ny > 1:
        if 'f' in m and y[0] > y[1]:
            v = np.concatenate([[y[0]], v])
            if x is not None:
                k = np.concatenate([[x[0]], k])
            else:
                k = np.concatenate([[1.0], k])

        if 'l' in m and y[ny - 2] < y[ny - 1]:
            v = np.concatenate([v, [y[ny - 1]]])
            if x is not None:
                k = np.concatenate([k, [x[ny - 1]]])
            else:
                k = np.concatenate([k, [float(ny)]])

        # Purge nearby peaks
        if 'm' in m:
            if len(v) > 0:
                iv = np.argmax(v)
                v = np.array([v[iv]])
                k = np.array([k[iv]])
        elif w is not None and np.isscalar(w) and w > 0:
            j = np.where(k[1:] - k[:-1] <= w)[0]
            while len(j) > 0:
                j = j + (v[j] >= v[j + 1]).astype(int)
                k = np.delete(k, j)
                v = np.delete(v, j)
                j = np.where(k[1:] - k[:-1] <= w)[0]
    elif ny > 0 and ('f' in m or 'l' in m):
        v = y.copy()
        if x is not None:
            k = x.copy()
        else:
            k = np.array([1.0])

    if 'v' in m:
        v = -v

    return k, v

v_zerocros

V_ZEROCROS - Find zero crossings in a signal.

v_zerocros

v_zerocros(y, m='b', x=None) -> tuple[ndarray, ndarray]

Find zero crossings in a signal.

Parameters:

Name	Type	Description	Default
y	array_like	Input waveform.	required
m	str	Mode string: 'p' - positive crossings only 'n' - negative crossings only 'b' - both (default) 'r' - round to sample values	'b'
x	array_like	X-axis values for y. Default: 1-based indices (MATLAB convention).	None

Returns:

Name	Type	Description
t	ndarray	X-axis positions of zero crossings (1-based if x not given).
s	ndarray	Estimated slope of y at the zero crossings.

Source code in pyvoicebox/v_zerocros.py

def v_zerocros(y, m='b', x=None) -> tuple[np.ndarray, np.ndarray]:
    """Find zero crossings in a signal.

    Parameters
    ----------
    y : array_like
        Input waveform.
    m : str, optional
        Mode string:
          'p' - positive crossings only
          'n' - negative crossings only
          'b' - both (default)
          'r' - round to sample values
    x : array_like, optional
        X-axis values for y. Default: 1-based indices (MATLAB convention).

    Returns
    -------
    t : ndarray
        X-axis positions of zero crossings (1-based if x not given).
    s : ndarray
        Estimated slope of y at the zero crossings.
    """
    y = np.asarray(y, dtype=float).ravel()
    s_sign = (y >= 0).astype(int)
    k = s_sign[1:] - s_sign[:-1]

    if 'p' in m:
        f = np.where(k > 0)[0]
    elif 'n' in m:
        f = np.where(k < 0)[0]
    else:
        f = np.where(k != 0)[0]

    s = y[f + 1] - y[f]
    # t uses 1-based indexing by default (MATLAB convention)
    t = (f + 1) - y[f] / s  # f is 0-based, add 1 for MATLAB convention

    if 'r' in m:
        t = np.round(t)

    if x is not None:
        x = np.asarray(x, dtype=float).ravel()
        tf = t - (f + 1)  # fractional sample
        t = x[f] * (1 - tf) + x[f + 1] * tf
        s = s / (x[f + 1] - x[f])

    return t, s

v_schmitt

V_SCHMITT - Pass input signal through a Schmitt trigger.

v_schmitt

v_schmitt(
    x, thresh=0.5, minwid=0, return_transitions=False
) -> ndarray

Pass input signal through a Schmitt trigger.

Parameters:

Name	Type	Description	Default
x	array_like	Input signal.	required
thresh	float or array_like	If scalar in [0,1]: hysteresis fraction. Thresholds are computed from max/min of x. If 2-element: [low, high] thresholds. Default is 0.5.	0.5
minwid	int	Minimum pulse width in samples. Default is 0.	0
return_transitions	bool	If True, return (y_transitions, t) where y_transitions contains alternating +1/-1 and t contains transition sample indices. If False (default), return full y signal.	False

Returns:

Name	Type	Description
y	ndarray	If return_transitions=False: signal of same length as x with values -1, 0, or +1. If return_transitions=True: alternating +1/-1 transition values.
t	ndarray (only if return_transitions=True)	Sample indices where x crossed the thresholds (1-based).

Source code in pyvoicebox/v_schmitt.py

def v_schmitt(x, thresh=0.5, minwid=0, return_transitions=False) -> np.ndarray:
    """Pass input signal through a Schmitt trigger.

    Parameters
    ----------
    x : array_like
        Input signal.
    thresh : float or array_like, optional
        If scalar in [0,1]: hysteresis fraction. Thresholds are computed from
        max/min of x. If 2-element: [low, high] thresholds. Default is 0.5.
    minwid : int, optional
        Minimum pulse width in samples. Default is 0.
    return_transitions : bool, optional
        If True, return (y_transitions, t) where y_transitions contains
        alternating +1/-1 and t contains transition sample indices.
        If False (default), return full y signal.

    Returns
    -------
    y : ndarray
        If return_transitions=False: signal of same length as x with
        values -1, 0, or +1.
        If return_transitions=True: alternating +1/-1 transition values.
    t : ndarray (only if return_transitions=True)
        Sample indices where x crossed the thresholds (1-based).
    """
    x = np.asarray(x, dtype=float).ravel()
    thresh = np.asarray(thresh, dtype=float).ravel()

    if len(thresh) < 2:
        xmax = np.max(x)
        xmin = np.min(x)
        low_offset = (xmax - xmin) * (1 - thresh[0]) / 2.0
        high = xmax - low_offset
        low = xmin + low_offset
    else:
        low = thresh[0]
        high = thresh[1]

    c = (x > high).astype(int) - (x < low).astype(int)
    # Zero out entries where the value equals the previous
    c[1:] = c[1:] * (c[1:] != c[:-1]).astype(int)

    t = np.where(c != 0)[0]
    if len(t) < 2:
        if return_transitions:
            if len(t) == 0:
                return np.array([], dtype=int), np.array([], dtype=int)
            return c[t], t + 1  # 1-based
        y = np.zeros(len(x), dtype=int)
        if len(t) > 0:
            y[t[0]:] = c[t[0]]
        return y

    # Remove duplicates (consecutive same-sign transitions)
    to_remove = []
    i = 0
    while i < len(t) - 1:
        if c[t[i + 1]] == c[t[i]]:
            to_remove.append(i + 1)
        i += 1
    t = np.delete(t, to_remove)

    # Apply minimum width constraint
    if minwid >= 1:
        # Remove transitions that are too close together
        to_remove = np.where(np.diff(t) < minwid)[0]
        t = np.delete(t, to_remove)
        # Remove duplicates again
        to_remove2 = []
        for i in range(len(t) - 1):
            if c[t[i + 1]] == c[t[i]]:
                to_remove2.append(i + 1)
        if to_remove2:
            t = np.delete(t, to_remove2)

    if return_transitions:
        return c[t], t + 1  # 1-based

    y = np.zeros(len(x), dtype=int)
    if len(t) > 0:
        y_cumul = np.zeros(len(x), dtype=int)
        y_cumul[t] = 2 * c[t]
        y_cumul[t[0]] = c[t[0]]
        y = np.cumsum(y_cumul)

    return y

v_teager

V_TEAGER - Calculate Teager energy waveform.

v_teager

v_teager(x, d=None, m='') -> ndarray

Calculate Teager energy waveform.

y(n) = abs(x(n))^2 - x(n+1)*conj(x(n-1))

Parameters:

Name	Type	Description	Default
x	array_like	Input signal.	required
d	int or None	Dimension to apply filter along (0-based). Default is first non-singleton dimension.	None
m	str	'x' to suppress extrapolation of end points. Output will be two samples shorter than input.	''

Returns:

Name	Type	Description
y	ndarray	Teager energy waveform.

Source code in pyvoicebox/v_teager.py

def v_teager(x, d=None, m='') -> np.ndarray:
    """Calculate Teager energy waveform.

    y(n) = abs(x(n))^2 - x(n+1)*conj(x(n-1))

    Parameters
    ----------
    x : array_like
        Input signal.
    d : int or None, optional
        Dimension to apply filter along (0-based).
        Default is first non-singleton dimension.
    m : str, optional
        'x' to suppress extrapolation of end points. Output will be
        two samples shorter than input.

    Returns
    -------
    y : ndarray
        Teager energy waveform.
    """
    x = np.asarray(x, dtype=float)
    e = x.shape
    p = x.size

    if d is None:
        d_found = None
        for i, si in enumerate(e):
            if si > 1:
                d_found = i
                break
        d = d_found if d_found is not None else 0

    k = e[d]  # size of active dimension
    q = p // k  # size of remainder

    # Move working dimension to front and reshape to 2D
    z = np.moveaxis(x, d, 0)
    r = z.shape
    z = z.reshape(k, q)

    if 'x' in m:
        y = z[1:k - 1, :] * np.conj(z[1:k - 1, :]) - z[2:k, :] * np.conj(z[0:k - 2, :])
        k_out = k - 2
    elif k >= 4:
        y = np.zeros((k, q))
        y[1:k - 1, :] = z[1:k - 1, :] * np.conj(z[1:k - 1, :]) - z[2:k, :] * np.conj(z[0:k - 2, :])
        y[0, :] = 2 * y[1, :] - y[2, :]
        y[k - 1, :] = 2 * y[k - 2, :] - y[k - 3, :]
        k_out = k
    elif k == 3:
        val = z[1, :] * np.conj(z[1, :]) - z[2, :] * np.conj(z[0, :])
        y = np.tile(val, (3, 1))
        k_out = k
    else:
        y = np.zeros((k, q))
        k_out = k

    # Reshape back
    r_out = list(r)
    r_out[0] = k_out
    y = y.reshape(r_out)
    y = np.moveaxis(y, 0, d)

    return np.real(y)

v_ditherq

V_DITHERQ - Add dither and quantize.

v_ditherq

v_ditherq(
    x, m="w", zi=None, rng=None
) -> tuple[ndarray, float]

Add dither and quantize.

Parameters:

Name	Type	Description	Default
x	array_like	Input signal.	required
m	str	Mode: 'w' white dither (default), 'h' high-pass dither, 'l' low-pass dither, 'n' no dither.	'w'
zi	float	Initial state for filtering. Default is random.	None
rng	Generator	Random number generator for reproducibility. Default uses np.random.	None

Returns:

Name	Type	Description
y	ndarray	Dithered and quantized signal.
zf	float	Output state.

Source code in pyvoicebox/v_ditherq.py

def v_ditherq(x, m='w', zi=None, rng=None) -> tuple[np.ndarray, float]:
    """Add dither and quantize.

    Parameters
    ----------
    x : array_like
        Input signal.
    m : str, optional
        Mode: 'w' white dither (default), 'h' high-pass dither,
        'l' low-pass dither, 'n' no dither.
    zi : float, optional
        Initial state for filtering. Default is random.
    rng : numpy.random.Generator, optional
        Random number generator for reproducibility. Default uses np.random.

    Returns
    -------
    y : ndarray
        Dithered and quantized signal.
    zf : float
        Output state.
    """
    x = np.asarray(x, dtype=float)
    was_row = (x.ndim == 1 or (x.ndim == 2 and x.shape[0] == 1))
    xflat = x.ravel()
    n = len(xflat)

    if rng is None:
        rng = np.random.default_rng()

    if zi is None:
        zi = rng.random()

    if 'n' in m:
        y = _matlab_round(xflat)
        zf = zi
    elif 'h' in m or 'l' in m:
        v = rng.random(n + 1)
        v[0] = zi
        zf = v[-1]
        if 'h' in m:
            y = _matlab_round(xflat + v[1:] - v[:-1])
        else:
            y = _matlab_round(xflat + v[1:] + v[:-1] - 1)
    else:
        # White dither (default)
        r = rng.random((n, 2))
        y = _matlab_round(xflat + r[:, 0] - r[:, 1])
        zf = rng.random()

    if was_row and x.ndim >= 2 and x.shape[0] == 1:
        y = y.reshape(1, -1)
    elif x.ndim == 0:
        y = y.reshape(())
    else:
        y = y.reshape(x.shape)

    return y, zf

Array helpers

v_nearnonz

V_NEARNONZ - Replace each zero element with nearest non-zero element.

v_nearnonz

v_nearnonz(x, d=None) -> tuple[ndarray, ndarray, ndarray]

Replace each zero element with the nearest non-zero element.

Parameters:

Name	Type	Description	Default
x	array_like	Input vector, matrix or larger array.	required
d	int or None	Dimension to apply filter along (0-based). Default is first non-singleton dimension.	None

Returns:

Name	Type	Description
v	ndarray	Same shape as x, with zeros replaced by nearest non-zero value.
y	ndarray	Same shape as x, index along dimension d from which v was taken (1-based).
w	ndarray	Same shape as x, distance to the nearest non-zero entry.

Source code in pyvoicebox/v_nearnonz.py

def v_nearnonz(x, d=None) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Replace each zero element with the nearest non-zero element.

    Parameters
    ----------
    x : array_like
        Input vector, matrix or larger array.
    d : int or None, optional
        Dimension to apply filter along (0-based).
        Default is first non-singleton dimension.

    Returns
    -------
    v : ndarray
        Same shape as x, with zeros replaced by nearest non-zero value.
    y : ndarray
        Same shape as x, index along dimension d from which v was taken (1-based).
    w : ndarray
        Same shape as x, distance to the nearest non-zero entry.
    """
    x = np.asarray(x, dtype=float)
    e = x.shape
    p = x.size

    if d is None:
        d_found = None
        for i, si in enumerate(e):
            if si > 1:
                d_found = i
                break
        d = d_found if d_found is not None else 0

    k = e[d]  # size of active dimension
    q = p // k  # size of remainder

    # Move working dimension to front and reshape to 2D
    z = np.moveaxis(x, d, 0)
    r = z.shape
    z = z.reshape(k, q)

    xx = (z != 0)  # boolean array of non-zero positions
    cx = np.cumsum(xx, axis=0)  # cumulative count of non-zeros

    v_out = np.zeros((k, q))
    y_out = np.zeros((k, q), dtype=int)
    w_out = np.zeros((k, q), dtype=int)

    for col in range(q):
        nz_indices = np.where(xx[:, col])[0]
        if len(nz_indices) == 0:
            # No non-zero entries: y=0, v=0, w=0
            continue

        # Build the position array: for each element find the nearest non-zero
        for row in range(k):
            # Find nearest non-zero
            dists = np.abs(nz_indices - row)
            # Tie-break: higher index wins (>=)
            min_dist = np.min(dists)
            candidates = nz_indices[dists == min_dist]
            nearest = candidates[-1]  # take higher index for tie-break
            y_out[row, col] = nearest + 1  # 1-based
            v_out[row, col] = z[nearest, col]
            w_out[row, col] = nearest - row

    # Reshape back
    v_out = v_out.reshape(r)
    y_out = y_out.reshape(r)
    w_out = w_out.reshape(r)

    v_out = np.moveaxis(v_out, 0, d)
    y_out = np.moveaxis(y_out, 0, d)
    w_out = np.moveaxis(w_out, 0, d)

    return v_out, y_out, w_out

v_rangelim

V_RANGELIM - Limit the range of matrix elements.

v_rangelim

v_rangelim(x, r, m='lp') -> ndarray

Limit the range of matrix elements.

Parameters:

Name	Type	Description	Default
x	array_like	Input data.	required
r	float or array_like	If 2-element: explicit [min, max] limits. If scalar: range specification depending on mode.	required
m	str	Mode string: 'd' - range r in dB: 20*log10(max/min) 'r' - range r is max/min ratio 'l' - range r is max-min difference (default) 'p' - max(x) is top of range (default) 't' - min(x) is bottom of range 'g' - geometric mean is centre of range 'u' - mean is centre of range 'm' - median is centre of range 'c' - clip out-of-range values (default) 'n' - set out-of-range values to NaN	'lp'

Returns:

Name	Type	Description
y	ndarray	Output data, same shape as x.

Source code in pyvoicebox/v_rangelim.py

def v_rangelim(x, r, m='lp') -> np.ndarray:
    """Limit the range of matrix elements.

    Parameters
    ----------
    x : array_like
        Input data.
    r : float or array_like
        If 2-element: explicit [min, max] limits.
        If scalar: range specification depending on mode.
    m : str, optional
        Mode string:
          'd' - range r in dB: 20*log10(max/min)
          'r' - range r is max/min ratio
          'l' - range r is max-min difference (default)
          'p' - max(x) is top of range (default)
          't' - min(x) is bottom of range
          'g' - geometric mean is centre of range
          'u' - mean is centre of range
          'm' - median is centre of range
          'c' - clip out-of-range values (default)
          'n' - set out-of-range values to NaN

    Returns
    -------
    y : ndarray
        Output data, same shape as x.
    """
    x = np.asarray(x, dtype=float)
    r = np.atleast_1d(np.asarray(r, dtype=float))

    if len(r) > 1:
        p = r[0]
        q = r[1]
    else:
        r_val = r[0]

        # Determine reference type: g=1, u=2, m=3, t=4, p=5
        ref_map = {'g': 1, 'u': 2, 'm': 3, 't': 4}
        ref_id = 5  # default 'p'
        for char, val in ref_map.items():
            if char in m:
                ref_id = val
                break

        # Determine range type: d=1, r=2, l=3
        if 'd' in m:
            ir = 1
        elif 'r' in m:
            ir = 2
        else:
            ir = 3

        if ir == 1:  # dB
            r_val = 10 ** (0.05 * r_val)

        xflat = x.ravel()

        if ir == 3:  # linear range
            if ref_id == 5:  # 'p': peak
                q = np.max(xflat)
                p = q - r_val
            elif ref_id == 4:  # 't': trough
                p = np.min(xflat)
                q = p + r_val
            elif ref_id == 3:  # 'm': median
                p = np.median(xflat) - 0.5 * r_val
                q = p + r_val
            elif ref_id == 2:  # 'u': mean
                p = np.mean(xflat) - 0.5 * r_val
                q = p + r_val
            elif ref_id == 1:  # 'g': geometric mean
                p = np.exp(np.mean(np.log(xflat))) - 0.5 * r_val
                q = p + r_val
        else:  # ratio range ('r' or 'd')
            if ref_id == 5:  # 'p': peak
                q = np.max(xflat)
                p = q / r_val
            elif ref_id == 4:  # 't': trough
                p = np.min(xflat)
                q = p * r_val
            elif ref_id == 3:  # 'm': median
                p = np.median(xflat) / np.sqrt(r_val)
                q = p * r_val
            elif ref_id == 2:  # 'u': mean
                p = np.mean(xflat) / np.sqrt(r_val)
                q = p * r_val
            elif ref_id == 1:  # 'g': geometric mean
                p = np.exp(np.mean(np.log(xflat))) / np.sqrt(r_val)
                q = p * r_val

    y = x.copy()
    if 'n' in m:
        y[(x < p) | (x > q)] = np.nan
    else:
        y[x < p] = p
        y[x > q] = q

    return y

v_horizdiff

V_HORIZDIFF - Estimate horizontal difference between two functions.

v_horizdiff

v_horizdiff(y, v, x=None, u=None, q='') -> ndarray

Estimate horizontal difference between two functions of x.

Approximately: y(x) = v(x+z).

Parameters:

Name	Type	Description	Default
y	(array_like, shape(n, m))	Each column is a function of x.	required
v	(array_like, shape(k))	Reference function.	required
x	(array_like, shape(n))	x values for y. Default (1:n).	None
u	(array_like, shape(k))	x values for v. Default same as x.	None
q	str	Interpolation mode: 'l' linear, 'p' pchip, 's' spline.	''

Returns:

Name	Type	Description
z	(ndarray, shape(n, m))	Horizontal difference.
zm	(ndarray, shape(m))	MMSE horizontal difference.

Source code in pyvoicebox/v_horizdiff.py

def v_horizdiff(y, v, x=None, u=None, q='') -> np.ndarray:
    """Estimate horizontal difference between two functions of x.

    Approximately: y(x) = v(x+z).

    Parameters
    ----------
    y : array_like, shape (n, m)
        Each column is a function of x.
    v : array_like, shape (k,)
        Reference function.
    x : array_like, shape (n,), optional
        x values for y. Default (1:n).
    u : array_like, shape (k,), optional
        x values for v. Default same as x.
    q : str, optional
        Interpolation mode: 'l' linear, 'p' pchip, 's' spline.

    Returns
    -------
    z : ndarray, shape (n, m)
        Horizontal difference.
    zm : ndarray, shape (m,)
        MMSE horizontal difference.
    """
    y = np.atleast_2d(np.asarray(y, dtype=float))
    if y.shape[0] == 1 and y.shape[1] > 1:
        y = y.T
    n, m_cols = y.shape
    v = np.asarray(v, dtype=float).ravel()

    if x is None:
        x = np.arange(1, n + 1, dtype=float)
    else:
        x = np.asarray(x, dtype=float).ravel()

    if u is None:
        u = x.copy()
    else:
        u = np.asarray(u, dtype=float).ravel()

    # Choose interpolation method
    if n >= 4 and 's' in q:
        kind = 'cubic'
    elif n >= 2 and ('s' in q or 'p' in q):
        kind = 'cubic'
    else:
        kind = 'linear'

    # Interpolate inverse function u(v)
    f_inv = interp1d(v, u, kind=kind, fill_value='extrapolate')
    z = f_inv(y) - x[:, np.newaxis]

    # MMSE horizontal difference
    zm = np.zeros(m_cols)
    f_fwd = interp1d(u, v, kind=kind, fill_value='extrapolate')
    for i in range(m_cols):
        def objective(zm_val):
            return np.sum((y[:, i] - f_fwd(x + zm_val)) ** 2)
        result = minimize_scalar(objective, bounds=(u[0] - x[-1], u[-1] - x[0]), method='bounded')
        zm[i] = result.x

    return z, zm

v_interval

V_INTERVAL - Classify X values into contiguous intervals.

v_interval

v_interval(x, y, m='') -> tuple[ndarray, ndarray]

Classify X values into contiguous intervals with boundaries from Y.

Parameters:

Name	Type	Description	Default
x	array_like	Vector of test values.	required
y	array_like	Vector of monotonically increasing interval boundaries. Interval i is [y[i], y[i+1]).	required
m	str	Mode options: For x < y[0]: 'e' - extrapolate: i=1, f<0 (default) 'c' - clip: i=1, f=0 'n' - NaN: i=NaN, f=NaN 'z' - zero: i=0, f<0 For x >= y[-1]: 'E' - extrapolate: i=ny-1, f>1 (default) 'C' - clip: i=ny-1, f=1 'N' - NaN: i=NaN, f=NaN 'Z' - zero: i=ny, f>1	''

Returns:

Name	Type	Description
i	ndarray	Interval indices (1-based, matching MATLAB convention).
f	ndarray	Fractional position within the interval.

Source code in pyvoicebox/v_interval.py

def v_interval(x, y, m='') -> tuple[np.ndarray, np.ndarray]:
    """Classify X values into contiguous intervals with boundaries from Y.

    Parameters
    ----------
    x : array_like
        Vector of test values.
    y : array_like
        Vector of monotonically increasing interval boundaries.
        Interval i is [y[i], y[i+1]).
    m : str, optional
        Mode options:
          For x < y[0]:
            'e' - extrapolate: i=1, f<0 (default)
            'c' - clip: i=1, f=0
            'n' - NaN: i=NaN, f=NaN
            'z' - zero: i=0, f<0
          For x >= y[-1]:
            'E' - extrapolate: i=ny-1, f>1 (default)
            'C' - clip: i=ny-1, f=1
            'N' - NaN: i=NaN, f=NaN
            'Z' - zero: i=ny, f>1

    Returns
    -------
    i : ndarray
        Interval indices (1-based, matching MATLAB convention).
    f : ndarray
        Fractional position within the interval.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float).ravel()
    x_shape = x.shape
    x_flat = x.ravel()
    ny = len(y)

    # Use searchsorted to find where each x would be inserted in y
    # This gives the index of the first element in y that is > x
    # k = number of y values <= x, equivalent to MATLAB's logic
    k = np.searchsorted(y, x_flat, side='right')  # k in range [0, ny]

    # Force i to lie in range [1, ny-1] (1-based)
    i = np.clip(k, 1, ny - 1)

    # Fractional position: f = (x - y[i-1]) / (y[i] - y[i-1])
    # i is 1-based, so y[i-1] and y[i] in 0-based indexing
    f = (x_flat - y[i - 1]) / (y[i] - y[i - 1])

    # Handle values below range
    klo = k < 1
    if np.any(klo):
        if 'c' in m:
            f[klo] = 0
        elif 'n' in m:
            i = i.astype(float)
            i[klo] = np.nan
            f[klo] = np.nan
        elif 'z' in m:
            i = i.astype(float)
            i[klo] = 0

    # Handle values above range
    khi = k >= ny
    if np.any(khi):
        if 'C' in m:
            f[khi] = 1
        elif 'N' in m:
            if not np.issubdtype(i.dtype, np.floating):
                i = i.astype(float)
            i[khi] = np.nan
            f[khi] = np.nan
        elif 'Z' in m:
            if not np.issubdtype(i.dtype, np.floating):
                i = i.astype(float)
            i[khi] = ny

    i = np.reshape(i, x_shape)
    f = np.reshape(f, x_shape)

    return i, f

v_modsym

V_MODSYM - Symmetric modulus function.

v_modsym

v_modsym(x, y=1, r=None) -> tuple[ndarray, ndarray]

Symmetric modulus function.

Adds an integer multiple of y onto x so that it lies in the range [r-y/2, r+y/2) if y is positive or (r-y/2, r+y/2] if y is negative.

Parameters:

Name	Type	Description	Default
x	array_like	Input data.	required
y	float or array_like	Modulus. Default is 1.	1
r	float or array_like	Reference data. Default is 0.	None

Returns:

Name	Type	Description
z	ndarray	Output data.
k	ndarray	Integer multiple of y that was added.

Source code in pyvoicebox/v_modsym.py

def v_modsym(x, y=1, r=None) -> tuple[np.ndarray, np.ndarray]:
    """Symmetric modulus function.

    Adds an integer multiple of y onto x so that it lies in the range
    [r-y/2, r+y/2) if y is positive or (r-y/2, r+y/2] if y is negative.

    Parameters
    ----------
    x : array_like
        Input data.
    y : float or array_like, optional
        Modulus. Default is 1.
    r : float or array_like, optional
        Reference data. Default is 0.

    Returns
    -------
    z : ndarray
        Output data.
    k : ndarray
        Integer multiple of y that was added.
    """
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)
    if r is None:
        v = 0.5 * y
    else:
        r = np.asarray(r, dtype=float)
        v = 0.5 * y - r
    z = np.mod(x + v, y) - v
    k = np.round((z - x) / y).astype(int)
    return z, k

v_zerotrim

V_ZEROTRIM - Remove trailing zero rows and columns.

v_zerotrim

v_zerotrim(x) -> ndarray

Remove trailing zero rows and columns from a matrix.

Parameters:

Name	Type	Description	Default
x	array_like	Input matrix.	required

Returns:

Name	Type	Description
z	ndarray	Trimmed matrix, or empty array if all zeros.

Source code in pyvoicebox/v_zerotrim.py

def v_zerotrim(x) -> np.ndarray:
    """Remove trailing zero rows and columns from a matrix.

    Parameters
    ----------
    x : array_like
        Input matrix.

    Returns
    -------
    z : ndarray
        Trimmed matrix, or empty array if all zeros.
    """
    x = np.asarray(x)
    if x.ndim == 1:
        nz = np.nonzero(x)[0]
        if len(nz) == 0:
            return np.array([])
        return x[:nz[-1] + 1]

    # Find last nonzero column
    col_any = np.any(x, axis=0)
    col_nz = np.nonzero(col_any)[0]
    if len(col_nz) == 0:
        return np.array([]).reshape(0, 0)
    c = col_nz[-1] + 1

    # Find last nonzero row
    row_any = np.any(x, axis=1)
    row_nz = np.nonzero(row_any)[0]
    r = row_nz[-1] + 1

    return x[:r, :c]