sidereon_core/

signal.rs

//! GPS L1 C/A code generation, coherent correlation, and acquisition.
//!
//! The functions here are deterministic time-domain signal primitives. They
//! intentionally keep the same simple operation order as the original Sidereon
//! implementation: C/A chips are generated by clocking the two IS-GPS-200 LFSRs,
//! sampled by zero-order hold, and correlation/acquisition are direct nested
//! sums over samples and Doppler bins.

pub mod analysis;

use crate::constants::F_L1_HZ;
use crate::tolerances::DOPPLER_GRID_EDGE_EPS_HZ;
use crate::validate;

/// GPS C/A code length in chips.
pub const CA_CODE_LENGTH: usize = 1023;

/// GPS C/A chipping rate in chips per second.
pub const CA_CHIP_RATE_HZ: f64 = 1_023_000.0;

const TWO_PI: f64 = 2.0 * std::f64::consts::PI;
const DEFAULT_DOPPLER_MIN_HZ: f64 = -2500.0;
const DEFAULT_DOPPLER_MAX_HZ: f64 = 2500.0;
const DEFAULT_DOPPLER_STEP_HZ: f64 = 500.0;
const DEFAULT_SAMPLE_RATE_HZ: f64 = 2.046e6;
const MAX_DOPPLER_BINS: usize = 4096;

/// Error returned by GPS C/A signal helpers.
#[derive(Debug, Clone, PartialEq)]
pub enum SignalError {
    /// Supported GPS space-vehicle PRNs are 1 through 32.
    UnsupportedPrn(i64),
    /// A boundary input was malformed before signal processing could start.
    InvalidInput {
        /// Name of the malformed field.
        field: &'static str,
        /// Stable validation reason.
        reason: &'static str,
    },
    /// A correlation/acquisition record had no samples.
    EmptySamples,
    /// The acquisition record was shorter than one sampled C/A-code period.
    TooShort,
}

impl core::fmt::Display for SignalError {
    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
        match self {
            Self::UnsupportedPrn(prn) => write!(f, "unsupported GPS C/A PRN {prn}"),
            Self::InvalidInput { field, reason } => {
                write!(f, "invalid signal input {field}: {reason}")
            }
            Self::EmptySamples => write!(f, "empty sample vector"),
            Self::TooShort => write!(f, "sample vector shorter than one C/A code period"),
        }
    }
}

impl std::error::Error for SignalError {}

/// One complex baseband sample.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct IqSample {
    /// In-phase component.
    pub i: f64,
    /// Quadrature component.
    pub q: f64,
}

impl IqSample {
    /// Construct a complex sample.
    pub const fn new(i: f64, q: f64) -> Self {
        Self { i, q }
    }

    /// Construct a real-valued sample with zero quadrature.
    pub const fn real(i: f64) -> Self {
        Self { i, q: 0.0 }
    }
}

/// Options for [`replica`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct ReplicaOptions {
    /// Sampling rate in hertz.
    pub sample_rate_hz: f64,
    /// Number of output samples.
    pub num_samples: usize,
    /// Initial C/A code phase in chips.
    pub code_phase_chips: f64,
    /// Code-rate Doppler in hertz.
    pub code_doppler_hz: f64,
}

impl ReplicaOptions {
    /// One C/A-code period at 2.046 MHz (two samples per chip).
    pub fn one_code_period() -> Self {
        let sample_rate_hz = DEFAULT_SAMPLE_RATE_HZ;
        let integration_time_s = CA_CODE_LENGTH as f64 / CA_CHIP_RATE_HZ;
        Self {
            sample_rate_hz,
            num_samples: (sample_rate_hz * integration_time_s).round() as usize,
            code_phase_chips: 0.0,
            code_doppler_hz: 0.0,
        }
    }
}

