Correlation and Signal Processing

def CORRELATION(x, y=None, maxlags=None, norm='unbiased'):
    """
    Correlation function computation (auto or cross-correlation).
    
    Parameters:
    - x: array-like, first input signal
    - y: array-like, second input signal (None for autocorrelation)
    - maxlags: int, maximum lag to compute (None for full range)
    - norm: str, normalization type ('biased', 'unbiased', 'coeff', None)
    
    Returns:
    array: Correlation values for positive lags only
    """

def xcorr(x, y=None, maxlags=None, norm='biased'):
    """
    Fast cross-correlation using numpy.correlate.
    
    Parameters:
    - x: array-like, first input signal
    - y: array-like, second input signal (None for autocorrelation)
    - maxlags: int, maximum lag to compute
    - norm: str, normalization type ('biased', 'unbiased', 'coeff', None)
    
    Returns:
    tuple: (correlation sequence array, lags array)
    """

def chirp(t, f0=0., t1=1., f1=100., form='linear', phase=0):
    """
    Generate chirp (frequency-swept) signal.
    
    Parameters:
    - t: array-like, time vector
    - f0: float, starting frequency
    - t1: float, end time for frequency sweep
    - f1: float, ending frequency
    - form: str, sweep type ('linear', 'quadratic', 'logarithmic', 'hyperbolic')
    - phase: float, initial phase in radians
    
    Returns:
    array: Generated chirp signal
    """

def morlet(lb, ub, n):
    """
    Generate Morlet wavelet for time-frequency analysis.
    
    Parameters:
    - lb: float, lower bound of time interval
    - ub: float, upper bound of time interval
    - n: int, number of samples
    
    Returns:
    array: Morlet wavelet samples
    """

def mexican(lb, ub, n):
    """
    Generate Mexican hat wavelet (second derivative of Gaussian).
    
    Parameters:
    - lb: float, lower bound of time interval
    - ub: float, upper bound of time interval  
    - n: int, number of samples
    
    Returns:
    array: Mexican hat wavelet samples
    """

def meyeraux(x):
    """
    Meyer wavelet auxiliary function.
    
    Parameters:
    - x: array-like, input values
    
    Returns:
    array: Meyer auxiliary function values
    """

def pow2db(x):
    """
    Convert power values to decibels.
    
    Parameters:
    - x: array-like, power values (linear scale)
    
    Returns:
    array: Power in decibels (10*log10(x))
    """

def db2pow(xdb):
    """
    Convert decibel values to linear power.
    
    Parameters:
    - xdb: array-like, power values in dB
    
    Returns:
    array: Linear power values (10^(xdb/10))
    """

def mag2db(x):
    """
    Convert magnitude values to decibels.
    
    Parameters:
    - x: array-like, magnitude values
    
    Returns:
    array: Magnitude in decibels (20*log10(x))
    """

def db2mag(xdb):
    """
    Convert decibel values to linear magnitude.
    
    Parameters:
    - xdb: array-like, magnitude values in dB
    
    Returns:
    array: Linear magnitude values (10^(xdb/20))
    """

def nextpow2(x):
    """
    Find the next power of 2 greater than or equal to x.
    
    Parameters:
    - x: float or int, input value
    
    Returns:
    int: Next power of 2
    """

def fftshift(x):
    """
    Wrapper to numpy.fft.fftshift for frequency domain centering.
    
    Parameters:
    - x: array-like, input array
    
    Returns:
    array: Shifted array with zero frequency at center
    """

def cshift(data, offset):
    """
    Circular shift of array elements.
    
    Parameters:
    - data: array-like, input data array
    - offset: int, number of positions to shift
    
    Returns:
    array: Circularly shifted array
    """

def twosided_2_onesided(data):
    """
    Convert two-sided PSD to one-sided PSD.
    
    Parameters:
    - data: array-like, two-sided PSD values
    
    Returns:
    array: One-sided PSD values (positive frequencies only)
    """

def onesided_2_twosided(data):
    """
    Convert one-sided PSD to two-sided PSD.
    
    Parameters:
    - data: array-like, one-sided PSD values
    
    Returns:
    array: Two-sided PSD values (negative and positive frequencies)
    """

def twosided_2_centerdc(data):
    """
    Convert two-sided PSD to center-DC format.
    
    Parameters:
    - data: array-like, two-sided PSD values
    
    Returns:
    array: Center-DC format PSD (DC at center)
    """

