tessl/pypi-spectrum

Spectrum Analysis Tools - contains tools to estimate Power Spectral Densities using methods based on Fourier transform, Parametric methods or eigenvalues analysis

—

Pending

Overview

Eval results

Files

Linear Prediction

Name: tessl/pypi-spectrum
Author: tessl

Linear prediction analysis and parameter conversion utilities for autoregressive modeling and speech processing applications. These functions provide comprehensive conversion between different parameter representations used in digital signal processing and linear predictive coding.

Capabilities

Levinson-Durbin Algorithm

Core algorithm for solving the Yule-Walker equations to obtain AR parameters from autocorrelation sequences.

def LEVINSON(r, order=None, allow_singularity=False):
    """
    Levinson-Durbin recursion for AR parameter estimation.
    
    Parameters:
    - r: array-like, autocorrelation sequence (r[0] is zero-lag)
    - order: int, AR model order (None uses len(r)-1)
    - allow_singularity: bool, allow singular autocorrelation matrix
    
    Returns:
    tuple: (AR coefficients array, prediction error float, reflection coefficients array)
    """

def rlevinson(a, efinal):
    """
    Reverse Levinson-Durbin (AR coefficients to autocorrelation).
    
    Parameters:
    - a: array-like, AR coefficients including a0=1
    - efinal: float, final prediction error
    
    Returns:
    array: Autocorrelation sequence
    """

def levdown(anxt, enxt=None):
    """
    Levinson downward recursion.
    
    Parameters:
    - anxt: array-like, AR coefficients at order n+1
    - enxt: float, prediction error at order n+1
    
    Returns:
    tuple: (AR coefficients at order n, prediction error at order n)
    """

def levup(acur, knxt, ecur=None):
    """
    Levinson upward recursion.
    
    Parameters:
    - acur: array-like, AR coefficients at current order
    - knxt: float, next reflection coefficient
    - ecur: float, current prediction error
    
    Returns:
    tuple: (AR coefficients at next order, prediction error at next order)
    """

Linear Predictive Coding

High-level LPC analysis function for speech and audio processing applications.

def lpc(x, N=None):
    """
    Linear Predictive Coding analysis.
    
    Parameters:
    - x: array-like, input time series data
    - N: int, LPC order (None for automatic selection)
    
    Returns:
    tuple: (LPC coefficients array, prediction error float)
    """

Autocorrelation Conversions

Functions for converting between autocorrelation sequences and polynomial coefficients.

def ac2poly(data):
    """
    Autocorrelation to polynomial coefficient conversion.
    
    Parameters:
    - data: array-like, autocorrelation sequence
    
    Returns:
    array: Polynomial coefficients
    """

def poly2ac(poly, efinal):
    """
    Polynomial coefficients to autocorrelation conversion.
    
    Parameters:
    - poly: array-like, polynomial coefficients
    - efinal: float, final prediction error
    
    Returns:
    array: Autocorrelation sequence
    """

def ac2rc(data):
    """
    Autocorrelation to reflection coefficients conversion.
    
    Parameters:
    - data: array-like, autocorrelation sequence
    
    Returns:
    array: Reflection coefficients
    """

def rc2ac(k, R0):
    """
    Reflection coefficients to autocorrelation conversion.
    
    Parameters:
    - k: array-like, reflection coefficients
    - R0: float, zero-lag autocorrelation
    
    Returns:
    array: Autocorrelation sequence
    """

Reflection Coefficient Conversions

Functions for converting between AR parameters and reflection coefficients (PARCOR coefficients).

def ar2rc(ar):
    """
    AR parameters to reflection coefficients conversion.
    
    Parameters:
    - ar: array-like, AR coefficients (without leading 1)
    
    Returns:
    array: Reflection coefficients
    """

def rc2poly(kr, r0=None):
    """
    Reflection coefficients to polynomial coefficients conversion.
    
    Parameters:
    - kr: array-like, reflection coefficients
    - r0: float, zero-lag autocorrelation (optional)
    
    Returns:
    array: Polynomial coefficients
    """

