///|
/// A frequency interval used by spectral quality and feature calculations.
pub struct SpectralBand {
  low : Double
  high : Double
}

///|
/// Construct a normalized spectral band.
pub fn SpectralBand::new(low : Double, high : Double) -> SpectralBand {
  let left = if low < 0.0 { 0.0 } else { low }
  let right = if high < left { left } else { high }
  { low: left, high: right }
}

///|
/// Return whether a frequency belongs to a band.
pub fn spectral_band_contains(band : SpectralBand, frequency : Double) -> Bool {
  frequency >= band.low && frequency <= band.high
}

///|
/// A compact set of features used by event classifiers and reports.
pub struct AdvancedSpectralFeatures {
  centroid : Double
  bandwidth : Double
  rolloff : Double
  flatness : Double
  entropy : Double
  skewness : Double
  kurtosis : Double
  peak_frequency : Double
  peak_power : Double
  slope : Double
}

///|
/// Sum non-negative powers in a selected frequency band.
pub fn advanced_band_power(spectrum : Spectrum, band : SpectralBand) -> Double {
  let mut total = 0.0
  let count = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < count; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) {
      total += spectrum.powers[index].max(0.0)
    }
  }
  total
}

///|
/// Return total positive spectral power.
pub fn advanced_total_power(spectrum : Spectrum) -> Double {
  advanced_band_power(spectrum, SpectralBand::new(0.0, 1.0e300))
}

///|
/// Compute the power-weighted spectral centroid.
pub fn advanced_spectral_centroid(
  spectrum : Spectrum,
  band : SpectralBand,
) -> Double {
  let mut weighted = 0.0
  let mut total = 0.0
  let count = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < count; index = index + 1 {
    let frequency = spectrum.frequencies[index]
    let power = spectrum.powers[index].max(0.0)
    if spectral_band_contains(band, frequency) {
      weighted += frequency * power
      total += power
    }
  }
  if total <= 0.0 {
    0.0
  } else {
    weighted / total
  }
}

///|
/// Compute the power-weighted spectral bandwidth.
pub fn advanced_spectral_bandwidth(
  spectrum : Spectrum,
  band : SpectralBand,
) -> Double {
  let centroid = advanced_spectral_centroid(spectrum, band)
  let mut weighted = 0.0
  let mut total = 0.0
  let count = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < count; index = index + 1 {
    let frequency = spectrum.frequencies[index]
    let power = spectrum.powers[index].max(0.0)
    if spectral_band_contains(band, frequency) {
      let distance = frequency - centroid
      weighted += distance * distance * power
      total += power
    }
  }
  if total <= 0.0 {
    0.0
  } else {
    (weighted / total).sqrt()
  }
}

///|
/// Find the frequency at which a power fraction is reached.
pub fn advanced_spectral_rolloff(
  spectrum : Spectrum,
  band : SpectralBand,
  fraction : Double,
) -> Double {
  let target_fraction = fraction.max(0.0).min(1.0)
  let total = advanced_band_power(spectrum, band)
  if total <= 0.0 {
    return band.low
  }
  let target = total * target_fraction
  let mut cumulative = 0.0
  let count = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < count; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) {
      cumulative += spectrum.powers[index].max(0.0)
      if cumulative >= target {
        return spectrum.frequencies[index]
      }
    }
  }
  band.high
}

///|
/// Compute the geometric-to-arithmetic mean ratio of power.
pub fn advanced_spectral_flatness(
  spectrum : Spectrum,
  band : SpectralBand,
) -> Double {
  let mut logarithmic = 0.0
  let mut arithmetic = 0.0
  let mut count = 0
  let limit = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < limit; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) {
      let power = spectrum.powers[index].max(1.0e-20)
      logarithmic += @math.ln(power)
      arithmetic += power
      count += 1
    }
  }
  if count == 0 || arithmetic <= 0.0 {
    0.0
  } else {
    let geometric = @math.exp(logarithmic / count.to_double())
    geometric / (arithmetic / count.to_double())
  }
}

