"""Noise analysis for LTspice .noise simulation results.

LTspice .noise analysis produces output with variables like 'onoise'
(output-referred noise spectral density in V/sqrt(Hz)) and 'inoise'
(input-referred noise spectral density). The data is complex-valued
in the .raw file; magnitude gives the spectral density.
"""

import numpy as np

# np.trapz was renamed to np.trapezoid in numpy 2.0
_trapz = getattr(np, "trapezoid", getattr(np, "trapz", None))

# Boltzmann constant (J/K)
_K_BOLTZMANN = 1.380649e-23


def compute_noise_spectral_density(frequency: np.ndarray, noise_signal: np.ndarray) -> dict:
    """Compute noise spectral density from raw noise simulation data.

    Takes the frequency array and complex noise signal directly from the
    .raw file and returns the noise spectral density in V/sqrt(Hz) and dB.

    Args:
        frequency: Frequency array in Hz (may be complex; real part is used)
        noise_signal: Complex noise signal from .raw file (magnitude = V/sqrt(Hz))

    Returns:
        Dict with frequency_hz, noise_density_v_per_sqrt_hz, noise_density_db
    """
    if len(frequency) == 0 or len(noise_signal) == 0:
        return {
            "frequency_hz": [],
            "noise_density_v_per_sqrt_hz": [],
            "noise_density_db": [],
        }

    freq = np.real(frequency).astype(np.float64)
    density = np.abs(noise_signal)

    # Single-point case: still return valid data
    density_db = 20.0 * np.log10(np.maximum(density, 1e-30))

    return {
        "frequency_hz": freq.tolist(),
        "noise_density_v_per_sqrt_hz": density.tolist(),
        "noise_density_db": density_db.tolist(),
    }


def compute_total_noise(
    frequency: np.ndarray,
    noise_signal: np.ndarray,
    f_low: float | None = None,
    f_high: float | None = None,
) -> dict:
    """Integrate noise spectral density over frequency to get total RMS noise.

    Computes total_rms = sqrt(integral(|noise|^2 * df)) using trapezoidal
    integration over the specified frequency range.

    Args:
        frequency: Frequency array in Hz (may be complex; real part is used)
        noise_signal: Complex noise signal from .raw file
        f_low: Lower integration bound in Hz (default: min frequency in data)
        f_high: Upper integration bound in Hz (default: max frequency in data)

    Returns:
        Dict with total_rms_v, integration_range_hz, equivalent_noise_bandwidth_hz
    """
    if len(frequency) < 2 or len(noise_signal) < 2:
        return {
            "total_rms_v": 0.0,
            "integration_range_hz": [0.0, 0.0],
            "equivalent_noise_bandwidth_hz": 0.0,
        }

    freq = np.real(frequency).astype(np.float64)
    density = np.abs(noise_signal)

    # Sort by frequency to ensure correct integration order
    sort_idx = np.argsort(freq)
    freq = freq[sort_idx]
    density = density[sort_idx]

    # Apply frequency bounds
    if f_low is None:
        f_low = float(freq[0])
    if f_high is None:
        f_high = float(freq[-1])

    mask = (freq >= f_low) & (freq <= f_high)
    freq_band = freq[mask]
    density_band = density[mask]

    if len(freq_band) < 2:
        return {
            "total_rms_v": 0.0,
            "integration_range_hz": [f_low, f_high],
            "equivalent_noise_bandwidth_hz": 0.0,
        }

    # Integrate |noise|^2 over frequency, then take sqrt for RMS
    noise_power = density_band**2
    integrated = float(_trapz(noise_power, freq_band))
    total_rms = float(np.sqrt(max(integrated, 0.0)))

    # Equivalent noise bandwidth: bandwidth of a brick-wall filter with the
    # same peak density that would pass the same total noise power
    peak_density = float(np.max(density_band))
    if peak_density > 1e-30:
        enbw = integrated / (peak_density**2)
    else:
        enbw = 0.0

    return {
        "total_rms_v": total_rms,
        "integration_range_hz": [f_low, f_high],
        "equivalent_noise_bandwidth_hz": float(enbw),
    }


