micarrsim/packages/acoustics-core/sampling.js

export const MAX_PRECISE_ANALOG_FREQUENCY_HZ = 1_000_000;
export const MAX_EXPLICIT_FOLD_SEGMENTS = 513;

function requireSampleRate(sampleRate) {
  if (!Number.isFinite(sampleRate) || sampleRate <= 0) {
    throw new RangeError("sampleRate must be a positive finite number.");
  }
  return sampleRate;
}

function requireFrequency(frequencyHz, label = "frequencyHz") {
  if (
    !Number.isFinite(frequencyHz) ||
    frequencyHz < 0 ||
    frequencyHz > MAX_PRECISE_ANALOG_FREQUENCY_HZ
  ) {
    throw new RangeError(
      `${label} must be between 0 and ${MAX_PRECISE_ANALOG_FREQUENCY_HZ} Hz.`,
    );
  }
  return frequencyHz;
}

function positiveModulo(value, modulus) {
  const remainder = value % modulus;
  return remainder < 0 ? remainder + modulus : remainder;
}

function normalizedZero(value) {
  return Object.is(value, -0) ? 0 : value;
}

/**
 * Describe the frequency that an ideal real-valued sampler records.
 *
 * No anti-alias filter is implied. A negative signed ADC frequency means that
 * the positive-frequency spectrum is mirrored when represented on 0..Nyquist.
 */
export function describeFrequencyAliasing(frequencyHz, sampleRate) {
  requireFrequency(frequencyHz);
  requireSampleRate(sampleRate);
  const nyquistHz = sampleRate / 2;
  const wrappedHz = positiveModulo(frequencyHz, sampleRate);
  const signedAdcFrequencyHz =
    wrappedHz <= nyquistHz ? wrappedHz : wrappedHz - sampleRate;
  const adcFrequencyHz = Math.abs(signedAdcFrequencyHz);
  const dcEdge = adcFrequencyHz === 0;
  const nyquistEdge = adcFrequencyHz === nyquistHz;
  return {
    model: "ideal-real-sampler-no-antialias-filter",
    sampleRate,
    nyquistHz,
    analogFrequencyHz: frequencyHz,
    adcFrequencyHz: normalizedZero(adcFrequencyHz),
    signedAdcFrequencyHz: normalizedZero(signedAdcFrequencyHz),
    aliased: frequencyHz > nyquistHz,
    spectralInversion: signedAdcFrequencyHz < 0,
    dcEdge,
    nyquistEdge,
    phaseDegenerate: dcEdge || nyquistEdge,
  };
}

export function foldFrequencyToNyquist(frequencyHz, sampleRate) {
  return describeFrequencyAliasing(frequencyHz, sampleRate).adcFrequencyHz;
}

function containsResidue(first, last, residue, modulus) {
  const firstResidue = positiveModulo(first, modulus);
  const offset = positiveModulo(residue - firstResidue, modulus);
  return first + offset <= last;
}

function contiguousBinRanges(bins, binWidthHz) {
  const ranges = [];
  let firstBin = bins[0];
  let lastBin = firstBin;
  for (let index = 1; index <= bins.length; index += 1) {
    const bin = bins[index];
    if (bin === lastBin + 1) {
      lastBin = bin;
      continue;
    }
    ranges.push({
      firstBin,
      lastBin,
      firstBinCenterHz: firstBin * binWidthHz,
      lastBinCenterHz: lastBin * binWidthHz,
    });
    firstBin = bin;
    lastBin = bin;
  }
  return ranges;
}

function foldBandSegments(minHz, maxHz, sampleRate) {
  const nyquistHz = sampleRate / 2;
  const firstZoneIndex = Math.floor(minHz / nyquistHz);
  const lastZoneIndex = Math.ceil(maxHz / nyquistHz) - 1;
  const segmentCount = lastZoneIndex - firstZoneIndex + 1;
  if (
    !Number.isSafeInteger(firstZoneIndex) ||
    !Number.isSafeInteger(lastZoneIndex) ||
    !Number.isSafeInteger(segmentCount) ||
    segmentCount < 1 ||
    segmentCount > MAX_EXPLICIT_FOLD_SEGMENTS
  ) {
    throw new RangeError(
      `band crosses ${segmentCount} Nyquist zones; the explicit limit is ` +
        `${MAX_EXPLICIT_FOLD_SEGMENTS}.`,
    );
  }
  const segments = [];
  let analogStartHz = minHz;
  while (analogStartHz < maxHz) {
    const zoneIndex = Math.floor(analogStartHz / nyquistHz);
    const nyquistZone = zoneIndex + 1;
    const nextBoundaryHz = nyquistZone * nyquistHz;
    const analogEndHz = Math.min(maxHz, nextBoundaryHz);
    if (!(analogEndHz > analogStartHz)) {
      throw new RangeError("band is too large to retain Nyquist-zone precision.");
    }
    const inverted = nyquistZone % 2 === 0;
    const signedAdcStartHz = inverted
      ? analogStartHz - nyquistZone * nyquistHz
      : analogStartHz - (nyquistZone - 1) * nyquistHz;
    const signedAdcEndHz = inverted
      ? analogEndHz - nyquistZone * nyquistHz
      : analogEndHz - (nyquistZone - 1) * nyquistHz;
    segments.push({
      analogMinHz: analogStartHz,
      analogMaxHz: analogEndHz,
      signedAdcStartHz: normalizedZero(signedAdcStartHz),
      signedAdcEndHz: normalizedZero(signedAdcEndHz),
      inverted,
      aliasOrder: Math.floor(nyquistZone / 2),
      nyquistZone,
    });
    analogStartHz = analogEndHz;
  }
  return segments;
}