def centerdc_2_twosided(data):
    """
    Convert center-DC format PSD to standard two-sided format.
    
    Parameters:
    - data: array-like, center-DC format PSD
    
    Returns:
    array: Two-sided PSD values
    """

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

# Generate test signals for correlation analysis
N = 200
n = np.arange(N)

# Signal 1: Sinusoid + noise
f1 = 0.1
signal1 = np.sin(2*np.pi*f1*n) + 0.3*np.random.randn(N)

# Signal 2: Delayed and noisy version of signal1
delay = 20
signal2 = np.zeros(N)
signal2[delay:] = signal1[:-delay] + 0.2*np.random.randn(N-delay)

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

# Plot signals
plt.subplot(3, 2, 1)
plt.plot(signal1, 'b-', label='Signal 1', alpha=0.8)
plt.plot(signal2, 'r-', label='Signal 2 (delayed)', alpha=0.8)
plt.title('Input Signals')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)

# Autocorrelation of signal 1
plt.subplot(3, 2, 2)
maxlags = 50
autocorr1 = spectrum.CORRELATION(signal1, maxlags=maxlags, norm='unbiased')
lags_pos = np.arange(maxlags+1)
plt.plot(lags_pos, autocorr1, 'b-', linewidth=2)
plt.title('Autocorrelation of Signal 1')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.grid(True)

# Cross-correlation using CORRELATION function
plt.subplot(3, 2, 3)
crosscorr = spectrum.CORRELATION(signal1, signal2, maxlags=maxlags, norm='unbiased')
plt.plot(lags_pos, crosscorr, 'g-', linewidth=2)
plt.title('Cross-correlation (CORRELATION function)')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.grid(True)

# Cross-correlation using xcorr function (returns both positive and negative lags)
plt.subplot(3, 2, 4)
crosscorr_full, lags_full = spectrum.xcorr(signal1, signal2, maxlags=maxlags, norm='unbiased')
plt.plot(lags_full, crosscorr_full, 'r-', linewidth=2)
plt.axvline(delay, color='k', linestyle='--', alpha=0.7, label=f'True delay={delay}')
peak_lag = lags_full[np.argmax(crosscorr_full)]
plt.axvline(peak_lag, color='g', linestyle='--', alpha=0.7, label=f'Detected delay={peak_lag}')
plt.title('Cross-correlation (xcorr function)')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.legend()
plt.grid(True)

# Compare different normalization methods
plt.subplot(3, 2, 5)
norms = ['biased', 'unbiased', 'coeff', None]
for norm_type in norms:
    if norm_type is not None:
        corr_norm = spectrum.CORRELATION(signal1, maxlags=30, norm=norm_type)
        label = f'{norm_type}'
    else:
        corr_norm = spectrum.CORRELATION(signal1, maxlags=30, norm=None)
        label = 'unnormalized'
    
    lags_norm = np.arange(len(corr_norm))
    plt.plot(lags_norm, corr_norm, linewidth=2, label=label)

plt.title('Autocorrelation: Different Normalizations')
plt.xlabel('Lag')
plt.ylabel('Correlation')
plt.legend()
plt.grid(True)

# Practical application: Time delay estimation
plt.subplot(3, 2, 6)
delays_to_test = range(5, 35, 5)
detected_delays = []
correlation_peaks = []

for test_delay in delays_to_test:
    # Create delayed version
    sig2_test = np.zeros(N)
    sig2_test[test_delay:] = signal1[:-test_delay] + 0.1*np.random.randn(N-test_delay)
    
    # Find peak in cross-correlation
    cc, lags = spectrum.xcorr(signal1, sig2_test, maxlags=maxlags)
    peak_idx = np.argmax(cc)
    detected_delay = lags[peak_idx]
    peak_value = cc[peak_idx]
    
    detected_delays.append(detected_delay)
    correlation_peaks.append(peak_value)

