///|
/// Spectral utilities for RR-derived tachograms.

///|
/// Named power-band definition.
pub(all) struct SpectralBand {
  name : String
  lower_hz : Double
  upper_hz : Double
} derive(FromJson, ToJson, Debug, Eq)

///|
/// Power and normalized power for one band.
pub(all) struct BandPower {
  name : String
  lower_hz : Double
  upper_hz : Double
  power : Double
  normalized_power : Double
  contribution_percent : Double
} derive(FromJson, ToJson, Debug, Eq)

///|
/// Spectral profile over a set of bands.
pub(all) struct SpectralProfile {
  bands : Array[BandPower]
  total_power : Double
  dominant_band : String
  centroid_hz : Double
  entropy : Double
  edge_95_hz : Double
  slope : Double
} derive(FromJson, ToJson, Debug, Eq)

///|
/// Standard frequency bands used by short resting recordings.
pub fn standard_spectral_bands() -> Array[SpectralBand] {
  [
    { name: "ulf", lower_hz: 0.0, upper_hz: 0.0033 },
    { name: "vlf", lower_hz: 0.0033, upper_hz: 0.04 },
    { name: "lf", lower_hz: 0.04, upper_hz: 0.15 },
    { name: "hf", lower_hz: 0.15, upper_hz: 0.40 },
    { name: "respiratory", lower_hz: 0.10, upper_hz: 0.50 },
  ]
}

///|
/// Clamp and normalize a user-provided spectral band.
pub fn normalize_spectral_band(band : SpectralBand) -> SpectralBand {
  let lower = if band.lower_hz < 0.0 { 0.0 } else { band.lower_hz }
  let upper = if band.upper_hz < lower { lower } else { band.upper_hz }
  { name: band.name, lower_hz: lower, upper_hz: upper }
}

///|
/// Calculate power in a band from an already-computed periodogram.
pub fn spectral_band_power(
  spectrum : Array[SpectrumBin],
  band : SpectralBand,
) -> Double {
  let normalized = normalize_spectral_band(band)
  let mut power = 0.0
  for bin in spectrum {
    if bin.frequency_hz >= normalized.lower_hz &&
      bin.frequency_hz < normalized.upper_hz {
      power += bin.power
    }
  }
  power
}

///|
/// Calculate all named band powers and relative contributions.
pub fn calculate_band_powers(
  spectrum : Array[SpectrumBin],
  bands : Array[SpectralBand],
) -> Array[BandPower] {
  let total = total_spectral_power(spectrum)
  let result = []
  for band in bands {
    let normalized = normalize_spectral_band(band)
    let power = spectral_band_power(spectrum, normalized)
    result.push({
      name: normalized.name,
      lower_hz: normalized.lower_hz,
      upper_hz: normalized.upper_hz,
      power,
      normalized_power: if total == 0.0 {
        0.0
      } else {
        power / total
      },
      contribution_percent: if total == 0.0 {
        0.0
      } else {
        power / total * 100.0
      },
    })
  }
  result
}

///|
/// Find the band with the largest absolute power.
pub fn dominant_band_name(bands : Array[BandPower]) -> String {
  let mut best = ""
  let mut power = -1.0
  for band in bands {
    if band.power > power {
      power = band.power
      best = band.name
    }
  }
  best
}

///|
/// Calculate a band power ratio with a safe zero denominator.
pub fn band_power_ratio(
  bands : Array[BandPower],
  numerator : String,
  denominator : String,
) -> Double {
  let mut top = 0.0
  let mut bottom = 0.0
  for band in bands {
    if band.name == numerator {
      top = band.power
    }
    if band.name == denominator {
      bottom = band.power
    }
  }
  if bottom == 0.0 {
    0.0
  } else {
    top / bottom
  }
}

///|
/// Return the frequency at which a cumulative power fraction is reached.
pub fn cumulative_power_frequency(
  spectrum : Array[SpectrumBin],
  fraction : Double,
) -> Double {
  if spectrum.length() == 0 {
    return 0.0
  }
  let target = fraction.clamp(min=0.0, max=1.0) * total_spectral_power(spectrum)
  let mut cumulative = 0.0
  for bin in spectrum {
    cumulative += bin.power
    if cumulative >= target {
      return bin.frequency_hz
    }
  }
  spectrum[spectrum.length() - 1].frequency_hz
}