function hasContinuousFoldOverlap(segments) {
  const adcRanges = segments.map((segment) => ({
    minHz: Math.min(
      Math.abs(segment.signedAdcStartHz),
      Math.abs(segment.signedAdcEndHz),
    ),
    maxHz: Math.max(
      Math.abs(segment.signedAdcStartHz),
      Math.abs(segment.signedAdcEndHz),
    ),
  })).sort((left, right) => left.minHz - right.minHz || left.maxHz - right.maxHz);
  let maximumEndHz = -Infinity;
  for (const range of adcRanges) {
    if (range.minHz < maximumEndHz) return true;
    maximumEndHz = Math.max(maximumEndHz, range.maxHz);
  }
  return false;
}

function adcBinSetKey(sampleRate, fftSize, ranges) {
  const rangeKey = ranges
    .map(({ firstBin, lastBin }) =>
      firstBin === lastBin ? String(firstBin) : `${firstBin}-${lastBin}`)
    .join(",");
  return `fs:${sampleRate}|fft:${fftSize}|bins:${rangeKey}`;
}

/**
 * Map a requested pre-sampling frequency interval to the observable one-sided
 * real FFT bins. Every matching bin is retained; this function does not choose
 * an alias order or steer with the original carrier.
 */
export function observeFrequencyBand(band, sampleRate, fftSize) {
  requireSampleRate(sampleRate);
  if (!Number.isInteger(fftSize) || fftSize < 2 || fftSize % 2 !== 0) {
    throw new RangeError("fftSize must be a positive even integer.");
  }
  if (!band || typeof band !== "object") {
    throw new TypeError("band must be an object.");
  }
  requireFrequency(band.minHz, "band.minHz");
  requireFrequency(band.maxHz, "band.maxHz");
  if (band.minHz >= band.maxHz) {
    throw new RangeError("band.minHz must be lower than band.maxHz.");
  }

  const binWidthHz = sampleRate / fftSize;
  const firstAnalogBin = Math.ceil(band.minHz / binWidthHz);
  const lastAnalogBin = Math.floor(band.maxHz / binWidthHz);
  if (
    !Number.isSafeInteger(firstAnalogBin) ||
    !Number.isSafeInteger(lastAnalogBin)
  ) {
    throw new RangeError("band frequency is too large for safe FFT-bin mapping.");
  }
  if (firstAnalogBin > lastAnalogBin) {
    throw new Error(
      `${band.minHz}–${band.maxHz} Hz contains no FFT-grid frequency at ${binWidthHz} Hz spacing.`,
    );
  }

  const bins = [];
  const halfSize = fftSize / 2;
  for (let adcBin = 0; adcBin <= halfSize; adcBin += 1) {
    const positiveResidue = adcBin;
    const negativeResidue = adcBin === 0 || adcBin === halfSize
      ? adcBin
      : fftSize - adcBin;
    if (
      containsResidue(firstAnalogBin, lastAnalogBin, positiveResidue, fftSize) ||
      containsResidue(firstAnalogBin, lastAnalogBin, negativeResidue, fftSize)
    ) {
      bins.push(adcBin);
    }
  }
  if (bins.length === 0) {
    throw new Error(`${band.minHz}–${band.maxHz} Hz maps to no observable FFT bin.`);
  }

  const nyquistHz = sampleRate / 2;
  const adcBinCenterRanges = contiguousBinRanges(bins, binWidthHz);
  const foldSegments = foldBandSegments(band.minHz, band.maxHz, sampleRate);
  const analogBinCount = lastAnalogBin - firstAnalogBin + 1;
  const continuousManyToOne = hasContinuousFoldOverlap(foldSegments);
  const gridManyToOne = analogBinCount > bins.length;
  return {
    model: "ideal-real-sampler-no-antialias-filter",
    sampleRate,
    nyquistHz,
    fftSize,
    binWidthHz,
    requestedMinHz: band.minHz,
    requestedMaxHz: band.maxHz,
    aliased: band.maxHz > nyquistHz,
    crossesFirstNyquist: band.minHz < nyquistHz && band.maxHz > nyquistHz,
    localizationInputModel: "observed-one-sided-real-fft-bins",
    selectionModel: "fft-bin-center-membership",
    originalCarrierUsed: false,
    foldBoundaryCount: Math.max(0, foldSegments.length - 1),
    continuousManyToOne,
    gridManyToOne,
    foldSegments,
    bins,
    adcBinCenterRanges,
    adcBinSetKey: adcBinSetKey(sampleRate, fftSize, adcBinCenterRanges),
  };
}