/// Options for [`correlate`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct CorrelateOptions {
    /// Sampling rate in hertz.
    pub sample_rate_hz: f64,
    /// Residual carrier Doppler to wipe off in hertz.
    pub doppler_hz: f64,
    /// Replica C/A code phase in chips.
    pub code_phase_chips: f64,
    /// Replica code-rate Doppler in hertz.
    pub code_doppler_hz: f64,
}

impl Default for CorrelateOptions {
    fn default() -> Self {
        Self {
            sample_rate_hz: DEFAULT_SAMPLE_RATE_HZ,
            doppler_hz: 0.0,
            code_phase_chips: 0.0,
            code_doppler_hz: 0.0,
        }
    }
}

/// Coherent correlation result.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct CorrelationResult {
    /// Real in-phase coherent sum.
    pub i: f64,
    /// Imaginary quadrature coherent sum.
    pub q: f64,
    /// Squared magnitude `i*i + q*q`.
    pub power: f64,
}

/// Options for [`acquire`].
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct AcquisitionOptions {
    /// Sampling rate in hertz.
    pub sample_rate_hz: f64,
    /// Minimum Doppler bin in hertz.
    pub doppler_min_hz: f64,
    /// Maximum Doppler bin in hertz.
    pub doppler_max_hz: f64,
    /// Doppler bin step in hertz.
    pub doppler_step_hz: f64,
}

impl Default for AcquisitionOptions {
    fn default() -> Self {
        Self {
            sample_rate_hz: DEFAULT_SAMPLE_RATE_HZ,
            doppler_min_hz: DEFAULT_DOPPLER_MIN_HZ,
            doppler_max_hz: DEFAULT_DOPPLER_MAX_HZ,
            doppler_step_hz: DEFAULT_DOPPLER_STEP_HZ,
        }
    }
}

/// Acquisition-grid metadata.
#[derive(Debug, Clone, PartialEq)]
pub struct AcquisitionGrid {
    /// Doppler bins searched, in hertz.
    pub doppler_hz: Vec<f64>,
    /// Number of code-phase bins searched.
    pub code_phase_bins: usize,
    /// Doppler step in hertz.
    pub doppler_step_hz: f64,
    /// Samples per C/A chip at the configured sample rate.
    pub samples_per_chip: f64,
}

/// Result of a 2D code-phase/Doppler acquisition search.
#[derive(Debug, Clone, PartialEq)]
pub struct AcquisitionResult {
    /// Recovered code phase in chips.
    pub code_phase_chips: f64,
    /// Recovered Doppler bin in hertz.
    pub doppler_hz: f64,
    /// Peak-to-mean-off-peak metric.
    pub peak_metric: f64,
    /// Alias for [`peak_metric`](Self::peak_metric).
    pub metric: f64,
    /// Peak correlator power.
    pub peak_power: f64,
    /// Search-grid metadata.
    pub grid: AcquisitionGrid,
}

/// Return the 1023 bipolar (`+1`/`-1`) GPS C/A chips for a PRN.
pub fn ca_code(prn: i64) -> Result<Vec<i8>, SignalError> {
    let taps = phase_select(prn)?;
    let raw = raw_code(taps);
    Ok(raw.into_iter().map(|bit| 1 - 2 * bit as i8).collect())
}

/// Return one bipolar chip at a wrapping zero-based index.
pub fn ca_chip(prn: i64, index: i64) -> Result<i8, SignalError> {
    let code = ca_code(prn)?;
    let idx = index.rem_euclid(CA_CODE_LENGTH as i64) as usize;
    Ok(code[idx])
}

/// Circular autocorrelation over all lags.
pub fn autocorrelation(code: &[i8]) -> Vec<i32> {
    (0..code.len())
        .map(|lag| correlation_at_equal_len(code, code, lag as i64))
        .collect()
}