///|
/// Compute Shannon entropy of normalized spectral power.
pub fn advanced_spectral_entropy(
  spectrum : Spectrum,
  band : SpectralBand,
) -> Double {
  let total = advanced_band_power(spectrum, band)
  if total <= 0.0 {
    return 0.0
  }
  let mut entropy = 0.0
  let limit = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < limit; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) {
      let probability = spectrum.powers[index].max(0.0) / total
      if probability > 0.0 {
        entropy -= probability * @math.ln(probability)
      }
    }
  }
  entropy
}

///|
/// Compute the third standardized central moment of power distribution.
pub fn advanced_spectral_skewness(
  spectrum : Spectrum,
  band : SpectralBand,
) -> Double {
  let mean = advanced_spectral_centroid(spectrum, band)
  let width = advanced_spectral_bandwidth(spectrum, band)
  let total = advanced_band_power(spectrum, band)
  if width <= 0.0 || total <= 0.0 {
    return 0.0
  }
  let mut moment = 0.0
  let limit = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < limit; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) {
      let z = (spectrum.frequencies[index] - mean) / width
      moment += z * z * z * spectrum.powers[index].max(0.0) / total
    }
  }
  moment
}

///|
/// Compute the fourth standardized central moment minus three.
pub fn advanced_spectral_kurtosis(
  spectrum : Spectrum,
  band : SpectralBand,
) -> Double {
  let mean = advanced_spectral_centroid(spectrum, band)
  let width = advanced_spectral_bandwidth(spectrum, band)
  let total = advanced_band_power(spectrum, band)
  if width <= 0.0 || total <= 0.0 {
    return 0.0
  }
  let mut moment = 0.0
  let limit = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < limit; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) {
      let z = (spectrum.frequencies[index] - mean) / width
      moment += z * z * z * z * spectrum.powers[index].max(0.0) / total
    }
  }
  moment - 3.0
}

///|
/// Return the strongest frequency and power in a band.
pub fn advanced_spectral_peak(
  spectrum : Spectrum,
  band : SpectralBand,
) -> (Double, Double) {
  let mut frequency = band.low
  let mut power = 0.0
  let limit = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < limit; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) &&
      spectrum.powers[index] > power {
      frequency = spectrum.frequencies[index]
      power = spectrum.powers[index]
    }
  }
  (frequency, power)
}

///|
/// Estimate a least-squares slope of log power against frequency.
pub fn advanced_spectral_slope(
  spectrum : Spectrum,
  band : SpectralBand,
) -> Double {
  let mut count = 0
  let mut sum_x = 0.0
  let mut sum_y = 0.0
  let mut sum_xx = 0.0
  let mut sum_xy = 0.0
  let limit = if spectrum.frequencies.length() < spectrum.powers.length() {
    spectrum.frequencies.length()
  } else {
    spectrum.powers.length()
  }
  for index = 0; index < limit; index = index + 1 {
    if spectral_band_contains(band, spectrum.frequencies[index]) {
      let x = spectrum.frequencies[index]
      let y = @math.ln(spectrum.powers[index].max(1.0e-20))
      count += 1
      sum_x += x
      sum_y += y
      sum_xx += x * x
      sum_xy += x * y
    }
  }
  let denominator = count.to_double() * sum_xx - sum_x * sum_x
  if count < 2 || denominator.abs() < 1.0e-20 {
    0.0
  } else {
    (count.to_double() * sum_xy - sum_x * sum_y) / denominator
  }
}

