Non-parametric Methods

def speriodogram(x, NFFT=None, detrend=True, sampling=1., scale_by_freq=True, window='hamming', axis=0):
    """
    Simple periodogram PSD estimation.
    
    Parameters:
    - x: array-like, input time series data
    - NFFT: int, FFT length (None uses len(x))
    - detrend: bool or str, detrending option (True, False, 'mean', 'linear')
    - sampling: float, sampling frequency
    - scale_by_freq: bool, scale PSD by frequency
    - window: str, window function name ('hamming', 'hanning', 'blackman', etc.)
    - axis: int, axis along which to compute PSD
    
    Returns:
    array: Power spectral density values
    """

class Periodogram(data, sampling=1., NFFT=None, window='hamming', detrend='mean', scale_by_freq=True):
    """
    Basic periodogram class for PSD estimation.
    
    Parameters:
    - data: array-like, input time series
    - sampling: float, sampling frequency
    - NFFT: int, FFT length for PSD computation
    - window: str, window function ('hanning', 'hamming', 'blackman', etc.)
    - detrend: str, detrending method ('mean', 'linear', None)
    - scale_by_freq: bool, scale PSD by frequency
    """

def WelchPeriodogram(data, NFFT=None, sampling=1., **kargs):
    """
    Welch's averaged periodogram method.
    
    Parameters:
    - data: array-like, input time series data
    - NFFT: int, segment length for averaging
    - sampling: float, sampling frequency
    - **kargs: additional parameters (window, overlap, detrend)
    
    Returns:
    tuple: (frequencies array, PSD values array)
    """

def DaniellPeriodogram(data, P, NFFT=None, detrend='mean', sampling=1., window_params={}):
    """
    Daniell's periodogram with smoothing.
    
    Parameters:
    - data: array-like, input time series data
    - P: int, smoothing parameter (Daniell kernel width)
    - NFFT: int, FFT length
    - detrend: str, detrending method
    - sampling: float, sampling frequency
    - window_params: dict, window function parameters
    
    Returns:
    tuple: (frequencies array, smoothed PSD array)
    """

class pdaniell(data, P, sampling=1., NFFT=None, detrend='mean', scale_by_freq=True):
    """
    Daniell periodogram class.
    
    Parameters:
    - data: array-like, input time series
    - P: int, Daniell smoothing parameter
    - sampling: float, sampling frequency
    - NFFT: int, FFT length for PSD computation
    - detrend: str, detrending method
    - scale_by_freq: bool, scale PSD by frequency
    """

class pcorrelogram(data, sampling=1., lag=-1, window='hamming', NFFT=None, scale_by_freq=True, detrend=None):
    """
    Correlogram PSD estimation class.
    
    Parameters:
    - data: array-like, input time series
    - sampling: float, sampling frequency
    - lag: int, maximum correlation lag
    - window: str, lag window function
    - NFFT: int, FFT length for PSD computation
    - scale_by_freq: bool, scale PSD by frequency
    - detrend: str, detrending method
    """

def pmtm(x, NW=None, k=None, NFFT=None, e=None, v=None, method="adapt", show=False):
    """
    Multi-taper method PSD estimation.
    
    Parameters:
    - x: array-like, input time series data
    - NW: float, time-bandwidth product (typically 2.5-4)
    - k: int, number of tapers (typically 2*NW - 1)
    - NFFT: int, FFT length
    - e: array-like, pre-computed tapers (DPSS sequences)
    - v: array-like, eigenvalues of tapers
    - method: str, averaging method ('adapt', 'unity', 'eigen')
    - show: bool, show intermediate plots
    
    Returns:
    tuple: (PSD array, frequencies array)
    """

def dpss(N, NW=None, k=None):
    """
    Discrete prolate spheroidal sequences (Slepian tapers).
    
    Parameters:
    - N: int, sequence length
    - NW: float, time-bandwidth product
    - k: int, number of sequences to compute
    
    Returns:
    tuple: (tapers array, eigenvalues array)
    """

class MultiTapering(Spectrum):
    """
    Multi-taper PSD estimation class.
    
    Parameters:
    - data: array-like, input time series
    - NW: float, time-bandwidth product
    - k: int, number of tapers
    - NFFT: int, FFT length
    - sampling: float, sampling frequency
    """

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

# Generate test signal with multiple components
fs = 1000  # Sampling frequency
T = 1.0    # Duration
t = np.linspace(0, T, int(fs*T), False)

# Signal: two sinusoids + noise
f1, f2 = 50, 120  # Hz
signal = (np.sin(2*np.pi*f1*t) + 0.5*np.sin(2*np.pi*f2*t) + 
         0.2*np.random.randn(len(t)))

# Different periodogram methods
p_basic = spectrum.Periodogram(signal, sampling=fs)
p_welch = spectrum.WelchPeriodogram(signal, NFFT=256, sampling=fs)
p_daniell = spectrum.pdaniell(signal, P=5, sampling=fs)