def poly2rc(a, efinal):
    """
    Polynomial coefficients to reflection coefficients conversion.
    
    Parameters:
    - a: array-like, polynomial coefficients
    - efinal: float, final prediction error
    
    Returns:
    array: Reflection coefficients
    """

Line Spectral Frequencies

Conversion between polynomial coefficients and line spectral frequencies (LSF), useful for quantization and stability analysis.

def poly2lsf(a):
    """
    Polynomial coefficients to line spectral frequencies conversion.
    
    Parameters:
    - a: array-like, polynomial coefficients
    
    Returns:
    array: Line spectral frequencies (normalized to [0, π])
    """

def lsf2poly(lsf):
    """
    Line spectral frequencies to polynomial coefficients conversion.
    
    Parameters:
    - lsf: array-like, line spectral frequencies
    
    Returns:
    array: Polynomial coefficients
    """

Alternative Parameterizations

Conversions involving inverse sine parameters and log area ratios for robust parameter representation.

def is2rc(inv_sin):
    """
    Inverse sine parameters to reflection coefficients conversion.
    
    Parameters:
    - inv_sin: array-like, inverse sine parameters
    
    Returns:
    array: Reflection coefficients
    """

def rc2is(k):
    """
    Reflection coefficients to inverse sine parameters conversion.
    
    Parameters:
    - k: array-like, reflection coefficients
    
    Returns:
    array: Inverse sine parameters
    """

def rc2lar(k):
    """
    Reflection coefficients to log area ratios conversion.
    
    Parameters:
    - k: array-like, reflection coefficients
    
    Returns:
    array: Log area ratios
    """

def lar2rc(g):
    """
    Log area ratios to reflection coefficients conversion.
    
    Parameters:
    - g: array-like, log area ratios
    
    Returns:
    array: Reflection coefficients
    """

Usage Examples

Basic Linear Prediction Analysis

import spectrum
import numpy as np
import matplotlib.pyplot as plt

# Generate AR signal for LPC analysis
N = 256
# Create AR(3) signal with known coefficients
ar_true = np.array([1.0, -2.2137, 2.9403, -1.7720])  # Includes leading 1
noise = np.random.randn(N)

# Generate signal using difference equation
signal = np.zeros(N)
signal[:3] = noise[:3]  # Initialize
for n in range(3, N):
    signal[n] = -sum(ar_true[1:] * signal[n-3:n][::-1]) + noise[n]

# LPC analysis
lpc_order = 3
lpc_coeffs, pred_error = spectrum.lpc(signal, N=lpc_order)

print(f"True AR coefficients: {-ar_true[1:]}")
print(f"LPC coefficients: {lpc_coeffs[1:]}")
print(f"Prediction error: {pred_error}")

# Compare with autocorrelation method
autocorr = spectrum.CORRELATION(signal, maxlags=lpc_order, norm='biased')
ar_levinson, err_levinson, rc_levinson = spectrum.LEVINSON(autocorr, order=lpc_order)

print(f"Levinson AR coeffs: {ar_levinson}")
print(f"Levinson error: {err_levinson}")
print(f"Reflection coeffs: {rc_levinson}")

Parameter Conversion Chain

import spectrum
import numpy as np

# Start with AR coefficients
ar_coeffs = np.array([0.5, -0.3, 0.1])  # AR(3) coefficients
print(f"Original AR coefficients: {ar_coeffs}")

# Convert AR to reflection coefficients
rc = spectrum.ar2rc(ar_coeffs)
print(f"Reflection coefficients: {rc}")

# Convert reflection coefficients back to AR
ar_reconstructed = spectrum.rc2poly(rc)[1:]  # Remove leading 1
print(f"Reconstructed AR: {ar_reconstructed}")

# Convert to line spectral frequencies
poly_with_one = np.concatenate([[1], ar_coeffs])
lsf = spectrum.poly2lsf(poly_with_one)
print(f"Line spectral frequencies: {lsf}")

# Convert LSF back to polynomial
poly_reconstructed = spectrum.lsf2poly(lsf)
ar_from_lsf = poly_reconstructed[1:]
print(f"AR from LSF: {ar_from_lsf}")

