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

Eigenvalue Methods

Name: tessl/pypi-spectrum
Author: tessl

High-resolution spectral estimation using eigenvalue decomposition of the data covariance matrix. These methods excel at resolving closely spaced frequency components and provide superior performance for short data records, particularly in the presence of noise.

Capabilities

MUSIC Method

Multiple Signal Classification (MUSIC) method provides high-resolution pseudospectrum estimation by exploiting the orthogonality between signal and noise subspaces.

def music(X, IP, NSIG=None, NFFT=4096, threshold=None, criteria='aic', verbose=False):
    """
    MUSIC pseudospectrum estimation.
    
    Parameters:
    - X: array-like, input time series data
    - IP: int, total number of sinusoids + noise sources
    - NSIG: int, number of sinusoidal sources (auto-detected if None)
    - NFFT: int, FFT length for pseudospectrum computation
    - threshold: float, threshold for automatic signal detection
    - criteria: str, model order selection criteria ('aic', 'mdl')
    - verbose: bool, enable verbose output
    
    Returns:
    array: MUSIC pseudospectrum values
    """

class pmusic(data, IP, NSIG=None, NFFT=None, sampling=1., scale_by_freq=False):
    """
    MUSIC pseudospectrum class.
    
    Parameters:
    - data: array-like, input time series
    - IP: int, total number of complex sinusoids + noise
    - NSIG: int, number of complex sinusoids (auto if None)
    - NFFT: int, FFT length for pseudospectrum computation
    - sampling: float, sampling frequency
    - scale_by_freq: bool, scale pseudospectrum by frequency
    """

Eigenvector Method

Eigenvector method provides an alternative eigenvalue-based approach with different weighting of the noise subspace eigenvectors.

def ev(X, IP, NSIG=None, NFFT=4096, threshold=None, criteria='aic', verbose=False):
    """
    Eigenvector method pseudospectrum estimation.
    
    Parameters:
    - X: array-like, input time series data
    - IP: int, total number of sinusoids + noise sources
    - NSIG: int, number of sinusoidal sources (auto-detected if None)
    - NFFT: int, FFT length for pseudospectrum computation
    - threshold: float, threshold for automatic signal detection
    - criteria: str, model order selection criteria ('aic', 'mdl')
    - verbose: bool, enable verbose output
    
    Returns:
    array: Eigenvector pseudospectrum values
    """

class pev(data, IP, NSIG=None, NFFT=None, sampling=1., scale_by_freq=False):
    """
    Eigenvector pseudospectrum class.
    
    Parameters:
    - data: array-like, input time series
    - IP: int, total number of complex sinusoids + noise  
    - NSIG: int, number of complex sinusoids (auto if None)
    - NFFT: int, FFT length for pseudospectrum computation
    - sampling: float, sampling frequency
    - scale_by_freq: bool, scale pseudospectrum by frequency
    """

Minimum Variance Method

Minimum variance distortionless response method provides adaptive beamforming-based spectral estimation.

def minvar(X, order, sampling=1., NFFT=4096):
    """
    Minimum variance pseudospectrum estimation.
    
    Parameters:
    - X: array-like, input time series data
    - order: int, order of the minimum variance filter
    - sampling: float, sampling frequency
    - NFFT: int, FFT length for pseudospectrum computation
    
    Returns:
    array: Minimum variance pseudospectrum values
    """

class pminvar(data, order, NFFT=None, sampling=1., scale_by_freq=False):
    """
    Minimum variance pseudospectrum class.
    
    Parameters:
    - data: array-like, input time series
    - order: int, order of the minimum variance filter
    - NFFT: int, FFT length for pseudospectrum computation
    - sampling: float, sampling frequency
    - scale_by_freq: bool, scale pseudospectrum by frequency
    """

General Eigenvalue Method

Unified interface for eigenvalue-based methods allowing selection between MUSIC and eigenvector methods.