plt.plot(delays_to_test, detected_delays, 'bo-', linewidth=2, label='Detected')
plt.plot(delays_to_test, delays_to_test, 'r--', linewidth=2, label='True delay')
plt.title('Time Delay Estimation Accuracy')
plt.xlabel('True Delay (samples)')
plt.ylabel('Detected Delay (samples)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Print delay estimation results
print("Delay Estimation Results:")
print("-" * 30)
for true_del, det_del, peak in zip(delays_to_test, detected_delays, correlation_peaks):
    error = abs(true_del - det_del)
    print(f"True: {true_del:2d}, Detected: {det_del:2d}, Error: {error:2d}, Peak: {peak:.3f}")

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

# Time parameters
fs = 1000  # Sampling frequency
T = 2.0    # Duration
t = np.linspace(0, T, int(fs*T), False)

# Generate different types of signals
plt.figure(figsize=(15, 12))

# 1. Linear chirp
plt.subplot(3, 3, 1)
chirp_linear = spectrum.chirp(t, f0=50, t1=T, f1=200, form='linear')
plt.plot(t, chirp_linear)
plt.title('Linear Chirp (50-200 Hz)')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

# 2. Logarithmic chirp
plt.subplot(3, 3, 2)
chirp_log = spectrum.chirp(t, f0=50, t1=T, f1=200, form='logarithmic')
plt.plot(t, chirp_log)
plt.title('Logarithmic Chirp (50-200 Hz)')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

# 3. Quadratic chirp
plt.subplot(3, 3, 3)
chirp_quad = spectrum.chirp(t, f0=50, t1=T, f1=200, form='quadratic')
plt.plot(t, chirp_quad)
plt.title('Quadratic Chirp (50-200 Hz)')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

# 4. Morlet wavelet
plt.subplot(3, 3, 4)
morlet_wave = spectrum.morlet(-3, 3, 200)
t_morlet = np.linspace(-3, 3, 200)
plt.plot(t_morlet, morlet_wave, 'r-', linewidth=2)
plt.title('Morlet Wavelet')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.grid(True)

# 5. Mexican hat wavelet
plt.subplot(3, 3, 5)
mexican_wave = spectrum.mexican(-3, 3, 200)
plt.plot(t_morlet, mexican_wave, 'g-', linewidth=2)
plt.title('Mexican Hat Wavelet')
plt.xlabel('Time')
plt.ylabel('Amplitude')
plt.grid(True)

# Spectrograms of chirp signals
plt.subplot(3, 3, 7)
from scipy import signal
f_spec, t_spec, Sxx = signal.spectrogram(chirp_linear, fs, nperseg=256, noverlap=128)
plt.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx), shading='gouraud')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
plt.title('Linear Chirp Spectrogram')
plt.ylim(0, 300)

plt.subplot(3, 3, 8)
f_spec, t_spec, Sxx = signal.spectrogram(chirp_log, fs, nperseg=256, noverlap=128)
plt.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx), shading='gouraud')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
plt.title('Logarithmic Chirp Spectrogram')
plt.ylim(0, 300)

plt.subplot(3, 3, 9)
f_spec, t_spec, Sxx = signal.spectrogram(chirp_quad, fs, nperseg=256, noverlap=128)
plt.pcolormesh(t_spec, f_spec, 10*np.log10(Sxx), shading='gouraud')
plt.ylabel('Frequency (Hz)')
plt.xlabel('Time (s)')
plt.title('Quadratic Chirp Spectrogram')
plt.ylim(0, 300)

plt.tight_layout()
plt.show()

# Compare chirp types in frequency domain
plt.figure(figsize=(12, 8))

signals = {
    'Linear': chirp_linear,
    'Logarithmic': chirp_log, 
    'Quadratic': chirp_quad
}

for i, (chirp_type, chirp_signal) in enumerate(signals.items()):
    plt.subplot(2, 2, i+1)
    
    # Compute PSD using Welch method
    f_psd, psd = signal.welch(chirp_signal, fs, nperseg=1024, noverlap=512)
    plt.semilogy(f_psd, psd, linewidth=2)
    plt.title(f'{chirp_type} Chirp - PSD')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('PSD')
    plt.grid(True)
    plt.xlim(0, 300)

# Wavelet frequency content
plt.subplot(2, 2, 4)
morlet_fft = np.fft.fft(morlet_wave, 1024)
freqs_wav = np.fft.fftfreq(1024, d=6/200)  # Adjust for time span
plt.semilogy(np.fft.fftshift(freqs_wav), np.fft.fftshift(np.abs(morlet_fft)**2), 'r-', linewidth=2, label='Morlet')

mexican_fft = np.fft.fft(mexican_wave, 1024)
plt.semilogy(np.fft.fftshift(freqs_wav), np.fft.fftshift(np.abs(mexican_fft)**2), 'g-', linewidth=2, label='Mexican Hat')
plt.title('Wavelet Frequency Content')
plt.xlabel('Frequency')
plt.ylabel('Power')
plt.legend()
plt.grid(True)
plt.xlim(-2, 2)

plt.tight_layout()
plt.show()

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

# Generate test data with different dynamic ranges
linear_powers = np.logspace(-3, 3, 100)  # 10^-3 to 10^3
linear_magnitudes = np.sqrt(linear_powers)