# Convert to log area ratios
lar = spectrum.rc2lar(rc)
print(f"Log area ratios: {lar}")

# Convert LAR back to reflection coefficients
rc_from_lar = spectrum.lar2rc(lar)
print(f"RC from LAR: {rc_from_lar}")

# Verify round-trip conversion accuracy
print(f"\nConversion errors:")
print(f"AR round-trip error: {np.max(np.abs(ar_coeffs - ar_reconstructed))}")
print(f"AR from LSF error: {np.max(np.abs(ar_coeffs - ar_from_lsf))}")
print(f"RC from LAR error: {np.max(np.abs(rc - rc_from_lar))}")

Autocorrelation-Based Analysis

import spectrum
import numpy as np
import matplotlib.pyplot as plt

# Generate test signal
fs = 1000
N = 500
t = np.arange(N) / fs

# Multi-component signal
signal = (np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*120*t) + 
         0.2*np.random.randn(N))

# Compute autocorrelation with different methods
max_lag = 20
autocorr_biased = spectrum.CORRELATION(signal, maxlags=max_lag, norm='biased')
autocorr_unbiased = spectrum.CORRELATION(signal, maxlags=max_lag, norm='unbiased')

# LPC analysis using autocorrelation
orders = [5, 10, 15, 20]
plt.figure(figsize=(15, 10))