///|
/// Calculate a simple log-log spectral slope.
pub fn spectral_log_slope(spectrum : Array[SpectrumBin]) -> Double {
  let points = []
  for bin in spectrum {
    if bin.frequency_hz > 0.0 && bin.power > 0.0 {
      points.push((@math.ln(bin.frequency_hz), @math.ln(bin.power)))
    }
  }
  if points.length() <= 1 {
    return 0.0
  }
  let mean_x = sum_values(points.map(fn(point) { point.0 })) /
    points.length().to_double()
  let mean_y = sum_values(points.map(fn(point) { point.1 })) /
    points.length().to_double()
  let mut numerator = 0.0
  let mut denominator = 0.0
  for point in points {
    let dx = point.0 - mean_x
    numerator += dx * (point.1 - mean_y)
    denominator += dx * dx
  }
  if denominator == 0.0 {
    0.0
  } else {
    numerator / denominator
  }
}

///|
/// Calculate the spectral profile for a tachogram.
pub fn calculate_spectral_profile(
  intervals : Array[Double],
  sample_rate_hz : Double,
  bands : Array[SpectralBand],
) -> SpectralProfile {
  let spectrum = calculate_periodogram(intervals, sample_rate_hz)
  let band_powers = calculate_band_powers(spectrum, bands)
  {
    bands: band_powers,
    total_power: total_spectral_power(spectrum),
    dominant_band: dominant_band_name(band_powers),
    centroid_hz: spectral_centroid(spectrum),
    entropy: spectral_entropy(spectrum),
    edge_95_hz: cumulative_power_frequency(spectrum, 0.95),
    slope: spectral_log_slope(spectrum),
  }
}

///|
/// Return a stable spectral feature vector.
pub fn spectral_profile_feature_vector(
  profile : SpectralProfile,
) -> Array[Double] {
  let result = []
  for band in profile.bands {
    result.push(band.normalized_power)
  }
  result.push(profile.total_power)
  result.push(profile.centroid_hz)
  result.push(profile.entropy)
  result.push(profile.edge_95_hz)
  result.push(profile.slope)
  result.push(band_power_ratio(profile.bands, "lf", "hf"))
  result
}

///|
/// Interpolate a spectrum onto an evenly-spaced frequency grid.
pub fn interpolate_spectrum(
  spectrum : Array[SpectrumBin],
  step_hz : Double,
  maximum_hz : Double,
) -> Array[SpectrumBin] {
  if spectrum.length() == 0 || step_hz <= 0.0 || maximum_hz <= 0.0 {
    return []
  }
  let result = []
  let mut frequency = 0.0
  while frequency <= maximum_hz {
    let mut nearest = spectrum[0]
    let mut distance = absolute_difference(nearest.frequency_hz, frequency)
    for bin in spectrum {
      let candidate_distance = absolute_difference(bin.frequency_hz, frequency)
      if candidate_distance < distance {
        distance = candidate_distance
        nearest = bin
      }
    }
    result.push({
      frequency_hz: frequency,
      power: nearest.power,
      amplitude: nearest.amplitude,
    })
    frequency += step_hz
  }
  result
}

///|
/// Smooth spectral power with a moving average.
pub fn smooth_spectrum(
  spectrum : Array[SpectrumBin],
  radius : Int,
) -> Array[SpectrumBin] {
  if radius < 0 {
    return []
  }
  let result = []
  for i in 0.. spectrum.length() {
      spectrum.length()
    } else {
      i + radius + 1
    }
    let values = []
    for j in start.. Bool {
  let bound = if minimum < 0.0 { 0.0 } else { minimum }
  spectral_band_power(spectrum, band) >= bound
}

///|
/// Estimate a respiratory peak from the respiratory band.
pub fn respiratory_peak_frequency(
  intervals : Array[Double],
  sample_rate_hz : Double,
) -> Double {
  let spectrum = calculate_periodogram(intervals, sample_rate_hz)
  let respiratory = { name: "respiratory", lower_hz: 0.10, upper_hz: 0.50 }
  let selected = []
  for bin in spectrum {
    if bin.frequency_hz >= respiratory.lower_hz &&
      bin.frequency_hz <= respiratory.upper_hz {
      selected.push(bin)
    }
  }
  if selected.length() == 0 {
    0.0
  } else {
    dominant_spectrum_bin(selected).frequency_hz
  }
}

///|
/// Convert a respiratory frequency to breaths per minute.
pub fn respiratory_rate_bpm(frequency_hz : Double) -> Double {
  if frequency_hz < 0.0 {
    0.0
  } else {
    frequency_hz * 60.0
  }
}

///|
/// Return whether the spectral profile contains valid finite values.
pub fn spectral_profile_is_usable(profile : SpectralProfile) -> Bool {
  profile.total_power > 0.0 &&
  profile.centroid_hz >= 0.0 &&
  profile.entropy >= 0.0 &&
  profile.edge_95_hz >= 0.0
}