def compute_spot_noise(frequency: np.ndarray, noise_signal: np.ndarray, target_freq: float) -> dict:
    """Interpolate noise spectral density at a specific frequency.

    Uses linear interpolation between adjacent data points to estimate
    the noise density at the requested frequency.

    Args:
        frequency: Frequency array in Hz (may be complex; real part is used)
        noise_signal: Complex noise signal from .raw file
        target_freq: Desired frequency in Hz

    Returns:
        Dict with spot_noise_v_per_sqrt_hz, spot_noise_db, actual_freq_hz
    """
    if len(frequency) == 0 or len(noise_signal) == 0:
        return {
            "spot_noise_v_per_sqrt_hz": 0.0,
            "spot_noise_db": float("-inf"),
            "actual_freq_hz": target_freq,
        }

    freq = np.real(frequency).astype(np.float64)
    density = np.abs(noise_signal)

    # Sort by frequency
    sort_idx = np.argsort(freq)
    freq = freq[sort_idx]
    density = density[sort_idx]

    # Clamp to data range
    if target_freq <= freq[0]:
        spot = float(density[0])
        actual = float(freq[0])
    elif target_freq >= freq[-1]:
        spot = float(density[-1])
        actual = float(freq[-1])
    else:
        # Linear interpolation
        spot = float(np.interp(target_freq, freq, density))
        actual = target_freq

    spot_db = 20.0 * np.log10(max(spot, 1e-30))

    return {
        "spot_noise_v_per_sqrt_hz": spot,
        "spot_noise_db": spot_db,
        "actual_freq_hz": actual,
    }


def compute_noise_figure(
    frequency: np.ndarray,
    noise_signal: np.ndarray,
    source_resistance: float = 50.0,
    temperature: float = 290.0,
) -> dict:
    """Compute noise figure from noise spectral density.

    Noise figure is the ratio of the measured output noise power to
    the thermal noise of the source resistance at the given temperature.
    NF(f) = 10*log10(|noise(f)|^2 / (4*k*T*R))

    Args:
        frequency: Frequency array in Hz (may be complex; real part is used)
        noise_signal: Complex noise signal from .raw file (output-referred)
        source_resistance: Source impedance in ohms (default 50)
        temperature: Temperature in Kelvin (default 290 K, IEEE standard)

    Returns:
        Dict with noise_figure_db (array), frequency_hz (array),
        min_nf_db, nf_at_1khz
    """
    if len(frequency) == 0 or len(noise_signal) == 0:
        return {
            "noise_figure_db": [],
            "frequency_hz": [],
            "min_nf_db": None,
            "nf_at_1khz": None,
        }

    freq = np.real(frequency).astype(np.float64)
    density = np.abs(noise_signal)

    # Thermal noise power spectral density of the source: 4*k*T*R (V^2/Hz)
    thermal_psd = 4.0 * _K_BOLTZMANN * temperature * source_resistance

    if thermal_psd < 1e-50:
        return {
            "noise_figure_db": [],
            "frequency_hz": freq.tolist(),
            "min_nf_db": None,
            "nf_at_1khz": None,
        }

    # NF = 10*log10(measured_noise_power / thermal_noise_power)
    # where noise_power = density^2 per Hz
    noise_power = density**2
    nf_ratio = noise_power / thermal_psd
    nf_db = 10.0 * np.log10(np.maximum(nf_ratio, 1e-30))

    min_nf_db = float(np.min(nf_db))

    # Noise figure at 1 kHz (interpolated)
    nf_at_1khz = None
    if freq[0] <= 1000.0 <= freq[-1]:
        sort_idx = np.argsort(freq)
        nf_at_1khz = float(np.interp(1000.0, freq[sort_idx], nf_db[sort_idx]))
    elif len(freq) == 1:
        nf_at_1khz = float(nf_db[0])

    return {
        "noise_figure_db": nf_db.tolist(),
        "frequency_hz": freq.tolist(),
        "min_nf_db": min_nf_db,
        "nf_at_1khz": nf_at_1khz,
    }


