Mathematical Tools

def CHOLESKY(A, B, method='scipy'):
    """
    Solve linear system AX=B using Cholesky decomposition.
    
    Parameters:
    - A: array-like, Hermitian positive definite coefficient matrix
    - B: array-like, right-hand side matrix or vector
    - method: str, computation method ('scipy', 'numpy', 'numpy_solver')
    
    Returns:
    array: Solution matrix X such that AX = B
    """

def pascal(n):
    """
    Generate Pascal matrix of size n x n.
    
    Parameters:
    - n: int, matrix size
    
    Returns:
    array: Pascal matrix (symmetric positive definite)
    """

def corrmtx(x_input, m, method='autocorrelation'):
    """
    Generate correlation matrix for linear prediction.
    
    Parameters:
    - x_input: array-like, input data vector
    - m: int, matrix order (number of coefficients)
    - method: str, matrix construction method ('autocorrelation', 'covariance', 'prewindowed', 'postwindowed')
    
    Returns:
    array: Correlation matrix for linear prediction analysis
    """

def csvd(A):
    """
    Complex singular value decomposition.
    
    Parameters:
    - A: array-like, input matrix (real or complex)
    
    Returns:
    tuple: (U array, singular values array, V^H array)
    """

def AIC(N, rho, k):
    """
    Akaike Information Criterion for model order selection.
    
    Parameters:
    - N: int, data length
    - rho: float, prediction error (noise variance)
    - k: int, model order (number of parameters)
    
    Returns:
    float: AIC value (lower is better)
    """

def AICc(N, rho, k, norm=True):
    """
    Corrected AIC for small sample sizes.
    
    Parameters:
    - N: int, data length
    - rho: float, prediction error
    - k: int, model order
    - norm: bool, apply normalization factor
    
    Returns:
    float: Corrected AIC value
    """

def MDL(N, rho, k):
    """
    Minimum Description Length criterion.
    
    Parameters:
    - N: int, data length
    - rho: float, prediction error
    - k: int, model order
    
    Returns:
    float: MDL value (lower is better)
    """

def FPE(N, rho, k=None):
    """
    Final Prediction Error criterion.
    
    Parameters:
    - N: int, data length
    - rho: float, prediction error
    - k: int, model order (None for simple FPE)
    
    Returns:
    float: FPE value
    """

def KIC(N, rho, k):
    """
    Kullback Information Criterion.
    
    Parameters:
    - N: int, data length
    - rho: float, prediction error
    - k: int, model order
    
    Returns:
    float: KIC value
    """

def AKICc(N, rho, k):
    """
    Asymmetric KIC corrected for small samples.
    
    Parameters:
    - N: int, data length
    - rho: float, prediction error
    - k: int, model order
    
    Returns:
    float: AKICc value
    """

def CAT(N, rho, k):
    """
    Criterion Autoregressive Transfer function.
    
    Parameters:
    - N: int, data length
    - rho: float, prediction error
    - k: int, model order
    
    Returns:
    float: CAT value
    """

class Criteria(name, N):
    """
    Automatic order selection for parametric methods.
    
    Parameters:
    - name: str, criterion name ('AIC', 'MDL', 'FPE', 'CAT', etc.)
    - N: int, data length
    
    Methods:
    - __call__(rho, k): compute criterion value
    - find_order(rho_values, orders): find optimal order
    """

def tf2zp(b, a):
    """
    Transfer function to zeros and poles conversion.
    
    Parameters:
    - b: array-like, numerator polynomial coefficients
    - a: array-like, denominator polynomial coefficients
    
    Returns:
    tuple: (zeros array, poles array)
    """

def zp2tf(z, p, k):
    """
    Zeros and poles to transfer function conversion.
    
    Parameters:
    - z: array-like, zeros of the transfer function
    - p: array-like, poles of the transfer function  
    - k: float, system gain
    
    Returns:
    tuple: (numerator coefficients array, denominator coefficients array)
    """