# Convert to dB scales
power_db = spectrum.pow2db(linear_powers)
magnitude_db = spectrum.mag2db(linear_magnitudes)

# Convert back to linear
power_recovered = spectrum.db2pow(power_db)
magnitude_recovered = spectrum.db2mag(magnitude_db)

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

# Linear vs dB power scales
plt.subplot(2, 3, 1)
plt.loglog(linear_powers, linear_powers, 'b-', linewidth=2, label='Linear scale')
plt.xlabel('Input Power')
plt.ylabel('Power')
plt.title('Linear Power Scale')
plt.grid(True)
plt.legend()

plt.subplot(2, 3, 2)
plt.semilogx(linear_powers, power_db, 'r-', linewidth=2)
plt.xlabel('Linear Power')
plt.ylabel('Power (dB)')
plt.title('Power to dB Conversion')
plt.grid(True)

plt.subplot(2, 3, 3)
plt.plot(power_db, power_recovered, 'g-', linewidth=2)
plt.xlabel('Power (dB)')
plt.ylabel('Recovered Linear Power')
plt.title('dB to Power Conversion')
plt.grid(True)

# Magnitude conversions
plt.subplot(2, 3, 4)
plt.loglog(linear_magnitudes, linear_magnitudes, 'b-', linewidth=2, label='Linear scale')
plt.xlabel('Input Magnitude')
plt.ylabel('Magnitude')
plt.title('Linear Magnitude Scale')
plt.grid(True)
plt.legend()

plt.subplot(2, 3, 5)
plt.semilogx(linear_magnitudes, magnitude_db, 'r-', linewidth=2)
plt.xlabel('Linear Magnitude')
plt.ylabel('Magnitude (dB)')
plt.title('Magnitude to dB Conversion')
plt.grid(True)

plt.subplot(2, 3, 6)
plt.plot(magnitude_db, magnitude_recovered, 'g-', linewidth=2)
plt.xlabel('Magnitude (dB)')
plt.ylabel('Recovered Linear Magnitude')
plt.title('dB to Magnitude Conversion')
plt.grid(True)

plt.tight_layout()
plt.show()

# Demonstrate relationship between power and magnitude dB
print("Power vs Magnitude dB Relationship:")
print("-" * 40)
test_values = [0.1, 1.0, 10.0, 100.0]

for val in test_values:
    pow_db = spectrum.pow2db(val)
    mag_db = spectrum.mag2db(np.sqrt(val))  # magnitude = sqrt(power)
    print(f"Power: {val:6.1f} -> {pow_db:6.1f} dB")
    print(f"  Mag: {np.sqrt(val):6.3f} -> {mag_db:6.1f} dB")
    print(f"  Ratio: {pow_db/mag_db:.1f} (should be 2.0)")
    print()

# Practical example: Signal analysis with dB scaling
fs = 1000
t = np.linspace(0, 1, fs)

# Multi-tone signal with different amplitudes
signal = (2.0 * np.sin(2*np.pi*50*t) +    # Strong component
         0.5 * np.sin(2*np.pi*150*t) +    # Medium component  
         0.1 * np.sin(2*np.pi*300*t) +    # Weak component
         0.02 * np.random.randn(fs))       # Noise

# Compute PSD and convert to dB
f, psd_linear = signal.welch(signal, fs, nperseg=256)
psd_db = spectrum.pow2db(psd_linear)

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

plt.subplot(2, 2, 1)
plt.plot(t[:200], signal[:200])
plt.title('Multi-tone Signal')
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

plt.subplot(2, 2, 2)
plt.plot(f, psd_linear)
plt.title('PSD - Linear Scale')
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.grid(True)

plt.subplot(2, 2, 3)
plt.semilogy(f, psd_linear)
plt.title('PSD - Log Scale') 
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD')
plt.grid(True)

plt.subplot(2, 2, 4)
plt.plot(f, psd_db)
plt.title('PSD - dB Scale')
plt.xlabel('Frequency (Hz)')
plt.ylabel('PSD (dB)')
plt.grid(True)

plt.tight_layout()
plt.show()

# Dynamic range comparison
linear_range = np.max(psd_linear) / np.min(psd_linear[psd_linear > 0])
db_range = np.max(psd_db) - np.min(psd_db[np.isfinite(psd_db)])