/// Circular cross-correlation over all lags.
pub fn cross_correlation(code_a: &[i8], code_b: &[i8]) -> Result<Vec<i32>, SignalError> {
    if code_a.len() != code_b.len() {
        return Err(invalid_signal_input("code_lengths", "length mismatch"));
    }
    Ok((0..code_a.len())
        .map(|lag| correlation_at_equal_len(code_a, code_b, lag as i64))
        .collect())
}

/// Single-lag circular correlation.
pub fn correlation_at(code_a: &[i8], code_b: &[i8], lag: i64) -> Result<i32, SignalError> {
    if code_a.len() != code_b.len() {
        return Err(invalid_signal_input("code_lengths", "length mismatch"));
    }
    validate_correlation_lag(code_a.len(), lag)?;
    Ok(correlation_at_equal_len(code_a, code_b, lag))
}

fn correlation_at_equal_len(code_a: &[i8], code_b: &[i8], lag: i64) -> i32 {
    let n = code_a.len() as i64;
    let mut acc = 0_i32;
    for (i, &chip_a) in code_a.iter().enumerate() {
        let j = (i as i64 + lag).rem_euclid(n) as usize;
        acc += i32::from(chip_a) * i32::from(code_b[j]);
    }
    acc
}

fn validate_correlation_lag(len: usize, lag: i64) -> Result<(), SignalError> {
    if len == 0 || lag <= 0 {
        return Ok(());
    }
    let max_index =
        i64::try_from(len - 1).map_err(|_| invalid_signal_input("code_lengths", "out of range"))?;
    if max_index > i64::MAX - lag {
        return Err(invalid_signal_input("lag", "out of range"));
    }
    Ok(())
}

/// Build a sampled C/A-code replica.
pub fn replica(prn: i64, options: ReplicaOptions) -> Result<Vec<i8>, SignalError> {
    let sample_rate_hz = signal_positive_step(options.sample_rate_hz, "sample_rate_hz")?;
    let code_phase_chips = signal_finite(options.code_phase_chips, "code_phase_chips")?;
    let code_doppler_hz = signal_finite(options.code_doppler_hz, "code_doppler_hz")?;
    let code = ca_code(prn)?;
    Ok(sample_code(
        &code,
        options.num_samples,
        sample_rate_hz,
        code_phase_chips,
        code_doppler_hz,
    ))
}

/// Coherently correlate a sample record against a PRN replica.
pub fn correlate(
    iq: &[IqSample],
    prn: i64,
    options: CorrelateOptions,
) -> Result<CorrelationResult, SignalError> {
    if iq.is_empty() {
        return Err(SignalError::EmptySamples);
    }
    validate_iq_samples(iq, "samples")?;
    let sample_rate_hz = signal_positive_step(options.sample_rate_hz, "sample_rate_hz")?;
    let doppler_hz = signal_finite(options.doppler_hz, "doppler_hz")?;
    let code_phase_chips = signal_finite(options.code_phase_chips, "code_phase_chips")?;
    let code_doppler_hz = signal_finite(options.code_doppler_hz, "code_doppler_hz")?;
    let code = ca_code(prn)?;
    let sampled = sample_code(
        &code,
        iq.len(),
        sample_rate_hz,
        code_phase_chips,
        code_doppler_hz,
    );
    let (i, q) = correlate_against(iq, &sampled, sample_rate_hz, doppler_hz)?;
    let power = signal_finite(i * i + q * q, "correlation_power")?;
    Ok(CorrelationResult { i, q, power })
}