def eigen(X, P, NSIG=None, method='music', threshold=None, NFFT=4096, criteria='aic', verbose=False):
    """
    General eigenvalue-based pseudospectrum estimation.
    
    Parameters:
    - X: array-like, input time series data
    - P: int, total number of sinusoids + noise sources
    - NSIG: int, number of sinusoidal sources (auto if None)
    - method: str, method selection ('music', 'ev')
    - threshold: float, threshold for signal detection
    - NFFT: int, FFT length for pseudospectrum
    - criteria: str, order selection criteria ('aic', 'mdl')
    - verbose: bool, enable verbose output
    
    Returns:
    array: Pseudospectrum values
    """

Model Order Selection

Automatic determination of signal and noise subspace dimensions using information-theoretic criteria.

def aic_eigen(s, N):
    """
    AIC criterion for eigenvalue methods.
    
    Parameters:
    - s: array-like, eigenvalues in descending order
    - N: int, data length
    
    Returns:
    array: AIC values for different model orders
    """

def mdl_eigen(s, N):
    """
    MDL criterion for eigenvalue methods.
    
    Parameters:
    - s: array-like, eigenvalues in descending order
    - N: int, data length
    
    Returns:
    array: MDL values for different model orders
    """

Usage Examples

MUSIC Spectral Analysis

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

# Generate signal with two closely spaced sinusoids in noise
N = 64
n = np.arange(N)
f1, f2 = 0.25, 0.27  # Closely spaced normalized frequencies
SNR_dB = 10

# Create signal
signal = (np.exp(1j * 2 * np.pi * f1 * n) + 
         np.exp(1j * 2 * np.pi * f2 * n))
noise_power = 10**(-SNR_dB/10) * np.mean(np.abs(signal)**2)
noise = np.sqrt(noise_power/2) * (np.random.randn(N) + 1j*np.random.randn(N))
noisy_signal = signal + noise

# MUSIC analysis
IP = 15  # Model order (number of complex exponentials + noise)
NSIG = 2  # Number of signal components

# Automatic signal number detection
p_music_auto = spectrum.pmusic(noisy_signal, IP=IP)
print(f"Auto-detected signals: {p_music_auto.NSIG}")

# Manual signal specification  
p_music = spectrum.pmusic(noisy_signal, IP=IP, NSIG=NSIG)

# Compare with eigenvector method
p_ev = spectrum.pev(noisy_signal, IP=IP, NSIG=NSIG)

# Plot results
freqs = np.linspace(0, 1, len(p_music.psd))
plt.figure(figsize=(12, 8))

plt.subplot(2, 1, 1)
plt.plot(n, np.real(noisy_signal), 'b-', alpha=0.7, label='Real part')
plt.plot(n, np.imag(noisy_signal), 'r-', alpha=0.7, label='Imag part')
plt.title(f'Noisy Signal: Two Complex Sinusoids (f1={f1}, f2={f2})')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(freqs, 10*np.log10(p_music.psd), 'b-', linewidth=2, label='MUSIC')
plt.plot(freqs, 10*np.log10(p_ev.psd), 'r--', linewidth=2, label='Eigenvector')
plt.axvline(f1, color='k', linestyle=':', alpha=0.7, label='True f1')
plt.axvline(f2, color='k', linestyle=':', alpha=0.7, label='True f2')
plt.xlabel('Normalized Frequency')
plt.ylabel('Pseudospectrum (dB)')
plt.title('MUSIC vs Eigenvector Method')
plt.legend()
plt.grid(True)
plt.xlim(0.2, 0.3)

plt.tight_layout()
plt.show()

Model Order Selection

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

# Generate test signal with known number of components
N = 100
n = np.arange(N)
frequencies = [0.1, 0.15, 0.3]  # Three sinusoids
amplitudes = [1.0, 0.8, 0.6]

signal = sum(A * np.exp(1j * 2 * np.pi * f * n) 
            for A, f in zip(amplitudes, frequencies))
