///|
/// One-sided periodogram bin.
pub(all) struct SpectrumBin {
frequency_hz : Double
power : Double
amplitude : Double
} derive(FromJson, ToJson, Debug, Eq)
///|
/// Power accumulated in a named frequency band.
pub(all) struct FrequencyBandPower {
name : String
lower_hz : Double
upper_hz : Double
power : Double
normalized_power : Double
} derive(FromJson, ToJson, Debug, Eq)
///|
/// Frequency-domain HRV descriptors.
pub(all) struct FrequencyMetrics {
sample_rate_hz : Double
total_power : Double
vlf : FrequencyBandPower
lf : FrequencyBandPower
hf : FrequencyBandPower
lf_hf_ratio : Double
spectral_centroid_hz : Double
spectral_entropy : Double
peak_frequency_hz : Double
spectrum : Array[SpectrumBin]
} derive(FromJson, ToJson, Debug, Eq)
///|
/// Calculate a real-valued DFT periodogram for a sampled sequence.
pub fn calculate_periodogram(
values : Array[Double],
sample_rate_hz : Double,
) -> Array[SpectrumBin] {
let result = []
let n = values.length()
if n == 0 || sample_rate_hz <= 0.0 {
return result
}
let max_k = n / 2
for k in 0..<=max_k {
let mut real = 0.0
let mut imaginary = 0.0
for i in 0.. Double {
let mut total = 0.0
for bin in spectrum {
if bin.frequency_hz >= lower_hz && bin.frequency_hz < upper_hz {
total += bin.power
}
}
total
}
///|
/// Return the total one-sided power.
pub fn total_spectral_power(spectrum : Array[SpectrumBin]) -> Double {
let mut total = 0.0
for bin in spectrum {
total += bin.power
}
total
}
///|
/// Return the largest non-DC spectrum bin.
pub fn dominant_spectrum_bin(spectrum : Array[SpectrumBin]) -> SpectrumBin {
if spectrum.length() == 0 {
return { frequency_hz: 0.0, power: 0.0, amplitude: 0.0 }
}
let mut best = spectrum[0]
for bin in spectrum {
if bin.power > best.power && bin.frequency_hz > 0.0 {
best = bin
}
}
best
}
///|
/// Calculate normalized spectral entropy from non-negative powers.
pub fn spectral_entropy(spectrum : Array[SpectrumBin]) -> Double {
let total = total_spectral_power(spectrum)
if total <= 0.0 {
return 0.0
}
let mut entropy = 0.0
let mut count = 0
for bin in spectrum {
if bin.power > 0.0 {
let probability = bin.power / total
entropy -= probability * @math.ln(probability)
count += 1
}
}
if count <= 1 {
0.0
} else {
entropy / @math.ln(count.to_double())
}
}
///|
/// Calculate the power-weighted average frequency.
pub fn spectral_centroid(spectrum : Array[SpectrumBin]) -> Double {
let total = total_spectral_power(spectrum)
if total <= 0.0 {
0.0
} else {
let mut weighted = 0.0
for bin in spectrum {
weighted += bin.frequency_hz * bin.power
}
weighted / total
}
}
///|
/// Return the first frequency whose cumulative power reaches a proportion.
pub fn spectral_edge_frequency(
spectrum : Array[SpectrumBin],
proportion : Double,
) -> Double {
let total = total_spectral_power(spectrum)
if total <= 0.0 {
return 0.0
}
let target = total * proportion.clamp(min=0.0, max=1.0)
let mut cumulative = 0.0
for bin in spectrum {
cumulative += bin.power
if cumulative >= target {
return bin.frequency_hz
}
}
if spectrum.length() == 0 {
0.0
} else {
spectrum[spectrum.length() - 1].frequency_hz
}
}
///|
/// Create a named frequency band from a spectrum.
pub fn make_frequency_band(
name : String,
lower_hz : Double,
upper_hz : Double,
spectrum : Array[SpectrumBin],
total_power : Double,
) -> FrequencyBandPower {
let power = integrate_band_power(spectrum, lower_hz, upper_hz)
{
name,
lower_hz,
upper_hz,
power,
normalized_power: if total_power == 0.0 {
0.0
} else {
power / total_power
},
}
}
///|
/// Compute conventional VLF, LF, and HF HRV bands.
pub fn calculate_frequency_metrics(
intervals : Array[Double],
sample_rate_hz : Double,
) -> FrequencyMetrics {
let tachogram = prepare_tachogram(intervals, sample_rate_hz, true, Hann)
let spectrum = calculate_periodogram(tachogram.values_ms, sample_rate_hz)
let total_power = total_spectral_power(spectrum)
let vlf = make_frequency_band("VLF", 0.0033, 0.04, spectrum, total_power)
let lf = make_frequency_band("LF", 0.04, 0.15, spectrum, total_power)
let hf = make_frequency_band("HF", 0.15, 0.40, spectrum, total_power)
let peak = dominant_spectrum_bin(spectrum)
{
sample_rate_hz,
total_power,
vlf,
lf,
hf,
lf_hf_ratio: if hf.power == 0.0 {
0.0
} else {
lf.power / hf.power
},
spectral_centroid_hz: spectral_centroid(spectrum),
spectral_entropy: spectral_entropy(spectrum),
peak_frequency_hz: peak.frequency_hz,
spectrum,
}
}
///|
/// Calculate the average spectrum across overlapping Welch frames.
pub fn welch_periodogram(
values : Array[Double],
sample_rate_hz : Double,
frame_size : Int,
hop_size : Int,
function : WindowFunction,
) -> Array[SpectrumBin] {
let frames = frame_values(values, frame_size, hop_size)
let result = []
if frames.length() == 0 || sample_rate_hz <= 0.0 {
return result
}
let first = calculate_periodogram(
prepare_spectral_values(frames[0], true, function),
sample_rate_hz,
)
for bin in first {
result.push({ frequency_hz: bin.frequency_hz, power: 0.0, amplitude: 0.0 })
}
for frame in frames {
let spectrum = calculate_periodogram(
prepare_spectral_values(frame, true, function),
sample_rate_hz,
)
for i in 0.. Double {
if lag < 0 || lag >= values.length() {
return 0.0
}
let average = mean_value(values)
let mut numerator = 0.0
let mut denominator = 0.0
for i in 0.. Int {
if values.length() <= 1 {
return 0
}
let low = if minimum_lag < 1 { 1 } else { minimum_lag }
let high = if maximum_lag >= values.length() {
values.length() - 1
} else {
maximum_lag
}
if low > high {
return 0
}
let mut best_lag = low
let mut best_value = autocorrelation(values, low)
for lag in (low + 1)..<=high {
let value = autocorrelation(values, lag)
if value > best_value {
best_lag = lag
best_value = value
}
}
best_lag
}
///|
/// Return a concise frequency feature vector for model inputs.
pub fn frequency_feature_vector(
intervals : Array[Double],
sample_rate_hz : Double,
) -> Array[Double] {
let metrics = calculate_frequency_metrics(intervals, sample_rate_hz)
[
metrics.total_power,
metrics.vlf.power,
metrics.lf.power,
metrics.hf.power,
metrics.vlf.normalized_power,
metrics.lf.normalized_power,
metrics.hf.normalized_power,
metrics.lf_hf_ratio,
metrics.spectral_centroid_hz,
metrics.spectral_entropy,
metrics.peak_frequency_hz,
]
}