/// Coherent correlation against an explicit sampled code.
///
/// The sum follows `zip(iq, code)` semantics and therefore uses the shorter of
/// the two non-empty input lengths.
pub fn correlate_against(
    iq: &[IqSample],
    code: &[i8],
    fs: f64,
    doppler_hz: f64,
) -> Result<(f64, f64), SignalError> {
    if iq.is_empty() {
        return Err(SignalError::EmptySamples);
    }
    validate_iq_samples(iq, "samples")?;
    if code.is_empty() {
        return Err(invalid_signal_input("code", "empty"));
    }
    let fs = signal_positive_step(fs, "sample_rate_hz")?;
    let doppler_hz = signal_finite(doppler_hz, "doppler_hz")?;
    let w = TWO_PI * doppler_hz / fs;
    let mut acc_i = 0.0;
    let mut acc_q = 0.0;
    for (n, (sample, &c)) in iq.iter().zip(code.iter()).enumerate() {
        let theta = w * n as f64;
        let cos = theta.cos();
        let sin = theta.sin();
        let cc = c as f64;
        let di = (sample.i * cos + sample.q * sin) * cc;
        let dq = (sample.q * cos - sample.i * sin) * cc;
        acc_i += di;
        acc_q += dq;
    }
    Ok((
        signal_finite(acc_i, "correlation_i")?,
        signal_finite(acc_q, "correlation_q")?,
    ))
}

/// Acquire a PRN by direct 2D code-phase/Doppler search.
pub fn acquire(
    samples: &[IqSample],
    prn: i64,
    options: AcquisitionOptions,
) -> Result<AcquisitionResult, SignalError> {
    if samples.is_empty() {
        return Err(SignalError::EmptySamples);
    }
    validate_iq_samples(samples, "samples")?;
    let sample_rate_hz = signal_positive_step(options.sample_rate_hz, "sample_rate_hz")?;
    let doppler_min_hz = signal_finite(options.doppler_min_hz, "doppler_min_hz")?;
    let doppler_max_hz = signal_finite(options.doppler_max_hz, "doppler_max_hz")?;
    let doppler_step_hz = signal_positive_step(options.doppler_step_hz, "doppler_step_hz")?;
    signal_range_order(doppler_min_hz, doppler_max_hz, "doppler_max_hz")?;
    let options = AcquisitionOptions {
        sample_rate_hz,
        doppler_min_hz,
        doppler_max_hz,
        doppler_step_hz,
    };

    let samples_per_chip = options.sample_rate_hz / CA_CHIP_RATE_HZ;
    let samples_per_code = (samples_per_chip * CA_CODE_LENGTH as f64).round() as usize;
    if samples_per_code == 0 {
        return Err(SignalError::InvalidInput {
            field: "sample_rate_hz",
            reason: "out of range",
        });
    }
    if samples.len() < samples_per_code {
        return Err(SignalError::TooShort);
    }

    let code = ca_code(prn)?;
    do_acquire(samples, &code, options, samples_per_chip, samples_per_code)
}

/// Coherent integration loss from residual frequency error.
pub fn coherent_loss(freq_error_hz: f64, integration_time_s: f64) -> Result<f64, SignalError> {
    let freq_error_hz = signal_finite(freq_error_hz, "freq_error_hz")?;
    let integration_time_s = signal_positive_step(integration_time_s, "integration_time_s")?;
    let x = signal_finite(
        std::f64::consts::PI * freq_error_hz * integration_time_s,
        "coherent_loss",
    )?;
    if x == 0.0 {
        Ok(1.0)
    } else {
        let s = x.sin() / x;
        signal_finite(s * s, "coherent_loss")
    }
}

/// Coherent integration loss in decibels.
pub fn coherent_loss_db(freq_error_hz: f64, integration_time_s: f64) -> Result<f64, SignalError> {
    let loss = coherent_loss(freq_error_hz, integration_time_s)?;
    if loss <= 0.0 {
        return Err(invalid_signal_input("coherent_loss_db", "out of range"));
    }
    let loss_db = 10.0 * loss.log10();
    if loss_db.is_finite() {
        Ok(loss_db)
    } else {
        Err(invalid_signal_input("coherent_loss_db", "out of range"))
    }
}

/// Post-correlation predetection SNR in dB.
pub fn snr_post_db(cn0_dbhz: f64, integration_time_s: f64) -> Result<f64, SignalError> {
    let cn0_dbhz = signal_finite(cn0_dbhz, "cn0_dbhz")?;
    let integration_time_s = signal_positive_step(integration_time_s, "integration_time_s")?;
    signal_finite(cn0_dbhz + 10.0 * integration_time_s.log10(), "snr_post_db")
}