def tf2zpk(b, a):
    """
    Transfer function to zeros, poles, and gain conversion.
    
    Parameters:
    - b: array-like, numerator polynomial coefficients
    - a: array-like, denominator polynomial coefficients
    
    Returns:
    tuple: (zeros array, poles array, gain float)
    """

def zpk2tf(z, p, k):
    """
    Zeros, poles, and gain to transfer function conversion.
    
    Parameters:
    - z: array-like, zeros of the transfer function
    - p: array-like, poles of the transfer function
    - k: float, system gain
    
    Returns:
    tuple: (numerator coefficients array, denominator coefficients array)
    """

def eqtflength(b, a):
    """
    Equalize lengths of transfer function numerator and denominator.
    
    Parameters:
    - b: array-like, numerator coefficients
    - a: array-like, denominator coefficients
    
    Returns:
    tuple: (equalized numerator array, equalized denominator array)
    """

def zpk2ss(z, p, k):
    """
    Convert zeros, poles, gain to state-space representation.
    
    Parameters:
    - z: array-like, zeros of the transfer function
    - p: array-like, poles of the transfer function
    - k: float, system gain
    
    Returns:
    tuple: (A matrix, B vector, C vector, D scalar) state-space matrices
    """

def ss2zpk(a, b, c, d, input=0):
    """
    Convert state-space to zeros, poles, gain representation.
    
    Parameters:
    - a: array-like, state matrix
    - b: array-like, input matrix
    - c: array-like, output matrix
    - d: array-like, feedthrough matrix
    - input: int, input channel for MIMO systems
    
    Returns:
    tuple: (zeros array, poles array, gain float)
    """

def tf2latc(num=[1.], den=[1.]):
    """
    Convert transfer function to lattice filter coefficients.
    
    Parameters:
    - num: array-like, numerator coefficients (default [1.])
    - den: array-like, denominator coefficients (default [1.])
    
    Returns:
    array: Lattice reflection coefficients
    """

def allpole2latc(num, den):
    """
    Convert all-pole transfer function to lattice form.
    
    Parameters:
    - num: array-like, numerator coefficients
    - den: array-like, denominator coefficients
    
    Returns:
    array: Lattice coefficients for all-pole system
    """

def HERMTOEP(T0, T, Z):
    """
    Solve Hermitian Toeplitz system efficiently.
    
    Parameters:
    - T0: complex, Toeplitz matrix diagonal element
    - T: array-like, first row/column of Toeplitz matrix (excluding T0)
    - Z: array-like, right-hand side vector
    
    Returns:
    array: Solution vector for Toeplitz system
    """

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

# Generate AR signal with known order
np.random.seed(42)
N = 200
true_order = 4

# AR coefficients for a 4th-order system
ar_true = np.array([1.0, -2.7607, 3.8106, -2.6535, 0.9238])

# Generate AR signal
noise = np.random.randn(N)
signal = np.zeros(N)
signal[:true_order] = noise[:true_order]

for n in range(true_order, N):
    signal[n] = -np.sum(ar_true[1:] * signal[n-true_order:n][::-1]) + noise[n]

# Test different model orders
max_order = 15
orders = range(1, max_order + 1)
criteria_values = {
    'AIC': [],
    'MDL': [],
    'FPE': [],
    'CAT': [],
    'AICc': []
}

# Compute AR parameters and criteria for each order
for order in orders:
    # Estimate AR parameters using Yule-Walker
    ar_coeffs, rho, _ = spectrum.aryule(signal, order)
    
    # Compute various criteria
    criteria_values['AIC'].append(spectrum.AIC(N, rho, order))
    criteria_values['MDL'].append(spectrum.MDL(N, rho, order))
    criteria_values['FPE'].append(spectrum.FPE(N, rho, order))
    criteria_values['CAT'].append(spectrum.CAT(N, rho, order))
    criteria_values['AICc'].append(spectrum.AICc(N, rho, order))