# Plot signal and autocorrelation
plt.subplot(2, 3, 1)
plt.plot(t[:200], signal[:200])
plt.title('Input Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

plt.subplot(2, 3, 2)
lags = np.arange(-max_lag, max_lag+1)
plt.plot(lags, autocorr_biased, 'b-', label='Biased')
plt.plot(lags, autocorr_unbiased, 'r--', label='Unbiased')
plt.title('Autocorrelation Function')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.legend()
plt.grid(True)

# LPC analysis for different orders
for i, order in enumerate(orders):
    plt.subplot(2, 3, 3+i)
    
    # Use positive lags for Levinson-Durbin
    autocorr_pos = autocorr_biased[max_lag:]  # Take positive lags including zero
    ar_coeffs, pred_error, rc = spectrum.LEVINSON(autocorr_pos, order=order)
    
    # Compute and plot frequency response
    w = np.linspace(0, np.pi, 512)
    # H(w) = 1/A(w) where A(w) is AR polynomial
    A_w = np.polyval(np.concatenate([[1], ar_coeffs]), np.exp(1j*w))
    H_w = 1 / A_w
    psd_ar = np.abs(H_w)**2 * pred_error
    
    freqs_hz = w * fs / (2*np.pi)
    plt.semilogy(freqs_hz, psd_ar, 'b-', linewidth=2)
    plt.title(f'LPC Spectrum (Order {order})')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('PSD')
    plt.xlim(0, fs/2)
    plt.grid(True)
    
    # Print some statistics
    print(f"Order {order}: Pred error = {pred_error:.6f}, RC = {rc[:3]}")

plt.tight_layout()
plt.show()

Stability Analysis Using Different Parameterizations

import spectrum
import numpy as np
import matplotlib.pyplot as plt

def check_stability(ar_coeffs):
    """Check if AR filter is stable (all poles inside unit circle)."""
    # Poles are roots of AR polynomial
    poly = np.concatenate([[1], ar_coeffs])
    poles = np.roots(poly)
    return np.all(np.abs(poles) < 1), poles

# Generate series of AR coefficients with varying stability
test_cases = [
    np.array([0.5, -0.2]),           # Stable
    np.array([0.8, -0.6]),           # Stable  
    np.array([1.2, -0.5]),           # Unstable
    np.array([-1.8, 0.9]),           # Marginally stable
    np.array([0.9, 0.9, -0.5])      # Higher order
]

plt.figure(figsize=(15, 12))

for i, ar_coeffs in enumerate(test_cases):
    # Check stability
    is_stable, poles = check_stability(ar_coeffs)
    
    # Convert to different representations
    rc = spectrum.ar2rc(ar_coeffs)
    poly_with_one = np.concatenate([[1], ar_coeffs])
    
    try:
        lsf = spectrum.poly2lsf(poly_with_one)
        lar = spectrum.rc2lar(rc)
        
        # Plot pole-zero diagram
        plt.subplot(len(test_cases), 4, 4*i + 1)
        plt.plot(np.real(poles), np.imag(poles), 'ro', markersize=8)
        circle = plt.Circle((0, 0), 1, fill=False, linestyle='--', color='blue')
        plt.gca().add_patch(circle)
        plt.axis('equal')
        plt.xlim(-1.5, 1.5)
        plt.ylim(-1.5, 1.5)
        plt.title(f'Case {i+1}: {"Stable" if is_stable else "Unstable"}')
        plt.xlabel('Real')
        plt.ylabel('Imag')
        plt.grid(True)
        
        # Plot reflection coefficients
        plt.subplot(len(test_cases), 4, 4*i + 2)
        plt.stem(range(len(rc)), rc, basefmt=' ')
        plt.axhline(1, color='r', linestyle='--', alpha=0.7)
        plt.axhline(-1, color='r', linestyle='--', alpha=0.7)
        plt.title('Reflection Coefficients')
        plt.ylabel('RC Value')
        plt.grid(True)
        
        # Plot line spectral frequencies
        plt.subplot(len(test_cases), 4, 4*i + 3)
        plt.stem(range(len(lsf)), lsf/np.pi, basefmt=' ')
        plt.title('Line Spectral Frequencies')
        plt.ylabel('LSF (normalized)')
        plt.ylim(0, 1)
        plt.grid(True)
        
        # Plot log area ratios
        plt.subplot(len(test_cases), 4, 4*i + 4)
        plt.stem(range(len(lar)), lar, basefmt=' ')
        plt.title('Log Area Ratios')
        plt.ylabel('LAR')
        plt.grid(True)
        
        # Print parameter ranges for stability analysis
        print(f"\nCase {i+1} ({'Stable' if is_stable else 'Unstable'}):")
        print(f"  AR coeffs: {ar_coeffs}")
        print(f"  RC range: [{np.min(rc):.3f}, {np.max(rc):.3f}]")
        print(f"  RC stability: {np.all(np.abs(rc) < 1)}")
        print(f"  Poles: {poles}")
        
    except Exception as e:
        print(f"Error processing case {i+1}: {e}")

plt.tight_layout()
plt.show()

# Demonstrate quantization robustness
print("\n" + "="*60)
print("Quantization Robustness Comparison")
print("="*60)

# Original AR coefficients  
ar_orig = np.array([0.7, -0.4, 0.15])
print(f"Original AR: {ar_orig}")

# Add quantization noise
quantization_levels = [8, 6, 4]  # bits
for bits in quantization_levels:
    # Quantize AR coefficients directly
    scale = 2**(bits-1) - 1
    ar_quant = np.round(ar_orig * scale) / scale
    
    # Quantize via reflection coefficients
    rc_orig = spectrum.ar2rc(ar_orig)
    rc_quant = np.round(rc_orig * scale) / scale
    ar_from_rc = spectrum.rc2poly(rc_quant)[1:]
    
    # Quantize via LSF
    poly_orig = np.concatenate([[1], ar_orig])
    lsf_orig = spectrum.poly2lsf(poly_orig)
    lsf_quant = np.round(lsf_orig * scale) / scale
    poly_from_lsf = spectrum.lsf2poly(lsf_quant)
    ar_from_lsf = poly_from_lsf[1:]
    
    # Check stability
    stable_direct = check_stability(ar_quant)[0]
    stable_rc = check_stability(ar_from_rc)[0]
    stable_lsf = check_stability(ar_from_lsf)[0]
    
    print(f"\n{bits}-bit quantization:")
    print(f"  Direct AR:    {ar_quant}, stable: {stable_direct}")
    print(f"  Via RC:       {ar_from_rc}, stable: {stable_rc}")
    print(f"  Via LSF:      {ar_from_lsf}, stable: {stable_lsf}")

Install with Tessl CLI