fn do_acquire(
    samples: &[IqSample],
    code: &[i8],
    options: AcquisitionOptions,
    samples_per_chip: f64,
    samples_per_code: usize,
) -> Result<AcquisitionResult, SignalError> {
    let doppler_bins = doppler_grid(
        options.doppler_min_hz,
        options.doppler_max_hz,
        options.doppler_step_hz,
    )?;

    let record = &samples[..samples_per_code];
    let base_code = sample_code(code, samples_per_code, options.sample_rate_hz, 0.0, 0.0);

    let mut grid = Vec::with_capacity(doppler_bins.len());
    for &d in &doppler_bins {
        let wiped = carrier_wipeoff(record, options.sample_rate_hz, d);
        validate_iq_samples(&wiped, "wiped samples")?;
        let powers = code_phase_powers(&wiped, &base_code);
        validate::finite_slice(&powers, "code phase powers").map_err(map_signal_input)?;
        grid.push((d, powers));
    }

    let mut peak_power = -1.0;
    let mut peak_doppler = 0.0;
    let mut peak_offset = 0_usize;
    for (d, powers) in &grid {
        for (off, &p) in powers.iter().enumerate() {
            if p > peak_power {
                peak_power = p;
                peak_doppler = *d;
                peak_offset = off;
            }
        }
    }

    let metric = peak_to_mean_off_peak(&grid, peak_power, peak_doppler, peak_offset);
    let code_phase_chips = peak_offset as f64 / samples_per_chip;

    Ok(AcquisitionResult {
        code_phase_chips,
        doppler_hz: peak_doppler,
        peak_metric: metric,
        metric,
        peak_power,
        grid: AcquisitionGrid {
            doppler_hz: doppler_bins,
            code_phase_bins: samples_per_code,
            doppler_step_hz: options.doppler_step_hz,
            samples_per_chip,
        },
    })
}

fn phase_select(prn: i64) -> Result<(usize, usize), SignalError> {
    match prn {
        1 => Ok((2, 6)),
        2 => Ok((3, 7)),
        3 => Ok((4, 8)),
        4 => Ok((5, 9)),
        5 => Ok((1, 9)),
        6 => Ok((2, 10)),
        7 => Ok((1, 8)),
        8 => Ok((2, 9)),
        9 => Ok((3, 10)),
        10 => Ok((2, 3)),
        11 => Ok((3, 4)),
        12 => Ok((5, 6)),
        13 => Ok((6, 7)),
        14 => Ok((7, 8)),
        15 => Ok((8, 9)),
        16 => Ok((9, 10)),
        17 => Ok((1, 4)),
        18 => Ok((2, 5)),
        19 => Ok((3, 6)),
        20 => Ok((4, 7)),
        21 => Ok((5, 8)),
        22 => Ok((6, 9)),
        23 => Ok((1, 3)),
        24 => Ok((4, 6)),
        25 => Ok((5, 7)),
        26 => Ok((6, 8)),
        27 => Ok((7, 9)),
        28 => Ok((8, 10)),
        29 => Ok((1, 6)),
        30 => Ok((2, 7)),
        31 => Ok((3, 8)),
        32 => Ok((4, 9)),
        _ => Err(SignalError::UnsupportedPrn(prn)),
    }
}

fn raw_code((tap_a, tap_b): (usize, usize)) -> Vec<u8> {
    let mut g1 = [1_u8; 10];
    let mut g2 = [1_u8; 10];
    let mut chips = Vec::with_capacity(CA_CODE_LENGTH);

    for _ in 0..CA_CODE_LENGTH {
        let g1_out = g1[9];
        let g2i = g2[tap_a - 1] ^ g2[tap_b - 1];
        chips.push(g1_out ^ g2i);
        step_g1(&mut g1);
        step_g2(&mut g2);
    }

    chips
}