///|
/// Compare two spectra using a normalized Euclidean distance.
pub fn advanced_spectral_distance(left : Spectrum, right : Spectrum) -> Double {
  let count = if left.powers.length() < right.powers.length() {
    left.powers.length()
  } else {
    right.powers.length()
  }
  if count == 0 {
    return 0.0
  }
  let left_total = advanced_total_power(left).max(1.0e-20)
  let right_total = advanced_total_power(right).max(1.0e-20)
  let mut sum = 0.0
  for index = 0; index < count; index = index + 1 {
    let a = left.powers[index] / left_total
    let b = right.powers[index] / right_total
    let difference = a - b
    sum += difference * difference
  }
  (sum / count.to_double()).sqrt()
}

///|
/// Smooth power bins with a centered moving average.
pub fn smooth_spectral_power(
  spectrum : Spectrum,
  radius : Int,
) -> Array[Double] {
  let result = Array::make(spectrum.powers.length(), 0.0)
  let window = if radius < 0 { 0 } else { radius }
  for index = 0; index < spectrum.powers.length(); index = index + 1 {
    let start = if index < window { 0 } else { index - window }
    let end = if index + window + 1 > spectrum.powers.length() {
      spectrum.powers.length()
    } else {
      index + window + 1
    }
    let mut sum = 0.0
    for cursor = start; cursor < end; cursor = cursor + 1 {
      sum += spectrum.powers[cursor]
    }
    result[index] = if end == start {
      0.0
    } else {
      sum / (end - start).to_double()
    }
  }
  result
}

///|
/// Normalize a power vector to unit sum.
pub fn normalize_spectral_power(power : Array[Double]) -> Array[Double] {
  let result = power.copy()
  let mut total = 0.0
  for value in result {
    total += value.max(0.0)
  }
  if total <= 0.0 {
    return Array::make(result.length(), 0.0)
  }
  for index = 0; index < result.length(); index = index + 1 {
    result[index] = result[index].max(0.0) / total
  }
  result
}

///|
/// Calculate spectral flux between adjacent spectra.
pub fn advanced_spectral_flux(
  previous : Spectrum,
  current : Spectrum,
) -> Double {
  let left = normalize_spectral_power(previous.powers)
  let right = normalize_spectral_power(current.powers)
  let count = if left.length() < right.length() {
    left.length()
  } else {
    right.length()
  }
  let mut sum = 0.0
  for index = 0; index < count; index = index + 1 {
    let delta = right[index] - left[index]
    if delta > 0.0 {
      sum += delta * delta
    }
  }
  sum.sqrt()
}

///|
/// A full feature record suitable for an event classifier.
pub fn advanced_spectral_features(
  spectrum : Spectrum,
  band : SpectralBand,
) -> AdvancedSpectralFeatures {
  let peak = advanced_spectral_peak(spectrum, band)
  {
    centroid: advanced_spectral_centroid(spectrum, band),
    bandwidth: advanced_spectral_bandwidth(spectrum, band),
    rolloff: advanced_spectral_rolloff(spectrum, band, 0.85),
    flatness: advanced_spectral_flatness(spectrum, band),
    entropy: advanced_spectral_entropy(spectrum, band),
    skewness: advanced_spectral_skewness(spectrum, band),
    kurtosis: advanced_spectral_kurtosis(spectrum, band),
    peak_frequency: peak.0,
    peak_power: peak.1,
    slope: advanced_spectral_slope(spectrum, band),
  }
}

///|
/// Format spectral features for stable logs and CSV-like reports.
pub fn advanced_spectral_features_text(
  features : AdvancedSpectralFeatures,
) -> String {
  "centroid=" +
  features.centroid.to_string() +
  " bandwidth=" +
  features.bandwidth.to_string() +
  " rolloff=" +
  features.rolloff.to_string() +
  " flatness=" +
  features.flatness.to_string() +
  " entropy=" +
  features.entropy.to_string() +
  " skewness=" +
  features.skewness.to_string() +
  " kurtosis=" +
  features.kurtosis.to_string() +
  " peak_hz=" +
  features.peak_frequency.to_string() +
  " peak_power=" +
  features.peak_power.to_string() +
  " slope=" +
  features.slope.to_string()
}