def _estimate_flicker_corner(frequency: np.ndarray, density: np.ndarray) -> float | None:
    """Estimate the 1/f noise corner frequency.

    The 1/f corner is where the noise transitions from 1/f (flicker) behavior
    to flat (white) noise. We find where the slope of log(density) vs log(freq)
    crosses -0.25 (midpoint between 0 for white and -0.5 for 1/f in V/sqrt(Hz)).

    Args:
        frequency: Sorted frequency array in Hz (positive, ascending)
        density: Noise spectral density magnitude (same order as frequency)

    Returns:
        Corner frequency in Hz, or None if not detectable
    """
    if len(frequency) < 4:
        return None

    # Work in log-log space
    pos_mask = (frequency > 0) & (density > 0)
    freq_pos = frequency[pos_mask]
    dens_pos = density[pos_mask]

    if len(freq_pos) < 4:
        return None

    log_f = np.log10(freq_pos)
    log_d = np.log10(dens_pos)

    # Compute local slope using central differences (smoothed)
    # Use a window of ~5 points for robustness
    n = len(log_f)
    slopes = np.zeros(n)
    half_win = min(2, (n - 1) // 2)

    for i in range(half_win, n - half_win):
        df = log_f[i + half_win] - log_f[i - half_win]
        dd = log_d[i + half_win] - log_d[i - half_win]
        if abs(df) > 1e-15:
            slopes[i] = dd / df

    # Fill edges with nearest valid slope
    slopes[:half_win] = slopes[half_win]
    slopes[n - half_win :] = slopes[n - half_win - 1]

    # Find where slope crosses the threshold (-0.25)
    # 1/f noise has slope ~ -0.5 in V/sqrt(Hz), white has slope ~ 0
    threshold = -0.25

    for i in range(len(slopes) - 1):
        if slopes[i] < threshold <= slopes[i + 1]:
            # Interpolate
            ds = slopes[i + 1] - slopes[i]
            if abs(ds) < 1e-15:
                return float(freq_pos[i])
            frac = (threshold - slopes[i]) / ds
            log_corner = log_f[i] + frac * (log_f[i + 1] - log_f[i])
            return float(10.0**log_corner)

    return None


def compute_noise_metrics(
    frequency: np.ndarray,
    noise_signal: np.ndarray,
    source_resistance: float = 50.0,
) -> dict:
    """Comprehensive noise analysis report.

    Combines spectral density, spot noise at standard frequencies, total
    integrated noise, noise figure, and 1/f corner estimation.

    Args:
        frequency: Frequency array in Hz (may be complex; real part is used)
        noise_signal: Complex noise signal from .raw file
        source_resistance: Source impedance in ohms for noise figure (default 50)

    Returns:
        Dict with spectral_density, spot_noise (at standard frequencies),
        total_noise, noise_figure, flicker_corner_hz
    """
    if len(frequency) < 2 or len(noise_signal) < 2:
        return {
            "spectral_density": compute_noise_spectral_density(frequency, noise_signal),
            "spot_noise": {},
            "total_noise": compute_total_noise(frequency, noise_signal),
            "noise_figure": compute_noise_figure(frequency, noise_signal, source_resistance),
            "flicker_corner_hz": None,
        }

    freq = np.real(frequency).astype(np.float64)

    # Spectral density
    spectral = compute_noise_spectral_density(frequency, noise_signal)

    # Spot noise at standard frequencies
    spot_freqs = [10.0, 100.0, 1000.0, 10000.0, 100000.0]
    spot_labels = ["10Hz", "100Hz", "1kHz", "10kHz", "100kHz"]
    spot_noise = {}

    f_min = float(np.min(freq))
    f_max = float(np.max(freq))

    for label, sf in zip(spot_labels, spot_freqs):
        if f_min <= sf <= f_max:
            spot_noise[label] = compute_spot_noise(frequency, noise_signal, sf)

    # Total noise over full bandwidth
    total = compute_total_noise(frequency, noise_signal)

    # Noise figure
    nf = compute_noise_figure(frequency, noise_signal, source_resistance)

    # 1/f corner frequency estimation
    sort_idx = np.argsort(freq)
    sorted_freq = freq[sort_idx]
    sorted_density = np.abs(noise_signal)[sort_idx]
    flicker_corner = _estimate_flicker_corner(sorted_freq, sorted_density)

    return {
        "spectral_density": spectral,
        "spot_noise": spot_noise,
        "total_noise": total,
        "noise_figure": nf,
        "flicker_corner_hz": float(flicker_corner) if flicker_corner is not None else None,
    }