fn step_g1(g1: &mut [u8; 10]) {
    let feedback = g1[2] ^ g1[9];
    shift(g1, feedback);
}

fn step_g2(g2: &mut [u8; 10]) {
    let feedback = g2[1] ^ g2[2] ^ g2[5] ^ g2[7] ^ g2[8] ^ g2[9];
    shift(g2, feedback);
}

fn shift(reg: &mut [u8; 10], feedback: u8) {
    for i in (1..reg.len()).rev() {
        reg[i] = reg[i - 1];
    }
    reg[0] = feedback;
}

fn sample_code(code: &[i8], n: usize, fs: f64, code_phase: f64, code_doppler: f64) -> Vec<i8> {
    let code_rate = CA_CHIP_RATE_HZ * (1.0 + code_doppler / F_L1_HZ);
    let per_sample = code_rate / fs;
    let len = code.len() as i64;

    (0..n)
        .map(|k| {
            let pos = code_phase + k as f64 * per_sample;
            let idx = (pos.floor() as i64).rem_euclid(len) as usize;
            code[idx]
        })
        .collect()
}

fn carrier_wipeoff(iq: &[IqSample], fs: f64, doppler_hz: f64) -> Vec<IqSample> {
    let w = TWO_PI * doppler_hz / fs;
    iq.iter()
        .enumerate()
        .map(|(k, sample)| {
            let theta = w * k as f64;
            let cos = theta.cos();
            let sin = theta.sin();
            IqSample {
                i: sample.i * cos + sample.q * sin,
                q: sample.q * cos - sample.i * sin,
            }
        })
        .collect()
}

fn code_phase_powers(wiped: &[IqSample], base_code: &[i8]) -> Vec<f64> {
    let n = wiped.len();
    (0..n)
        .map(|offset| {
            let mut i = 0.0;
            let mut q = 0.0;
            for k in 0..n {
                let sample = wiped[k];
                let c = base_code[(k + offset) % n] as f64;
                i += sample.i * c;
                q += sample.q * c;
            }
            i * i + q * q
        })
        .collect()
}

fn peak_to_mean_off_peak(
    grid: &[(f64, Vec<f64>)],
    peak_power: f64,
    peak_doppler: f64,
    peak_offset: usize,
) -> f64 {
    let n = grid.first().map_or(0, |(_, powers)| powers.len());
    let mut sum = 0.0;
    let mut count = 0_usize;

    for (d, powers) in grid {
        for (off, power) in powers.iter().enumerate() {
            if *d == peak_doppler && abs_circular_diff(off, peak_offset, n) <= 1 {
                continue;
            }
            sum += *power;
            count += 1;
        }
    }

    if count == 0 {
        0.0
    } else if sum <= 0.0 && peak_power > 0.0 {
        1.0e12
    } else if sum <= 0.0 {
        0.0
    } else {
        peak_power / (sum / count as f64)
    }
}

fn doppler_grid(dmin: f64, dmax: f64, dstep: f64) -> Result<Vec<f64>, SignalError> {
    let dstep = signal_positive_step(dstep, "doppler_step_hz")?;
    let last_bin_index = doppler_last_bin_index(dmin, dmax, dstep)?;
    Ok((0..=last_bin_index)
        .map(|k| dmin + k as f64 * dstep)
        .filter(|d| *d <= dmax + DOPPLER_GRID_EDGE_EPS_HZ)
        .collect())
}

fn doppler_last_bin_index(dmin: f64, dmax: f64, dstep: f64) -> Result<usize, SignalError> {
    let last_bin_index = ((dmax - dmin) / dstep).round();
    if !last_bin_index.is_finite() || last_bin_index < 0.0 {
        return Err(invalid_signal_input("doppler_grid", "out of range"));
    }
    let bin_count = last_bin_index + 1.0;
    if bin_count > MAX_DOPPLER_BINS as f64 {
        return Err(invalid_signal_input("doppler_grid", "out of range"));
    }
    Ok(last_bin_index as usize)
}