# Find optimal orders
optimal_orders = {}
for name, values in criteria_values.items():
    optimal_orders[name] = orders[np.argmin(values)]

# Plot results
plt.figure(figsize=(15, 10))

plt.subplot(2, 3, 1)
plt.plot(signal)
plt.title(f'Generated AR({true_order}) Signal')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

# Plot criteria
for i, (name, values) in enumerate(criteria_values.items()):
    plt.subplot(2, 3, 2 + i)
    plt.plot(orders, values, 'o-', linewidth=2, markersize=6)
    plt.axvline(true_order, color='r', linestyle='--', alpha=0.7, label=f'True order = {true_order}')
    plt.axvline(optimal_orders[name], color='g', linestyle=':', alpha=0.7, 
               label=f'Optimal = {optimal_orders[name]}')
    plt.title(f'{name} Criterion')
    plt.xlabel('Model Order')
    plt.ylabel(f'{name} Value')
    plt.legend()
    plt.grid(True)

plt.tight_layout()
plt.show()

# Print results summary
print("Model Order Selection Results:")
print("-" * 40)
print(f"True AR order: {true_order}")
print("\nCriterion    Optimal Order")
print("-" * 25)
for name, optimal in optimal_orders.items():
    status = "✓" if optimal == true_order else "✗"
    print(f"{name:8s}     {optimal:2d}          {status}")

# Demonstrate Criteria class usage
print("\nUsing Criteria class:")
print("-" * 25)
aic_criterion = spectrum.Criteria('AIC', N)
for order in [3, 4, 5]:
    ar_coeffs, rho, _ = spectrum.aryule(signal, order)
    aic_val = aic_criterion(rho, order)  # This should work if implemented
    print(f"Order {order}: AIC = {aic_val:.3f}")

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

# Demonstrate Cholesky solver for linear systems
print("Cholesky Decomposition Example:")
print("-" * 35)

# Create a positive definite matrix (correlation matrix)
N = 100
n = np.arange(N)
signal = np.sin(2*np.pi*0.1*n) + 0.1*np.random.randn(N)

# Create autocorrelation matrix using different methods
methods = ['autocorrelation', 'covariance']
order = 10

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

for i, method in enumerate(methods):
    # Generate correlation matrix
    R = spectrum.corrmtx(signal, order-1, method=method)
    
    plt.subplot(2, 4, 1 + i)
    plt.imshow(np.real(R), cmap='coolwarm')
    plt.colorbar()
    plt.title(f'Correlation Matrix ({method})')
    
    plt.subplot(2, 4, 3 + i)
    plt.imshow(np.imag(R), cmap='coolwarm')
    plt.colorbar()
    plt.title(f'Imaginary Part ({method})')
    
    # Test Cholesky solver
    b = np.random.randn(R.shape[0])  # Random right-hand side
    
    try:
        # Solve using Cholesky decomposition
        x_chol = spectrum.CHOLESKY(R, b, method='scipy')
        
        # Solve using standard method for comparison
        x_std = np.linalg.solve(R, b)
        
        error = np.linalg.norm(x_chol - x_std)
        print(f"{method} method:")
        print(f"  Matrix size: {R.shape}")
        print(f"  Condition number: {np.linalg.cond(R):.2e}")
        print(f"  Cholesky error: {error:.2e}")
        print(f"  Matrix is positive definite: {np.all(np.linalg.eigvals(R) > 0)}")
        
    except Exception as e:
        print(f"Cholesky failed for {method}: {e}")

# Pascal matrix example
plt.subplot(2, 4, 5)
pascal_mat = spectrum.pascal(8)
plt.imshow(pascal_mat, cmap='Blues')
plt.colorbar()
plt.title('Pascal Matrix')

# Complex SVD demonstration
plt.subplot(2, 4, 6)
# Create complex test matrix
A_complex = np.random.randn(6, 4) + 1j * np.random.randn(6, 4)
U, s, Vh = spectrum.csvd(A_complex)