noise = 0.1 * (np.random.randn(N) + 1j * np.random.randn(N))
noisy_signal = signal + noise

# Test different model orders
orders = range(5, 20)
aic_values = []
mdl_values = []

for order in orders:
    # Compute correlation matrix  
    R = spectrum.corrmtx(noisy_signal, order-1, method='autocorrelation')
    eigenvals = np.linalg.eigvals(R)
    eigenvals = np.sort(eigenvals)[::-1]  # Descending order
    
    # Compute criteria
    aic_vals = spectrum.aic_eigen(eigenvals, N)
    mdl_vals = spectrum.mdl_eigen(eigenvals, N)
    
    aic_values.append(np.min(aic_vals))
    mdl_values.append(np.min(mdl_vals))

# Plot model order selection
plt.figure(figsize=(10, 6))
plt.subplot(1, 2, 1)
plt.plot(orders, aic_values, 'bo-', label='AIC')
plt.xlabel('Model Order')
plt.ylabel('AIC Value')
plt.title('AIC Model Order Selection')
plt.grid(True)
optimal_aic = orders[np.argmin(aic_values)]
plt.axvline(optimal_aic, color='r', linestyle='--', 
           label=f'Optimal: {optimal_aic}')
plt.legend()

plt.subplot(1, 2, 2)  
plt.plot(orders, mdl_values, 'ro-', label='MDL')
plt.xlabel('Model Order')
plt.ylabel('MDL Value')
plt.title('MDL Model Order Selection')
plt.grid(True)
optimal_mdl = orders[np.argmin(mdl_values)]
plt.axvline(optimal_mdl, color='b', linestyle='--', 
           label=f'Optimal: {optimal_mdl}')
plt.legend()

plt.tight_layout()
plt.show()

print(f"True number of signals: {len(frequencies)}")
print(f"AIC optimal order: {optimal_aic}")
print(f"MDL optimal order: {optimal_mdl}")

Minimum Variance Spectral Analysis

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

# Generate signal with multiple frequency components
N = 128
t = np.arange(N)
fs = 1.0  # Sampling frequency

# Multiple sinusoids with different amplitudes
signal = (2.0 * np.sin(2*np.pi*0.1*t) +   # Strong component
         1.0 * np.sin(2*np.pi*0.25*t) +    # Medium component  
         0.5 * np.sin(2*np.pi*0.4*t))      # Weak component
noise = 0.2 * np.random.randn(N)
noisy_signal = signal + noise

# Compare different eigenvalue methods
orders = [10, 15, 20]
methods = ['MUSIC', 'Eigenvector', 'Minimum Variance']

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

for i, order in enumerate(orders):
    plt.subplot(len(orders), 3, 3*i + 1)
    # MUSIC
    p_music = spectrum.pmusic(noisy_signal, IP=order, NSIG=3)
    freqs = np.linspace(0, fs/2, len(p_music.psd))
    plt.semilogy(freqs, p_music.psd)
    plt.title(f'MUSIC (Order={order})')
    plt.ylabel('Pseudospectrum')
    if i == len(orders)-1:
        plt.xlabel('Frequency')
    plt.grid(True)
    
    plt.subplot(len(orders), 3, 3*i + 2)
    # Eigenvector  
    p_ev = spectrum.pev(noisy_signal, IP=order, NSIG=3)
    plt.semilogy(freqs, p_ev.psd)
    plt.title(f'Eigenvector (Order={order})')
    if i == len(orders)-1:
        plt.xlabel('Frequency')
    plt.grid(True)
    
    plt.subplot(len(orders), 3, 3*i + 3)
    # Minimum Variance
    p_minvar = spectrum.pminvar(noisy_signal, order=order)
    plt.semilogy(freqs, p_minvar.psd)
    plt.title(f'Min Variance (Order={order})')
    if i == len(orders)-1:
        plt.xlabel('Frequency')
    plt.grid(True)

plt.tight_layout()
plt.show()

Install with Tessl CLI