fn signal_positive_step(x: f64, field: &'static str) -> Result<f64, SignalError> {
    validate::positive_step(x, field).map_err(map_signal_input)
}

fn signal_finite(x: f64, field: &'static str) -> Result<f64, SignalError> {
    validate::finite(x, field).map_err(map_signal_input)
}

fn signal_range_order(lo: f64, hi: f64, field: &'static str) -> Result<(), SignalError> {
    validate::range_order(lo, hi, field).map_err(map_signal_input)
}

fn validate_iq_samples(samples: &[IqSample], field: &'static str) -> Result<(), SignalError> {
    for sample in samples {
        if !sample.i.is_finite() || !sample.q.is_finite() {
            return Err(invalid_signal_input(field, "not finite"));
        }
    }
    Ok(())
}

fn map_signal_input(error: validate::FieldError) -> SignalError {
    invalid_signal_input(error.field(), error.reason())
}

fn invalid_signal_input(field: &'static str, reason: &'static str) -> SignalError {
    SignalError::InvalidInput { field, reason }
}

fn abs_circular_diff(a: usize, b: usize, n: usize) -> usize {
    let d = a.abs_diff(b) % n;
    d.min(n - d)
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn unsupported_prn_is_tagged() {
        assert_eq!(ca_code(33), Err(SignalError::UnsupportedPrn(33)));
        assert_eq!(ca_chip(0, 0), Err(SignalError::UnsupportedPrn(0)));
    }

    #[test]
    fn code_balance_and_correlation_shape_are_pinned() {
        for prn in 1..=32 {
            let code = ca_code(prn).unwrap();
            assert_eq!(code.len(), CA_CODE_LENGTH);
            assert_eq!(code.iter().filter(|&&chip| chip == -1).count(), 512);
            assert_eq!(code.iter().filter(|&&chip| chip == 1).count(), 511);
            assert_eq!(code.iter().map(|&chip| i32::from(chip)).sum::<i32>(), -1);
        }

        let code = ca_code(1).unwrap();
        let corr = autocorrelation(&code);
        assert_eq!(corr[0], 1023);
        assert!(!corr[1..].contains(&1023));
        let mut values = corr[1..].to_vec();
        values.sort_unstable();
        values.dedup();
        assert_eq!(values, vec![-65, -1, 63]);
    }

    #[test]
    fn loss_and_snr_primitives_are_deterministic() {
        assert_eq!(
            coherent_loss(0.0, 1.0e-3).unwrap().to_bits(),
            1.0_f64.to_bits()
        );
        assert_eq!(
            snr_post_db(40.0, 1.0e-3).unwrap().to_bits(),
            10.0_f64.to_bits()
        );
        assert_eq!(
            coherent_loss(f64::MAX, 1.0),
            Err(invalid_signal_input("coherent_loss", "not finite"))
        );
    }

    #[test]
    fn correlation_rejects_nonfinite_derived_outputs() {
        let samples = [IqSample::real(f64::MAX), IqSample::real(f64::MAX)];
        let code = [1_i8, 1_i8];

        assert_eq!(
            correlate_against(&samples, &code, DEFAULT_SAMPLE_RATE_HZ, 0.0),
            Err(invalid_signal_input("correlation_i", "not finite"))
        );
        assert_eq!(
            correlate(
                &samples[..1],
                1,
                CorrelateOptions {
                    sample_rate_hz: DEFAULT_SAMPLE_RATE_HZ,
                    doppler_hz: 0.0,
                    code_phase_chips: 0.0,
                    code_doppler_hz: 0.0,
                }
            ),
            Err(invalid_signal_input("correlation_power", "not finite"))
        );
    }
}