plt.semilogy(s, 'o-', linewidth=2, markersize=8)
plt.title('Singular Values (Complex SVD)')
plt.xlabel('Index')
plt.ylabel('Singular Value')
plt.grid(True)

plt.tight_layout()
plt.show()

# Verify SVD reconstruction
A_reconstructed = U @ np.diag(s) @ Vh
reconstruction_error = np.linalg.norm(A_complex - A_reconstructed)
print(f"\nComplex SVD verification:")
print(f"Original matrix shape: {A_complex.shape}")
print(f"Reconstruction error: {reconstruction_error:.2e}")

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

# Define a test transfer function (2nd order low-pass filter)
# H(s) = 1 / (s^2 + 2*ζ*ωn*s + ωn^2)
wn = 2*np.pi*10  # Natural frequency = 10 Hz
zeta = 0.7       # Damping ratio

# Continuous-time coefficients
num_s = [wn**2]
den_s = [1, 2*zeta*wn, wn**2]

# Convert to discrete-time using bilinear transform (approximation)
fs = 100  # Sampling frequency
T = 1/fs

# Bilinear transform: s = 2/T * (z-1)/(z+1)
# This is a simplified example - normally use scipy.signal.bilinear
num_z = num_s
den_z = den_s  # Simplified for demonstration

# Convert to zeros, poles, gain representation
zeros, poles, gain = spectrum.tf2zpk(num_z, den_z)

print("Transfer Function Analysis:")
print("-" * 30)
print(f"Numerator coefficients: {num_z}")
print(f"Denominator coefficients: {den_z}")
print(f"Zeros: {zeros}")
print(f"Poles: {poles}")
print(f"Gain: {gain}")

# Convert back to transfer function
num_recovered, den_recovered = spectrum.zpk2tf(zeros, poles, gain)
print(f"Recovered numerator: {num_recovered}")
print(f"Recovered denominator: {den_recovered}")

# Verify conversion accuracy
num_error = np.max(np.abs(np.array(num_z) - num_recovered[:len(num_z)]))
den_error = np.max(np.abs(np.array(den_z) - den_recovered[:len(den_z)]))
print(f"Numerator error: {num_error:.2e}")
print(f"Denominator error: {den_error:.2e}")

# Pole-zero plot and frequency response
plt.figure(figsize=(15, 8))

# Pole-zero plot
plt.subplot(2, 3, 1)
plt.plot(np.real(zeros), np.imag(zeros), 'bo', markersize=10, label='Zeros')
plt.plot(np.real(poles), np.imag(poles), 'rx', markersize=10, label='Poles')

# Unit circle for discrete-time systems
theta = np.linspace(0, 2*np.pi, 100)
plt.plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.5, label='Unit Circle')

plt.axis('equal')
plt.xlabel('Real Part')
plt.ylabel('Imaginary Part')
plt.title('Pole-Zero Plot')
plt.legend()
plt.grid(True)

# Frequency response
w = np.linspace(0, np.pi, 512)  # Normalized frequency
H = np.polyval(num_recovered, np.exp(1j*w)) / np.polyval(den_recovered, np.exp(1j*w))

plt.subplot(2, 3, 2)
plt.plot(w/(2*np.pi)*fs, 20*np.log10(np.abs(H)))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.title('Frequency Response - Magnitude')
plt.grid(True)

plt.subplot(2, 3, 3)
plt.plot(w/(2*np.pi)*fs, np.angle(H)*180/np.pi)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Phase (degrees)')
plt.title('Frequency Response - Phase')
plt.grid(True)

# Test lattice conversion for all-pole system (AR filter)
ar_coeffs = [1, -1.5, 0.7]  # AR polynomial
lattice_coeffs = spectrum.tf2latc(num=[1], den=ar_coeffs)