# Plot comparison
plt.figure(figsize=(12, 10))

plt.subplot(3, 1, 1)
plt.plot(t[:200], signal[:200])
plt.title('Input Signal (first 200 samples)')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

plt.subplot(3, 1, 2)
freqs = np.linspace(0, fs/2, len(p_basic.psd))
plt.semilogy(freqs, p_basic.psd, label='Basic Periodogram')
plt.axvline(f1, color='r', linestyle='--', alpha=0.7, label=f'f1={f1}Hz')
plt.axvline(f2, color='r', linestyle='--', alpha=0.7, label=f'f2={f2}Hz')
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.title('Basic Periodogram')
plt.legend()
plt.grid(True)
plt.xlim(0, 200)

plt.subplot(3, 1, 3)
freqs_welch, psd_welch = p_welch
plt.semilogy(freqs_welch, psd_welch, 'g-', label='Welch Method')
freqs_daniell = np.linspace(0, fs/2, len(p_daniell.psd))
plt.semilogy(freqs_daniell, p_daniell.psd, 'r-', label='Daniell Method')
plt.axvline(f1, color='k', linestyle='--', alpha=0.7)
plt.axvline(f2, color='k', linestyle='--', alpha=0.7)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.title('Welch vs Daniell Methods')
plt.legend()
plt.grid(True)
plt.xlim(0, 200)

plt.tight_layout()
plt.show()

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

# Generate signal with closely spaced frequencies
fs = 500
N = 500
t = np.arange(N) / fs

# Two closely spaced sinusoids + noise
f1, f2 = 100, 105  # 5 Hz apart
A1, A2 = 1.0, 0.8
noise_level = 0.3

signal = (A1*np.sin(2*np.pi*f1*t) + A2*np.sin(2*np.pi*f2*t) + 
         noise_level*np.random.randn(N))

# Multi-taper analysis with different parameters
NW_values = [2.5, 4.0, 6.0]
colors = ['blue', 'red', 'green']

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

# Plot signal
plt.subplot(2, 2, 1)
plt.plot(t, signal)
plt.title(f'Signal: {f1}Hz + {f2}Hz + noise')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

# Compare with basic periodogram
plt.subplot(2, 2, 2)
p_period = spectrum.Periodogram(signal, sampling=fs)
freqs = np.linspace(0, fs/2, len(p_period.psd))
plt.semilogy(freqs, p_period.psd, 'k-', alpha=0.7, label='Periodogram')
plt.axvline(f1, color='r', linestyle=':', alpha=0.8)
plt.axvline(f2, color='r', linestyle=':', alpha=0.8)
plt.xlim(90, 115)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.title('Basic Periodogram')
plt.legend()
plt.grid(True)

# Multi-taper with different NW values
plt.subplot(2, 2, 3)
for i, NW in enumerate(NW_values):
    k = int(2*NW - 1)  # Number of tapers
    psd_mt, freqs_mt = spectrum.pmtm(signal, NW=NW, k=k, NFFT=1024)
    plt.semilogy(freqs_mt*fs, psd_mt, color=colors[i], 
                label=f'NW={NW}, k={k}', linewidth=2)

plt.axvline(f1, color='r', linestyle=':', alpha=0.8)
plt.axvline(f2, color='r', linestyle=':', alpha=0.8)
plt.xlim(90, 115)
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.title('Multi-Taper Method (Different NW)')
plt.legend()
plt.grid(True)

# Taper visualization
plt.subplot(2, 2, 4)
NW = 4.0
k = int(2*NW - 1)
tapers, eigenvals = spectrum.dpss(N, NW=NW, k=k)
for i in range(min(4, k)):  # Show first 4 tapers
    plt.plot(tapers[:, i], label=f'Taper {i+1} (λ={eigenvals[i]:.3f})')
plt.xlabel('Sample')
plt.ylabel('Amplitude')  
plt.title(f'DPSS Tapers (N={N}, NW={NW})')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Print taper properties
print(f"Multi-taper parameters: NW={NW}, k={k}")
print(f"Frequency resolution: {2*NW/N*fs:.2f} Hz")
print(f"Eigenvalue concentration: {np.sum(eigenvals > 0.9)}/{k} tapers > 0.9")

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

# Generate AR signal for correlogram analysis
N = 256
# AR(2) coefficients creating spectral peak around 0.2-0.3 normalized freq
ar_coeffs = [1, -1.6, 0.81]  # Poles near unit circle
noise = np.random.randn(N)

# Generate AR signal
from scipy import signal
ar_signal = signal.lfilter([1], ar_coeffs, noise)

# Different correlation lags for correlogram
lag_values = [20, 40, 80]
windows = ['rectangular', 'hamming', 'blackman']

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

# Original signal
plt.subplot(2, 3, 1)
plt.plot(ar_signal)
plt.title('AR(2) Signal')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.grid(True)