print(f"Dynamic Range Analysis:")
print(f"Linear scale: {linear_range:.1e} (ratio)")
print(f"dB scale: {db_range:.1f} dB")
print(f"Equivalent ratio: {spectrum.db2pow(db_range):.1e}")

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

# Demonstrate utility functions
print("Utility Functions Demonstration:")
print("-" * 35)

# Next power of 2
test_values = [100, 255, 256, 513, 1000, 1024]
for val in test_values:
    next_pow2 = spectrum.nextpow2(val)
    print(f"nextpow2({val:4d}) = {next_pow2:4d} (2^{int(np.log2(next_pow2))})")

# Generate test PSD for format conversions
N = 512
freqs = np.linspace(0, 0.5, N//2 + 1)  # One-sided frequencies
# Create artificial PSD with multiple peaks
psd_onesided = (np.exp(-(freqs-0.1)**2/0.01) + 
               0.5*np.exp(-(freqs-0.3)**2/0.005) + 
               0.1*np.ones_like(freqs))  # Noise floor

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

# Original one-sided PSD
plt.subplot(2, 3, 1)
plt.plot(freqs, psd_onesided, 'b-', linewidth=2)
plt.title('Original One-sided PSD')
plt.xlabel('Normalized Frequency')
plt.ylabel('PSD')
plt.grid(True)

# Convert to two-sided
psd_twosided = spectrum.onesided_2_twosided(psd_onesided)
freqs_twosided = np.linspace(-0.5, 0.5, len(psd_twosided))

plt.subplot(2, 3, 2)
plt.plot(freqs_twosided, psd_twosided, 'r-', linewidth=2)
plt.title('Two-sided PSD')
plt.xlabel('Normalized Frequency')  
plt.ylabel('PSD')
plt.grid(True)

# Convert two-sided to center-DC
psd_centerdc = spectrum.twosided_2_centerdc(psd_twosided)

plt.subplot(2, 3, 3)
plt.plot(freqs_twosided, psd_centerdc, 'g-', linewidth=2)
plt.title('Center-DC Format PSD')
plt.xlabel('Normalized Frequency')
plt.ylabel('PSD')
plt.grid(True)

# Convert back to two-sided
psd_twosided_recovered = spectrum.centerdc_2_twosided(psd_centerdc)

plt.subplot(2, 3, 4)
plt.plot(freqs_twosided, psd_twosided_recovered, 'm--', linewidth=2, label='Recovered')
plt.plot(freqs_twosided, psd_twosided, 'r:', alpha=0.7, label='Original')
plt.title('Recovered Two-sided PSD')
plt.xlabel('Normalized Frequency')
plt.ylabel('PSD')
plt.legend()
plt.grid(True)

# Convert back to one-sided
psd_onesided_recovered = spectrum.twosided_2_onesided(psd_twosided_recovered)

plt.subplot(2, 3, 5)
plt.plot(freqs, psd_onesided_recovered, 'c--', linewidth=2, label='Recovered')
plt.plot(freqs, psd_onesided, 'b:', alpha=0.7, label='Original')
plt.title('Recovered One-sided PSD')
plt.xlabel('Normalized Frequency')
plt.ylabel('PSD')
plt.legend()
plt.grid(True)

# Circular shift demonstration
plt.subplot(2, 3, 6)
test_signal = np.sin(2*np.pi*0.1*np.arange(100))
shifts = [0, 10, -15, 25]
for i, shift in enumerate(shifts):
    shifted = spectrum.cshift(test_signal, shift)
    plt.plot(shifted[:50], label=f'Shift {shift}', alpha=0.8)
plt.title('Circular Shift Examples')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

# Verify conversion accuracy
print("\nFormat Conversion Accuracy:")
print("-" * 30)
error_onesided = np.max(np.abs(psd_onesided - psd_onesided_recovered))
error_twosided = np.max(np.abs(psd_twosided - psd_twosided_recovered))

print(f"One-sided round-trip error: {error_onesided:.2e}")
print(f"Two-sided round-trip error: {error_twosided:.2e}")

# Demonstrate FFT shift
fft_example = np.fft.fft(np.random.randn(16))
print("\nFFT Shift Example:")
print("Original FFT order:", np.arange(16))
print("After fftshift:   ", np.arange(16)[np.fft.fftshift(np.arange(16)) == np.arange(16)])

# Show frequency ordering
N_fft = 16
original_freqs = np.fft.fftfreq(N_fft)
shifted_freqs = spectrum.fftshift(original_freqs)

print(f"Original freq bins: {original_freqs}")
print(f"Shifted freq bins:  {shifted_freqs}")