plt.subplot(2, 3, 4)
plt.stem(range(len(lattice_coeffs)), lattice_coeffs, basefmt=' ')
plt.title('Lattice Coefficients')
plt.xlabel('Stage')
plt.ylabel('Reflection Coefficient')
plt.grid(True)

# System stability analysis
plt.subplot(2, 3, 5)
# Plot poles in z-plane for stability analysis
pole_magnitudes = np.abs(poles)
stable = np.all(pole_magnitudes < 1)

colors = ['green' if mag < 1 else 'red' for mag in pole_magnitudes]
plt.scatter(np.real(poles), np.imag(poles), c=colors, s=100, alpha=0.7)
plt.plot(np.cos(theta), np.sin(theta), 'k--', alpha=0.5)
plt.axis('equal')
plt.xlim(-1.5, 1.5)
plt.ylim(-1.5, 1.5)
plt.xlabel('Real Part')
plt.ylabel('Imaginary Part')
plt.title(f'System Stability ({"Stable" if stable else "Unstable"})')
plt.grid(True)

# Group delay
plt.subplot(2, 3, 6)
# Approximate group delay calculation
w_gd = w[1:-1]
phase = np.angle(H)
group_delay = -np.diff(phase) / np.diff(w)
group_delay = np.append(group_delay, group_delay[-1])  # Extend for plotting

plt.plot(w_gd/(2*np.pi)*fs, group_delay[:-1])
plt.xlabel('Frequency (Hz)')
plt.ylabel('Group Delay (samples)')
plt.title('Group Delay')
plt.grid(True)

plt.tight_layout()
plt.show()

# Print stability analysis
print(f"\nStability Analysis:")
print(f"All poles inside unit circle: {stable}")
for i, (pole, mag) in enumerate(zip(poles, pole_magnitudes)):
    status = "Stable" if mag < 1 else "Unstable"
    print(f"Pole {i+1}: {pole:.3f}, |pole| = {mag:.3f} ({status})")

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

# Demonstrate Toeplitz solver
print("Toeplitz System Solver:")
print("-" * 25)

# Create Hermitian Toeplitz matrix
n = 8
T0 = 2.0  # Diagonal element
T = np.array([1.5, 0.8, 0.3, 0.1, 0.05, 0.02, 0.01])  # First row/column

# Construct full Toeplitz matrix for verification
from scipy.linalg import toeplitz
full_T = toeplitz([T0] + list(T), [T0] + list(T))

# Random right-hand side
Z = np.random.randn(n) + 1j * np.random.randn(n)

# Solve using specialized Toeplitz solver
try:
    X_toep = spectrum.HERMTOEP(T0, T, Z)
    print(f"Toeplitz solver successful")
    
    # Verify solution
    residual = np.linalg.norm(full_T @ X_toep - Z)
    print(f"Residual error: {residual:.2e}")
    
except Exception as e:
    print(f"Toeplitz solver failed: {e}")
    # Fall back to standard solver
    X_toep = np.linalg.solve(full_T, Z)

# Compare computational complexity demonstration
sizes = [8, 16, 32, 64, 128]
times_standard = []
times_toeplitz = []

import time

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

plt.subplot(2, 2, 1)
plt.imshow(np.real(full_T), cmap='coolwarm')
plt.colorbar()
plt.title('Toeplitz Matrix (Real Part)')

plt.subplot(2, 2, 2)
solution_error = np.abs(X_toep - np.linalg.solve(full_T, Z))
plt.semilogy(solution_error, 'o-')
plt.title('Solution Error (Toeplitz vs Standard)')
plt.xlabel('Element Index')
plt.ylabel('Absolute Error')
plt.grid(True)

# Model selection criteria comparison on real data
# Generate test signals with different characteristics
plt.subplot(2, 2, 3)

# Generate AR signals with different orders
true_orders = [2, 4, 6, 8]
N = 150
test_orders = range(1, 12)