# Autocorrelation function
plt.subplot(2, 3, 2)
max_lag = 50
autocorr = spectrum.CORRELATION(ar_signal, maxlags=max_lag, norm='unbiased')
lags = np.arange(-max_lag, max_lag+1)
plt.plot(lags, autocorr)
plt.title('Autocorrelation Function')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.grid(True)

# Different lag lengths
plt.subplot(2, 3, 4)
for i, lag in enumerate(lag_values):
    p_corr = spectrum.pcorrelogram(ar_signal, lag=lag, window='hamming')
    freqs = np.linspace(0, 0.5, len(p_corr.psd))
    plt.semilogy(freqs, p_corr.psd, label=f'Lag={lag}', linewidth=2)
plt.xlabel('Normalized Frequency')
plt.ylabel('PSD')
plt.title('Correlogram: Different Lag Lengths')
plt.legend()
plt.grid(True)

# Different windows
plt.subplot(2, 3, 5)
lag = 40
for i, window in enumerate(windows):
    p_corr = spectrum.pcorrelogram(ar_signal, lag=lag, window=window)
    freqs = np.linspace(0, 0.5, len(p_corr.psd))
    plt.semilogy(freqs, p_corr.psd, label=f'{window.capitalize()}', linewidth=2)
plt.xlabel('Normalized Frequency')  
plt.ylabel('PSD')
plt.title(f'Correlogram: Different Windows (lag={lag})')
plt.legend()
plt.grid(True)

# Comparison with parametric method
plt.subplot(2, 3, 6)
p_corr = spectrum.pcorrelogram(ar_signal, lag=40, window='hamming')
p_burg = spectrum.pburg(ar_signal, order=10)
freqs = np.linspace(0, 0.5, len(p_corr.psd))
plt.semilogy(freqs, p_corr.psd, 'b-', label='Correlogram', linewidth=2)
plt.semilogy(freqs, p_burg.psd, 'r--', label='Burg AR(10)', linewidth=2)
plt.xlabel('Normalized Frequency')
plt.ylabel('PSD')
plt.title('Correlogram vs Parametric')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

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

def compare_methods(signal, title, fs=1.0):
    """Compare different non-parametric methods on a given signal."""
    
    # Method parameters
    methods = {
        'Periodogram': spectrum.Periodogram(signal, sampling=fs),
        'Correlogram': spectrum.pcorrelogram(signal, lag=len(signal)//4, 
                                           window='hamming'),
        'Multi-taper': None  # Will compute separately
    }
    
    # Multi-taper method
    NW = 3.0
    k = int(2*NW - 1)
    psd_mt, freqs_mt = spectrum.pmtm(signal, NW=NW, k=k, NFFT=512)
    
    plt.figure(figsize=(12, 8))
    
    # Plot signal
    plt.subplot(2, 2, 1)
    t = np.arange(len(signal)) / fs
    plt.plot(t, signal)
    plt.title(f'{title} - Time Domain')
    plt.xlabel('Time (s)')
    plt.ylabel('Amplitude')
    plt.grid(True)
    
    # Plot spectra
    plt.subplot(2, 2, 2)
    
    # Periodogram
    freqs = np.linspace(0, fs/2, len(methods['Periodogram'].psd))
    plt.semilogy(freqs, methods['Periodogram'].psd, 'b-', 
                label='Periodogram', linewidth=2)
    
    # Correlogram  
    freqs_corr = np.linspace(0, fs/2, len(methods['Correlogram'].psd))
    plt.semilogy(freqs_corr, methods['Correlogram'].psd, 'g--', 
                label='Correlogram', linewidth=2)
    
    # Multi-taper
    plt.semilogy(freqs_mt*fs, psd_mt, 'r:', 
                label='Multi-taper', linewidth=2)
    
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('PSD')
    plt.title(f'{title} - PSD Comparison')
    plt.legend()
    plt.grid(True)
    
    plt.tight_layout()
    plt.show()

# Test different signal types
fs = 100

# 1. Narrow-band signal
N1 = 300
t1 = np.arange(N1) / fs
signal1 = np.sin(2*np.pi*10*t1) + 0.3*np.random.randn(N1)
compare_methods(signal1, "Narrow-band Signal (10 Hz)", fs)

# 2. Wide-band signal  
N2 = 400
signal2 = np.random.randn(N2)  # White noise
# Add some color with filtering
from scipy import signal as sp_signal
b, a = sp_signal.butter(4, 0.3, 'low')
signal2 = sp_signal.lfilter(b, a, signal2)
compare_methods(signal2, "Wide-band Colored Noise", fs)

# 3. Multi-component signal
N3 = 500
t3 = np.arange(N3) / fs  
signal3 = (np.sin(2*np.pi*5*t3) + 0.8*np.sin(2*np.pi*15*t3) + 
          0.6*np.sin(2*np.pi*25*t3) + 0.2*np.random.randn(N3))
compare_methods(signal3, "Multi-component Signal", fs)