criteria_comparison = {}
for criterion in ['AIC', 'MDL', 'FPE']:
    criteria_comparison[criterion] = []

for true_order in true_orders:
    # Generate random AR coefficients
    ar_coeffs = np.random.randn(true_order + 1)
    ar_coeffs[0] = 1  # Normalize
    
    # Ensure stability by scaling if needed
    poles = np.roots(ar_coeffs)
    if np.any(np.abs(poles) >= 1):
        ar_coeffs[1:] *= 0.9 / np.max(np.abs(poles))
    
    # Generate signal
    noise = np.random.randn(N)
    signal = np.zeros(N)
    signal[:true_order] = noise[:true_order]
    
    for n in range(true_order, N):
        signal[n] = -np.sum(ar_coeffs[1:] * signal[n-true_order:n][::-1]) + noise[n]
    
    # Find optimal orders using different criteria
    for criterion in criteria_comparison:
        criterion_values = []
        for order in test_orders:
            _, rho, _ = spectrum.aryule(signal, order)
            if criterion == 'AIC':
                val = spectrum.AIC(N, rho, order)
            elif criterion == 'MDL':
                val = spectrum.MDL(N, rho, order)
            elif criterion == 'FPE':
                val = spectrum.FPE(N, rho, order)
            criterion_values.append(val)
        
        optimal_order = test_orders[np.argmin(criterion_values)]
        criteria_comparison[criterion].append(optimal_order)

# Plot criterion performance
x_pos = np.arange(len(true_orders))
width = 0.25

for i, criterion in enumerate(criteria_comparison):
    plt.bar(x_pos + i*width, criteria_comparison[criterion], width, 
           label=criterion, alpha=0.7)

plt.bar(x_pos + len(criteria_comparison)*width, true_orders, width, 
       label='True Order', alpha=0.7, color='red')

plt.xlabel('Test Signal')
plt.ylabel('Selected Order')
plt.title('Criterion Performance Comparison')
plt.xticks(x_pos + width, [f'AR({order})' for order in true_orders])
plt.legend()
plt.grid(True, alpha=0.3)

# Statistical analysis of criterion performance
plt.subplot(2, 2, 4)
accuracy_scores = {}
for criterion in criteria_comparison:
    correct_selections = sum(1 for pred, true in zip(criteria_comparison[criterion], true_orders) 
                           if pred == true)
    accuracy_scores[criterion] = correct_selections / len(true_orders) * 100

criteria_names = list(accuracy_scores.keys())
accuracy_values = list(accuracy_scores.values())

plt.bar(criteria_names, accuracy_values, color=['blue', 'green', 'orange'])
plt.ylabel('Accuracy (%)')
plt.title('Order Selection Accuracy')
plt.ylim(0, 100)
for i, v in enumerate(accuracy_values):
    plt.text(i, v + 2, f'{v:.1f}%', ha='center')
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\nCriterion Performance Summary:")
print(f"-" * 35)
for criterion, accuracy in accuracy_scores.items():
    print(f"{criterion:8s}: {accuracy:5.1f}% correct selections")

# Demonstrate eigenvalue computations for correlation matrices
print(f"\nCorrelation Matrix Properties:")
print(f"-" * 32)

# Generate signal and correlation matrix
test_signal = np.random.randn(100)
R = spectrum.corrmtx(test_signal, 10, method='autocorrelation')

eigenvals = np.linalg.eigvals(R)
condition_number = np.max(eigenvals) / np.min(eigenvals)
trace = np.trace(R)
determinant = np.linalg.det(R)

print(f"Matrix size: {R.shape}")
print(f"Condition number: {condition_number:.2e}")
print(f"Trace: {trace:.3f}")
print(f"Determinant: {determinant:.3e}")
print(f"Positive definite: {np.all(eigenvals > 0)}")
print(f"Eigenvalue range: [{np.min(eigenvals):.3e}, {np.max(eigenvals):.3e}]")