diff --git a/.gitignore b/.gitignore index cb65245..6dd5510 100644 --- a/.gitignore +++ b/.gitignore @@ -4,3 +4,4 @@ dump978-fa skyaware978 faup978 fec_tests +adaptive_tests diff --git a/Makefile b/Makefile index 47532a6..6bbfff5 100644 --- a/Makefile +++ b/Makefile @@ -14,7 +14,7 @@ LIBS_SDR=-lSoapySDR all: dump978-fa skyaware978 -dump978-fa: dump978_main.o socket_output.o message_dispatch.o fec.o libs/fec/init_rs_char.o libs/fec/decode_rs_char.o sample_source.o soapy_source.o convert.o demodulator.o uat_message.o stratux_serial.o +dump978-fa: dump978_main.o socket_output.o message_dispatch.o fec.o libs/fec/init_rs_char.o libs/fec/decode_rs_char.o sample_source.o soapy_source.o convert.o demodulator.o uat_message.o stratux_serial.o adaptive.o $(CXX) $(CXXFLAGS) $(LDFLAGS) $^ -o $@ $(LIBS) $(LIBS_SDR) faup978: faup978_main.o socket_input.o uat_message.o track.o faup978_reporter.o @@ -23,14 +23,18 @@ faup978: faup978_main.o socket_input.o uat_message.o track.o faup978_reporter.o skyaware978: skyaware978_main.o socket_input.o uat_message.o track.o skyaware_writer.o $(CXX) $(CXXFLAGS) $(LDFLAGS) $^ -o $@ $(LIBS) -test: fec_tests +test: fec_tests adaptive_tests ./fec_tests + ./adaptive_tests fec_tests: fec_tests.o libs/fec/init_rs_char.o libs/fec/decode_rs_char.o libs/fec/encode_rs_char.o $(CXX) $(CXXFLAGS) $(LDFLAGS) $^ -o $@ -lstdc++ +adaptive_tests: adaptive_tests.o adaptive.o + $(CXX) $(CXXFLAGS) $(LDFLAGS) $^ -o $@ -lstdc++ -lm + format: clang-format -style=file -i *.cc *.h clean: - rm -f *.o libs/fec/*.o dump978-fa faup978 skyaware978 fec_tests + rm -f *.o libs/fec/*.o dump978-fa faup978 skyaware978 fec_tests adaptive_tests diff --git a/README.md b/README.md index c2ee410..68a7ac1 100644 --- a/README.md +++ b/README.md @@ -65,6 +65,39 @@ The main options are: Pass `--help` for a full list of options. +## Adaptive gain control + +When receiving samples from an SDR (`--sdr` only), dump978-fa can adjust the +SDR gain automatically. This is a port of +[dump1090-fa's adaptive gain feature](https://github.com/flightaware/dump1090/blob/master/README.adaptive-gain.md) +and works the same way. Two independent modes are available; both are off by +default and can be enabled together: + + * `--adaptive-range` enables **dynamic range mode**: dump978 measures the + noise floor and searches for the highest gain that still provides a target + dynamic range (by default, 60% of the SDR's gain range in dB - + approximately 30dB for an RTLSDR). This suits most installs. + * `--adaptive-burst` enables **burst mode**: dump978 watches for loud, + message-length bursts of signal that failed to decode - a sign that a + nearby transmitter is overloading the receiver - and reduces gain when + they arrive too often, raising it again once things quieten down. + +Useful tuning options: + + * `--adaptive-min-gain ` / `--adaptive-max-gain ` restrict the gain + range that adaptive gain control will use. If you know an approximately + correct gain for your installation, setting limits around it speeds up + convergence. + * `--adaptive-duty-cycle ` controls the fraction of received + samples inspected (default: 50). Lower values reduce CPU use on slow + machines at a small cost in measurement accuracy. + * `--adaptive-range-target ` overrides the target dynamic range. + +Further tunables mirroring dump1090-fa's options are available; pass `--help` +for the full list. Adaptive gain control requires manual gain control of the +SDR and cannot be combined with `--sdr-auto-gain`. Gain changes and a +once-a-minute status line are logged to stderr. + ## Third-party code Third-party source code included in libs/: diff --git a/adaptive.cc b/adaptive.cc new file mode 100644 index 0000000..dddc61c --- /dev/null +++ b/adaptive.cc @@ -0,0 +1,622 @@ +// adaptive.cc: adaptive gain control +// +// Ported from dump1090's adaptive.c +// (https://github.com/flightaware/dump1090) +// +// Copyright (c) 2021 FlightAware, LLC +// +// This file is free software: you may copy, redistribute and/or modify it +// under the terms of the GNU General Public License as published by the +// Free Software Foundation, either version 2 of the License, or (at your +// option) any later version. +// +// This file is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#include "adaptive.h" + +#include +#include +#include + +using namespace flightaware::uat; + +// +// block handling +// +// 1 block = approx 1 second of samples. Control updates are done at the end of each block only. +// Each block is made up of an integer number of subblocks (currently 20) +// +// 1 subblock = approx 50ms of samples. Duty cycle decisions are made at the subblock level; +// either the whole subblock is processed, or the whole subblock is skipped. +// Each subblock is made up of an integer number of windows (currently 1250) +// +// 1 window = approx 40us of samples. Burst measurements are made by counting samples within each window. +// +// All three levels are aligned, i.e. every block boundary is also a subblock boundary; +// every subblock boundary is also a window boundary. +// +// Duty cycle is expressed as N/D where D = kSubblocksPerBlock and +// N = round(kSubblocksPerBlock * duty_cycle); i.e. within each block, there are exactly +// N active subblocks out of D total subblocks. The active subblocks are distributed +// evenly across the block by increasing a counter by N on each subblock, modulo D, and +// marking the subblock as active each time the counter rolls over. + +static const unsigned kSubblocksPerBlock = 20; // subblocks per block +static const unsigned kReportInterval = 60; // blocks (~seconds) between periodic reports + +// Return the dB value of gain step 'step', clamping out-of-range steps to the table bounds +double AdaptiveGainControl::StepDb(int step) const { + if (gain_steps_.empty()) + return 0; + if (step < 0) + step = 0; + if (step >= GainStepCount()) + step = GainStepCount() - 1; + return gain_steps_[step]; +} + +// Return the gain step closest to 'db' +int AdaptiveGainControl::NearestStep(double db) const { + int best = 0; + for (int i = 1; i < GainStepCount(); ++i) { + if (std::fabs(gain_steps_[i] - db) < std::fabs(gain_steps_[best] - db)) + best = i; + } + return best; +} + +// Try to change the SDR gain to 'step' and tell the user about it, +// with 'why' as the reason to show. Return true if the gain actually changed. +bool AdaptiveGainControl::SetGainStep(int step, const char *why) { + if (step < gain_min_) + step = gain_min_; + if (step > gain_max_) + step = gain_max_; + + if (current_step_ == step) + return false; + + fprintf(stderr, "adaptive: changing gain from %.1fdB (step %d) to %.1fdB (step %d) because: %s\n", StepDb(current_step_), current_step_, StepDb(step), step, why); + + double reported = sdr_->SetGainDb(StepDb(step)); + if (std::fabs(reported - StepDb(step)) >= 0.5) + fprintf(stderr, "adaptive: device reports gain %.1fdB after requesting %.1fdB\n", reported, StepDb(step)); + + current_step_ = step; + ++report_gain_changes_; + return true; +} + +// Update internal state to reflect a gain change +// (usually after SetGainStep returns true, but also called during init) +void AdaptiveGainControl::GainChanged() { + gain_up_db_ = StepDb(current_step_ + 1) - StepDb(current_step_); + gain_down_db_ = StepDb(current_step_) - StepDb(current_step_ - 1); + + double loud_threshold_dbfs = 0 - gain_up_db_ - 3.0; + burst_loud_threshold_ = std::pow(10, loud_threshold_dbfs / 10.0); + + range_change_timer_ = params_.range_change_delay; + burst_change_timer_ = params_.burst_change_delay; + burst_loud_blocks_ = 0; + burst_quiet_blocks_ = 0; +} + +void AdaptiveGainControl::IncreaseGain(const char *why) { + if (SetGainStep(current_step_ + 1, why)) + GainChanged(); +} + +void AdaptiveGainControl::DecreaseGain(const char *why) { + if (SetGainStep(current_step_ - 1, why)) + GainChanged(); +} + +AdaptiveGainControl::AdaptiveGainControl(const AdaptiveParams ¶ms, SdrGainControl::Pointer sdr, unsigned sample_rate) : params_(params), sdr_(sdr) { + double range_min = 0, range_max = 0, range_step = 0; + sdr_->GetGainLimits(range_min, range_max, range_step); + + // If the SDR doesn't support gain control, disable ourselves + if (!(range_max > range_min)) { + if (params_.burst_control || params_.range_control) { + fprintf(stderr, "warning: adaptive gain control requested, but SDR gain control not available, ignored.\n"); + } + params_.burst_control = false; + params_.range_control = false; + } + + // If we're disabled, do nothing + if (!params_.burst_control && !params_.range_control) + return; + + // Build the gain step table. Unlike dump1090, which sees the tuner's discrete + // gain table directly, SoapySDR reports gain as a min/max/step range; if the + // device does not report a usable step size, probe in 1dB steps. + double step_db = range_step; + if (step_db <= 0 || (range_max - range_min) / step_db > 512) + step_db = 1.0; + int step_count = (int)((range_max - range_min) / step_db + 1e-6) + 1; + for (int i = 0; i < step_count; ++i) + gain_steps_.push_back(range_min + i * step_db); + + // Autoselect a target dynamic range if none was given + // (matches dump1090's default for SoapySDR devices; approximately 30dB for an RTLSDR) + if (params_.range_target <= 0) + params_.range_target = (range_max - range_min) * 0.6; + + // Set up window, subblock, and block sizes + // Look for 40us bursts + samples_per_window_ = sample_rate / 25000; + + // Use ~50ms subblocks; ensure it's an exact multiple of window size + samples_per_subblock_ = samples_per_window_ * 1250; + + subblocks_remaining_ = kSubblocksPerBlock; + subblock_samples_remaining_ = samples_per_subblock_; + subblock_active_ = false; + + double N = std::round(kSubblocksPerBlock * params_.duty_cycle); + if (N <= 0) + N = 1; + if (N > kSubblocksPerBlock) + N = kSubblocksPerBlock; + fprintf(stderr, "adaptive: using %.0f%% duty cycle\n", 100.0 * N / kSubblocksPerBlock); + dutycycle_N_ = (unsigned)N; + + burst_window_remaining_ = samples_per_window_; + burst_window_counter_ = 0; + + range_radix_.assign(65536, 0); + range_state_ = RangeState::RESCAN_UP; + + // select and enforce gain limits + int maxgain = GainStepCount() - 1; + for (gain_min_ = 0; gain_min_ < maxgain; ++gain_min_) { + if (StepDb(gain_min_) >= params_.min_gain_db) + break; + } + + for (gain_max_ = maxgain; gain_max_ > gain_min_; --gain_max_) { + if (StepDb(gain_max_) <= params_.max_gain_db) + break; + } + + current_step_ = NearestStep(sdr_->GetGainDb()); + + fprintf(stderr, "adaptive: enabled adaptive gain control with gain limits %.1fdB (step %d) .. %.1fdB (step %d)\n", StepDb(gain_min_), gain_min_, StepDb(gain_max_), gain_max_); + if (params_.range_control) + fprintf(stderr, "adaptive: enabled dynamic range control, target dynamic range %.1fdB\n", params_.range_target); + if (params_.burst_control) + fprintf(stderr, "adaptive: enabled burst control\n"); + SetGainStep(current_step_, "constraining gain to adaptive gain limits"); + GainChanged(); + + range_gain_limit_ = current_step_; + report_timer_ = kReportInterval; +} + +// Feed some samples into the adaptive system. Any number of samples might be passed in. +void AdaptiveGainControl::Update(const std::uint16_t *buf, unsigned length, const double *decoded_power) { + if (!params_.burst_control && !params_.range_control) + return; + + // process complete subblocks + while (length >= subblock_samples_remaining_) { + if (subblock_active_) + UpdateSubblock(buf, subblock_samples_remaining_, decoded_power); + + buf += subblock_samples_remaining_; + length -= subblock_samples_remaining_; + subblock_samples_remaining_ = samples_per_subblock_; + + dutycycle_counter_ += dutycycle_N_; + if (dutycycle_counter_ >= kSubblocksPerBlock) { + dutycycle_counter_ -= kSubblocksPerBlock; + subblock_active_ = true; + } else { + subblock_active_ = false; + // fake a quiet window to reset any existing run + BurstEndOfWindow(0); + } + + if (!--subblocks_remaining_) { + // Block completed, do a control update + subblocks_remaining_ = kSubblocksPerBlock; + EndOfBlock(); + } + } + + // process final samples that don't complete a subblock + if (length > 0) { + if (subblock_active_) + UpdateSubblock(buf, length, decoded_power); + subblock_samples_remaining_ -= length; + } +} + +// Feed some samples into the adaptive system. The samples are guaranteed to not cross a subblock boundary. +// The samples should be processsed (i.e. duty cycle is in the active part) +void AdaptiveGainControl::UpdateSubblock(const std::uint16_t *buf, unsigned length, const double *decoded_power) { + if (decoded_power) { + if (*decoded_power >= burst_loud_threshold_) + ++burst_block_loud_decoded_; + BurstSkip(length); + } else { + BurstUpdate(buf, length); + RangeUpdate(buf, length); + } +} + +// Burst measurement: ignore the next 'length' samples (they are a successfully decoded message) +void AdaptiveGainControl::BurstSkip(unsigned length) { + if (!params_.burst_control) + return; + + // first window + if (length < burst_window_remaining_) { + // partial fill + burst_window_remaining_ -= length; + return; + } + + // skip remainder of first window, dispatch it + BurstEndOfWindow(burst_window_counter_); + length -= burst_window_remaining_; + + // skip remaining windows, dispatch them + unsigned windows = length / samples_per_window_; + unsigned samples = windows * samples_per_window_; + while (windows--) + BurstEndOfWindow(0); + + length -= samples; + + // final partial window + burst_window_counter_ = 0; + burst_window_remaining_ = samples_per_window_ - length; +} + +// Burst measurement: process 'length' samples from 'buf', look for loud bursts; +// the samples might cross burst window boundaries; +// the samples will not cross a block boundary. +void AdaptiveGainControl::BurstUpdate(const std::uint16_t *buf, unsigned length) { + if (!params_.burst_control) + return; + + // first window + if (length < burst_window_remaining_) { + // partial fill + burst_window_counter_ += BurstCountSamples(buf, length); + burst_window_remaining_ -= length; + return; + } + + // complete fill of first partial window + unsigned n = burst_window_remaining_; + unsigned counter = burst_window_counter_ + BurstCountSamples(buf, n); + BurstEndOfWindow(counter); + buf += n; + length -= n; + + // remaining windows + unsigned windows = length / samples_per_window_; + unsigned samples = windows * samples_per_window_; + BurstScanWindows(buf, windows); + buf += samples; + length -= samples; + + // final partial window + burst_window_counter_ = BurstCountSamples(buf, length); + burst_window_remaining_ = samples_per_window_ - length; +} + +// Burst measurement: process 'windows' complete burst windows starting at 'buf'; +// 'buf' is aligned to the start of a burst window +void AdaptiveGainControl::BurstScanWindows(const std::uint16_t *buf, unsigned windows) { + while (windows--) { + unsigned counter = BurstCountSamples(buf, samples_per_window_); + buf += samples_per_window_; + BurstEndOfWindow(counter); + } +} + +// Burst measurement: process 'n' samples from 'buf', look for loud samples; +// the samples are guaranteed not to cross window boundaries; +// return the number of loud samples seen +// (dump1090 uses a SIMD helper here; a scalar loop is fine at 2.08MS/s with duty cycling) +unsigned AdaptiveGainControl::BurstCountSamples(const std::uint16_t *buf, unsigned n) { + unsigned counter = 0; + while (n--) { + if (*buf++ > 46395 /* -3dBFS */) + ++counter; + } + return counter; +} + +// Burst measurement: we reached the end of a burst window with 'counter' +// loud samples seen, handle that window. +void AdaptiveGainControl::BurstEndOfWindow(unsigned counter) { + if (counter > samples_per_window_ / 4) { + // This window is loud, extend any existing run of loud windows + ++burst_runlength_; + } else { + // Quiet window. If we saw a run of loud windows about as long as a + // UAT downlink message, count that as a candidate for an + // over-amplified message that was not decoded. Downlink frames are + // 276 bits (short ADS-B) or 420 bits (long ADS-B) at 1.041667Mbps, + // i.e. ~265us or ~403us, spanning ~7 or ~11 windows of ~40us; accept + // runs of 5..13 windows (~200us..~520us). This excludes ground + // uplink frames (~4.2ms) and continuous interference. (dump1090 uses + // 2..5 windows here, matching Mode S message lengths of 64..120us.) + if (burst_runlength_ >= 5 && burst_runlength_ <= 13) + ++burst_block_loud_undecoded_; + burst_runlength_ = 0; + } +} + +// Noise measurement: process 'length' samples from 'buf'. +// The samples will not cross a block boundary. +void AdaptiveGainControl::RangeUpdate(const std::uint16_t *buf, unsigned length) { + if (!params_.range_control) + return; + + range_radix_counter_ += length; + while (length--) { + // do a very simple radix sort of sample magnitudes + // so we can later find the Nth percentile value + ++range_radix_[buf[0]]; + ++buf; + } +} + +// Noise measurement: we reached the end of a block, update +// our noise estimate +void AdaptiveGainControl::RangeEndOfBlock() { + if (!params_.range_control) + return; + + unsigned n = 0, i = 0; + + // measure Nth percentile magnitude + unsigned count_n = range_radix_counter_ * params_.range_percentile / 100; + while (i < 65536 && n <= count_n) + n += range_radix_[i++]; + std::uint16_t percentile_n = i - 1; + + // maintain an EMA of the Nth percentile + range_smoothed_ = range_smoothed_ * (1 - params_.range_alpha) + percentile_n * params_.range_alpha; + + // reset radix sort for the next block + std::fill(range_radix_.begin(), range_radix_.end(), 0); + range_radix_counter_ = 0; +} + +// Burst measurement: we reached the end of a block, update our burst rate estimate +void AdaptiveGainControl::BurstEndOfBlock() { + if (!params_.burst_control) + return; + + // scale rates based on the actual duty cycle fraction + // (e.g. if we are only inspecting 2/5 of samples, then scale the rate by 5/2) + double scale = (double)kSubblocksPerBlock / dutycycle_N_; + + // maintain an EMA of the number of undecoded loud bursts seen per block + report_loud_undecoded_ += burst_block_loud_undecoded_; + burst_loud_undecoded_smoothed_ = burst_loud_undecoded_smoothed_ * (1 - params_.burst_alpha) + scale * burst_block_loud_undecoded_ * params_.burst_alpha; + burst_block_loud_undecoded_ = 0; + + // maintain an EMA of the number of decoded, but loud, messages seen per block + report_loud_decoded_ += burst_block_loud_decoded_; + burst_loud_decoded_smoothed_ = burst_loud_decoded_smoothed_ * (1 - params_.burst_alpha) + scale * burst_block_loud_decoded_ * params_.burst_alpha; + burst_block_loud_decoded_ = 0; +} + +// Adaptive gain: we reached a block boundary. Update measurements and act on them. +void AdaptiveGainControl::EndOfBlock() { + RangeEndOfBlock(); + BurstEndOfBlock(); + + ControlUpdate(); + + if (report_timer_ > 0 && --report_timer_ == 0) { + PeriodicReport(); + report_timer_ = kReportInterval; + } +} + +// Periodic status line so that convergence is observable in logs +// (replaces dump1090's stats subsystem) +void AdaptiveGainControl::PeriodicReport() { + fprintf(stderr, "adaptive: gain %.1fdB (step %d)", StepDb(current_step_), current_step_); + if (params_.range_control) { + if (range_smoothed_ > 0) + fprintf(stderr, ", noise floor %.1fdBFS, dynamic range %.1fdB", 20 * std::log10(range_smoothed_ / 65536.0), -20 * std::log10(range_smoothed_ / 65536.0)); + else + fprintf(stderr, ", noise floor not yet measured"); + } + if (params_.burst_control) + fprintf(stderr, ", loud undecoded %u, loud decoded %u", report_loud_undecoded_, report_loud_decoded_); + fprintf(stderr, ", gain changes %u (last %us)\n", report_gain_changes_, kReportInterval); + + report_gain_changes_ = 0; + report_loud_undecoded_ = 0; + report_loud_decoded_ = 0; +} + +void AdaptiveGainControl::ControlUpdate() { + // votes for what to do with the gain + // "gain_not_up" overlaps somewhat with "gain_down", but they are not identical; + // burst control may want to prevent gain from increasing, but not necessarily + // decrease gain. + + bool gain_up = false; + const char *gain_up_reason = nullptr; + bool gain_down = false; + const char *gain_down_reason = nullptr; + bool gain_not_up = false; + + int current_gain = current_step_; + + if (burst_change_timer_) + --burst_change_timer_; + if (range_change_timer_ > 0) + --range_change_timer_; + if (range_rescan_timer_ > 0) + --range_rescan_timer_; + + if (params_.burst_control && !burst_change_timer_) { + if (burst_loud_undecoded_smoothed_ > params_.burst_loud_rate) { + burst_quiet_blocks_ = 0; + ++burst_loud_blocks_; + } else if (burst_loud_decoded_smoothed_ < params_.burst_quiet_rate) { + burst_loud_blocks_ = 0; + ++burst_quiet_blocks_; + } else { + burst_loud_blocks_ = 0; + burst_quiet_blocks_ = 0; + } + + if (burst_loud_blocks_ >= params_.burst_loud_runlength) { + // we need to reduce gain (further) + gain_down = gain_not_up = true; + gain_down_reason = "high rate of loud undecoded messages"; + + // if we're currently doing a downward scan, reducing gain further may confuse it; + // stop that scan and restart it once we are no longer in a reduced-gain state + if (range_state_ == RangeState::SCAN_DOWN || range_state_ == RangeState::RESCAN_DOWN) { + range_state_ = RangeState::SCAN_IDLE; + range_rescan_timer_ = 0; + } + } else if (burst_quiet_blocks_ < params_.burst_quiet_runlength) { + // we're OK at the current gain, but should not increase it + gain_not_up = true; + } else if (current_gain < range_gain_limit_) { + // we're OK at the current gain, and can increase gain to the previously discovered + // dynamic range limit + gain_up = true; + gain_up_reason = "low loud message rate and gain below dynamic range limit"; + } + } + + if (params_.range_control && !range_change_timer_) { + double available_range = -20 * std::log10(range_smoothed_ / 65536.0); + // allow the gain limit to increase if this gain setting is acceptable + // (decreasing the limit is done separately depending on the current state as we make slightly different decisions in IDLE + // to provide hysteresis) + if (available_range >= params_.range_target && current_gain > range_gain_limit_) { + range_gain_limit_ = current_gain; + } + switch (range_state_) { + case RangeState::SCAN_UP: + case RangeState::RESCAN_UP: + if (available_range < params_.range_target) { + // Current gain fails to meet our target. Switch to downward scanning. + fprintf(stderr, "adaptive: available dynamic range (%.1fdB) < required dynamic range (%.1fdB), switching to downward scan\n", available_range, params_.range_target); + gain_down = gain_not_up = true; + gain_down_reason = "probing dynamic range gain lower bound"; + range_state_ = (range_state_ == RangeState::RESCAN_UP ? RangeState::RESCAN_DOWN : RangeState::SCAN_DOWN); + if (range_gain_limit_ >= current_gain) { + range_gain_limit_ = current_gain - 1; + } + break; + } + + if (current_step_ >= gain_max_) { + // We have reached our upper gain limit + fprintf(stderr, "adaptive: reached upper gain limit, halting dynamic range scan here\n"); + range_state_ = RangeState::SCAN_IDLE; + range_rescan_timer_ = params_.range_rescan_delay; + break; + } + + // This gain step is OK and we have more to try, try the next gain step up. + // (But if burst detection has inhibited increasing gain, don't do anything yet, just try again next block) + if (!gain_not_up) { + fprintf(stderr, "adaptive: available dynamic range (%.1fdB) >= required dynamic range (%.1fdB), continuing upward scan\n", available_range, params_.range_target); + gain_up = true; + gain_up_reason = "probing dynamic range gain upper bound"; + } + break; + + case RangeState::SCAN_DOWN: + case RangeState::RESCAN_DOWN: + if (available_range >= params_.range_target) { + // Current gain meets our target; we are done with the scan. + fprintf(stderr, "adaptive: available dynamic range (%.1fdB) >= required dynamic range (%.1fdB), stopping downwards scan here\n", available_range, params_.range_target); + range_state_ = RangeState::SCAN_IDLE; + // Note: quirk ported as-is from dump1090's adaptive.c for behavior parity: + // the state is updated before this comparison, so the rescan delay is + // always selected here regardless of whether this was a SCAN or RESCAN. + range_rescan_timer_ = (range_state_ == RangeState::SCAN_DOWN ? params_.range_scan_delay : params_.range_rescan_delay); + break; + } + + if (range_gain_limit_ >= current_gain) { + range_gain_limit_ = current_gain - 1; + } + + if (current_step_ <= gain_min_) { + fprintf(stderr, "adaptive: reached lower gain limit, halting dynamic range scan here\n"); + range_state_ = RangeState::SCAN_IDLE; + range_rescan_timer_ = params_.range_rescan_delay; + break; + } + + // This gain step is too loud and we have more to try, try the next gain step down + fprintf(stderr, "adaptive: available dynamic range (%.1fdB) < required dynamic range (%.1fdB), continuing downwards scan\n", available_range, params_.range_target); + gain_down = gain_not_up = true; + gain_down_reason = "probing dynamic range gain lower bound"; + break; + + case RangeState::SCAN_IDLE: + // Look for increased noise that could be compensated for by decreasing gain. + // Do this even if we're waiting to rescan or if burst control is also active + if (available_range + gain_down_db_ / 2 < params_.range_target && current_step_ > gain_min_) { + fprintf(stderr, "adaptive: available dynamic range (%.1fdB) + half gain step down (%.1fdB) < required dynamic range (%.1fdB), starting downward scan\n", available_range, gain_down_db_ / 2, params_.range_target); + if (range_gain_limit_ >= current_gain) { + range_gain_limit_ = current_gain - 1; + } + range_state_ = RangeState::SCAN_DOWN; + gain_down = gain_not_up = true; + gain_down_reason = "dynamic range fell below target value"; + break; + } + + // Infrequently consider increasing gain to handle the case where we've selected a too-low gain where the noise floor is dominated by noise unrelated to the gain setting. + // But don't do this while burst control is preventing gain increases. + if (!range_rescan_timer_ && !gain_not_up) { + if (available_range >= params_.range_target && current_step_ < gain_max_) { + fprintf(stderr, "adaptive: start periodic scan for acceptable dynamic range at increased gain\n"); + gain_up = true; + gain_up_reason = "periodic re-probing of dynamic range gain upper bound"; + range_state_ = RangeState::RESCAN_UP; + break; + } + + // Nothing to do for a while. + range_rescan_timer_ = params_.range_rescan_delay; + } + + break; + + default: + fprintf(stderr, "adaptive: in a weird state (%u), trying to fix it\n", (unsigned)range_state_); + range_state_ = RangeState::SCAN_IDLE; + range_rescan_timer_ = params_.range_rescan_delay; + break; + } + } + + // now actually perform any gain changes + + if (gain_down) + DecreaseGain(gain_down_reason); + else if (gain_up && !gain_not_up) + IncreaseGain(gain_up_reason); +} diff --git a/adaptive.h b/adaptive.h new file mode 100644 index 0000000..0f97217 --- /dev/null +++ b/adaptive.h @@ -0,0 +1,176 @@ +// -*- c++ -*- + +// adaptive.h: adaptive gain control +// +// Ported from dump1090's adaptive.c / adaptive.h +// (https://github.com/flightaware/dump1090) +// +// Copyright (c) 2021 FlightAware, LLC +// +// This file is free software: you may copy, redistribute and/or modify it +// under the terms of the GNU General Public License as published by the +// Free Software Foundation, either version 2 of the License, or (at your +// option) any later version. +// +// This file is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#ifndef DUMP978_ADAPTIVE_H +#define DUMP978_ADAPTIVE_H + +#include +#include +#include + +namespace flightaware { + namespace uat { + // Tunables for adaptive gain control; defaults match dump1090. + struct AdaptiveParams { + bool burst_control = false; // enable burst / overload control? + bool range_control = false; // enable dynamic range control? + double min_gain_db = 0; // never step gain below this + double max_gain_db = 99999; // never step gain above this + double duty_cycle = 0.5; // fraction of samples to inspect + double burst_alpha = 2.0 / (5 + 1); // EMA factor for burst rates + unsigned burst_change_delay = 5; // blocks to wait after a gain change before acting on burst measurements + double burst_loud_rate = 5.0; // undecoded loud burst rate that provokes a gain decrease + unsigned burst_loud_runlength = 10; // consecutive loud blocks needed before decreasing gain + double burst_quiet_rate = 5.0; // loud decoded message rate that must not be exceeded before increasing gain + unsigned burst_quiet_runlength = 10; // consecutive quiet blocks needed before increasing gain + double range_alpha = 2.0 / (5 + 1); // EMA factor for the noise floor estimate + unsigned range_percentile = 40; // percentile of sample magnitude to use as the noise floor + double range_target = 0; // target dynamic range in dB; <= 0 means autoselect + unsigned range_change_delay = 10; // blocks to wait after a gain change before acting on noise measurements + unsigned range_scan_delay = 300; // blocks to wait before rescanning after an abortive scan + unsigned range_rescan_delay = 3600; // blocks to wait between periodic rescans for increased gain + }; + + // Abstract interface for changing SDR gain; implemented by SoapySampleSource. + // Kept free of SoapySDR types so that adaptive gain control is unit-testable. + class SdrGainControl { + public: + typedef std::shared_ptr Pointer; + + virtual ~SdrGainControl() {} + + // Return the device's current gain in dB + virtual double GetGainDb() = 0; + + // Try to set the device gain to `db` dB; return the gain the device reports afterwards + virtual double SetGainDb(double db) = 0; + + // Return the device's overall gain range; `step` may be zero if the device + // does not report a step size + virtual void GetGainLimits(double &min, double &max, double &step) = 0; + }; + + // Software adaptive gain control, a port of dump1090's adaptive.c. + // + // Feed magnitude samples via Update(); the controller measures the noise + // floor and loud undecoded bursts, and steps the SDR gain up or down at + // ~1 second intervals to hold a target dynamic range and back off from + // overload. Not thread-safe: all calls must come from the same thread + // (the SDR rx/demodulator thread). + class AdaptiveGainControl { + public: + typedef std::shared_ptr Pointer; + + static Pointer Create(const AdaptiveParams ¶ms, SdrGainControl::Pointer sdr, unsigned sample_rate) { return Pointer(new AdaptiveGainControl(params, sdr, sample_rate)); } + + // false if adaptive control was requested but the SDR's gain is not controllable + bool Enabled() const { return params_.burst_control || params_.range_control; } + + // Feed samples into the adaptive system. Any number of samples might be passed in. + // `buf` holds one uint16 magnitude per sample (65536 = full scale). + // `decoded_power` is nullptr for samples that were not part of a successfully + // decoded message; for decoded message samples, it points to the mean power of + // the message (linear, 1.0 = full scale, i.e. mean magnitude-squared). + void Update(const std::uint16_t *buf, unsigned length, const double *decoded_power); + + private: + AdaptiveGainControl(const AdaptiveParams ¶ms, SdrGainControl::Pointer sdr, unsigned sample_rate); + + int GainStepCount() const { return (int)gain_steps_.size(); } + double StepDb(int step) const; + int NearestStep(double db) const; + bool SetGainStep(int step, const char *why); + void GainChanged(); + void IncreaseGain(const char *why); + void DecreaseGain(const char *why); + + void UpdateSubblock(const std::uint16_t *buf, unsigned length, const double *decoded_power); + void EndOfBlock(); + void ControlUpdate(); + void PeriodicReport(); + + void BurstUpdate(const std::uint16_t *buf, unsigned length); + void BurstSkip(unsigned length); + unsigned BurstCountSamples(const std::uint16_t *buf, unsigned n); + void BurstScanWindows(const std::uint16_t *buf, unsigned windows); + void BurstEndOfWindow(unsigned counter); + void BurstEndOfBlock(); + + void RangeUpdate(const std::uint16_t *buf, unsigned length); + void RangeEndOfBlock(); + + AdaptiveParams params_; + SdrGainControl::Pointer sdr_; + + // gain step table and limits + std::vector gain_steps_; // dB value for each gain step + int current_step_ = 0; // gain step we last set (logical state; the hardware may quantize) + int gain_min_ = 0; // lowest gain step we will use + int gain_max_ = 0; // highest gain step we will use + + // gain steps relative to current gain + double gain_up_db_ = 0; + double gain_down_db_ = 0; + + // block handling + unsigned samples_per_window_ = 0; // samples per burst window (~40us) + unsigned samples_per_subblock_ = 0; // samples per subblock (~50ms) + unsigned subblocks_remaining_ = 0; // subblocks remaining in the current block + + unsigned dutycycle_N_ = 0; // subblock duty cycle numerator N (out of D = kSubblocksPerBlock) + unsigned dutycycle_counter_ = 0; // subblock duty cycle counter (modulo D) + bool subblock_active_ = false; // is the current subblock active i.e. samples should be processed, not skipped? + unsigned subblock_samples_remaining_ = 0; // samples remaining in the current subblock + + // burst handling + unsigned burst_window_remaining_ = 0; // samples remaining in the current burst window + unsigned burst_window_counter_ = 0; // loud samples seen in current burst window + unsigned burst_runlength_ = 0; // consecutive loud burst windows seen + unsigned burst_block_loud_undecoded_ = 0; // loud undecoded bursts seen in this block so far + unsigned burst_block_loud_decoded_ = 0; // loud decoded messages seen in this block so far + double burst_loud_undecoded_smoothed_ = 0; // smoothed rate of loud misdecodes per block + double burst_loud_decoded_smoothed_ = 0; // smoothed rate of loud successful decodes per block + unsigned burst_change_timer_ = 0; // countdown inhibiting control after changing gain + double burst_loud_threshold_ = 0; // current signal level threshold for a "loud decode" + unsigned burst_loud_blocks_ = 0; // consecutive blocks with loud rate + unsigned burst_quiet_blocks_ = 0; // consecutive blocks with quiet rate + + // noise floor measurement (adaptive dynamic range) + std::vector range_radix_; // radix-sort buckets for current block + unsigned range_radix_counter_ = 0; // sum of all radix-sort buckets (= number of samples sorted) + double range_smoothed_ = 0; // smoothed noise floor estimate (u16 magnitude units) + enum class RangeState { SCAN_IDLE, SCAN_UP, SCAN_DOWN, RESCAN_UP, RESCAN_DOWN }; + RangeState range_state_ = RangeState::SCAN_UP; + unsigned range_change_timer_ = 0; // countdown inhibiting control after changing gain + unsigned range_rescan_timer_ = 0; // countdown to next upwards gain reprobe + int range_gain_limit_ = 0; // probed maximum gain step with acceptable dynamic range + + // periodic reporting (replaces dump1090's stats subsystem) + unsigned report_timer_ = 0; // blocks until the next periodic report + unsigned report_gain_changes_ = 0; // gain changes since the last report + unsigned report_loud_undecoded_ = 0; // loud undecoded bursts since the last report + unsigned report_loud_decoded_ = 0; // loud decoded messages since the last report + }; + }; // namespace flightaware::uat +}; // namespace flightaware + +#endif diff --git a/adaptive_tests.cc b/adaptive_tests.cc new file mode 100644 index 0000000..cdbb112 --- /dev/null +++ b/adaptive_tests.cc @@ -0,0 +1,241 @@ +// adaptive_tests.cc: tests for adaptive gain control +// +// Copyright (c) 2021 FlightAware, LLC +// +// This file is free software: you may copy, redistribute and/or modify it +// under the terms of the GNU General Public License as published by the +// Free Software Foundation, either version 2 of the License, or (at your +// option) any later version. +// +// This file is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#include "adaptive.h" + +#include +#include +#include + +using namespace flightaware::uat; + +static int failures = 0; + +#define CHECK(cond) \ + do { \ + if (!(cond)) { \ + std::fprintf(stderr, "FAIL at line %d: %s\n", __LINE__, #cond); \ + ++failures; \ + } \ + } while (0) + +// A fake SDR that behaves like an RTLSDR seen through SoapySDR: +// gain range 0 .. 49.6 dB, no step size reported. +class FakeGainControl : public SdrGainControl { + public: + FakeGainControl() : gain_(49.6) {} + + double GetGainDb() override { return gain_; } + + double SetGainDb(double db) override { + gain_ = db; + history.push_back(db); + return db; + } + + void GetGainLimits(double &min, double &max, double &step) override { + min = 0; + max = 49.6; + step = 0; + } + + std::vector history; + + private: + double gain_; +}; + +// Use a low sample rate so the tests can feed many ~1 second blocks quickly. +// The adaptive block structure scales with the sample rate (window = rate/25000). +static const unsigned kTestRate = 250000; +static const unsigned kSamplesPerWindow = kTestRate / 25000; +static const unsigned kSamplesPerSubblock = kSamplesPerWindow * 1250; +static const unsigned kSubblocksPerBlock = 20; + +// Feed `blocks` whole blocks built from copies of `subblock` +static void FeedBlocks(AdaptiveGainControl::Pointer adaptive, const std::vector &subblock, unsigned blocks, const double *decoded_power = nullptr) { + for (unsigned b = 0; b < blocks; ++b) { + for (unsigned s = 0; s < kSubblocksPerBlock; ++s) { + adaptive->Update(subblock.data(), subblock.size(), decoded_power); + } + } +} + +static std::vector ConstantSubblock(std::uint16_t magnitude) { return std::vector(kSamplesPerSubblock, magnitude); } + +// A subblock containing repeated full-scale bursts of the given window count +// separated by silence; at 8 windows (~320us), what an overloaded receiver +// failing to decode UAT downlink messages looks like +static std::vector BurstySubblock(unsigned burst_windows = 8) { + std::vector subblock(kSamplesPerSubblock, 0); + const unsigned period = kSamplesPerWindow * 40; + for (unsigned i = 0; i < subblock.size(); ++i) { + if (i % period < kSamplesPerWindow * burst_windows) + subblock[i] = 65535; + } + return subblock; +} + +// uint16 magnitude for a given dBFS level +static std::uint16_t MagnitudeAtDbfs(double dbfs) { return (std::uint16_t)std::round(std::pow(10.0, dbfs / 20.0) * 65536.0); } + +// The gain limits from --adaptive-min-gain/--adaptive-max-gain must be honored, +// including the initial constraining of the starting gain. +static void TestGainLimits() { + std::fprintf(stderr, "--- TestGainLimits\n"); + auto fake = std::make_shared(); + AdaptiveParams params; + params.range_control = true; + params.min_gain_db = 20; + params.max_gain_db = 40; + auto adaptive = AdaptiveGainControl::Create(params, fake, kTestRate); + + CHECK(adaptive->Enabled()); + // the fake device starts at 49.6dB; init must pull it down to <= 40dB + CHECK(!fake->history.empty()); + CHECK(fake->GetGainDb() <= 40.0); + CHECK(fake->GetGainDb() >= 20.0); +} + +// With plenty of dynamic range available, the range scan should hold the gain +// at the upper limit and never step down. +static void TestRangeQuiet() { + std::fprintf(stderr, "--- TestRangeQuiet\n"); + auto fake = std::make_shared(); + AdaptiveParams params; + params.range_control = true; + params.duty_cycle = 1.0; + auto adaptive = AdaptiveGainControl::Create(params, fake, kTestRate); + + // constant -40dBFS noise floor; default target is 0.6 * 49.6 = 29.8dB + FeedBlocks(adaptive, ConstantSubblock(MagnitudeAtDbfs(-40.0)), 30); + + CHECK(fake->GetGainDb() > 49.0); +} + +// With too little dynamic range, the range scan should walk the gain down +// until it hits the configured lower limit. +static void TestRangeScanDown() { + std::fprintf(stderr, "--- TestRangeScanDown\n"); + auto fake = std::make_shared(); + AdaptiveParams params; + params.range_control = true; + params.duty_cycle = 1.0; + params.min_gain_db = 40; + auto adaptive = AdaptiveGainControl::Create(params, fake, kTestRate); + + // constant -20dBFS noise floor: only 20dB of dynamic range available, + // below the ~29.8dB target no matter the gain, so the scan runs to the lower limit + FeedBlocks(adaptive, ConstantSubblock(MagnitudeAtDbfs(-20.0)), 150); + + CHECK(!fake->history.empty()); + CHECK(std::fabs(fake->GetGainDb() - 40.0) < 0.5); + // the scan must step down one step at a time + for (std::size_t i = 1; i < fake->history.size(); ++i) { + CHECK(fake->history[i] < fake->history[i - 1]); + CHECK(fake->history[i - 1] - fake->history[i] < 1.5); + } +} + +// Loud, message-length, undecoded bursts should push the gain down; after the +// bursts stop, the gain should recover towards the previously known-good limit. +static void TestBurstBackoffAndRecovery() { + std::fprintf(stderr, "--- TestBurstBackoffAndRecovery\n"); + auto fake = std::make_shared(); + AdaptiveParams params; + params.burst_control = true; + auto adaptive = AdaptiveGainControl::Create(params, fake, kTestRate); + + FeedBlocks(adaptive, BurstySubblock(), 40); + + CHECK(!fake->history.empty()); + double backed_off = fake->GetGainDb(); + CHECK(backed_off < 49.0); + + // quiet input: burst control should raise the gain back towards the limit + FeedBlocks(adaptive, ConstantSubblock(0), 40); + CHECK(fake->GetGainDb() > backed_off); +} + +// Loud bursts much shorter than a UAT downlink message (e.g. DME pulses or +// Mode S messages from the adjacent band) must not trigger a backoff. +static void TestBurstWrongLength() { + std::fprintf(stderr, "--- TestBurstWrongLength\n"); + auto fake = std::make_shared(); + AdaptiveParams params; + params.burst_control = true; + auto adaptive = AdaptiveGainControl::Create(params, fake, kTestRate); + + FeedBlocks(adaptive, BurstySubblock(2), 40); + + CHECK(fake->history.empty()); + CHECK(fake->GetGainDb() > 49.0); +} + +// The same bursts fed as successfully decoded messages must not trigger a backoff. +static void TestBurstDecodedSkip() { + std::fprintf(stderr, "--- TestBurstDecodedSkip\n"); + auto fake = std::make_shared(); + AdaptiveParams params; + params.burst_control = true; + auto adaptive = AdaptiveGainControl::Create(params, fake, kTestRate); + + // decoded messages with a quiet mean power (below the loud-decode threshold) + double decoded_power = 1e-4; + FeedBlocks(adaptive, BurstySubblock(), 40, &decoded_power); + + CHECK(fake->history.empty()); + CHECK(fake->GetGainDb() > 49.0); +} + +// An uncontrollable SDR (no gain range) must disable adaptive control cleanly. +class FixedGainControl : public SdrGainControl { + public: + double GetGainDb() override { return 0; } + double SetGainDb(double db) override { return 0; } + void GetGainLimits(double &min, double &max, double &step) override { min = max = step = 0; } +}; + +static void TestUncontrollableSdr() { + std::fprintf(stderr, "--- TestUncontrollableSdr\n"); + AdaptiveParams params; + params.range_control = true; + params.burst_control = true; + auto adaptive = AdaptiveGainControl::Create(params, std::make_shared(), kTestRate); + CHECK(!adaptive->Enabled()); +} + +int main(int argc, char **argv) { + (void)argc; + (void)argv; + + TestGainLimits(); + TestRangeQuiet(); + TestRangeScanDown(); + TestBurstBackoffAndRecovery(); + TestBurstWrongLength(); + TestBurstDecodedSkip(); + TestUncontrollableSdr(); + + if (failures) { + std::fprintf(stderr, "%d test(s) FAILED\n", failures); + return 1; + } + + std::fprintf(stderr, "all tests passed\n"); + return 0; +} diff --git a/convert.cc b/convert.cc index 53a314b..d10399d 100644 --- a/convert.cc +++ b/convert.cc @@ -32,6 +32,12 @@ static inline std::uint16_t scaled_atan(double x) { static inline double magsq(double i, double q) { return i * i + q * q; } +// convert a magnitude-squared value (1.0 = full scale) to a uint16 magnitude (65536 = full scale, clamped) +static inline std::uint16_t scaled_mag16(double msq) { + double mag = std::round(std::sqrt(msq) * 65536.0); + return mag > 65535 ? 65535 : (std::uint16_t)mag; +} + SampleConverter::Pointer SampleConverter::Create(SampleFormat format) { switch (format) { case SampleFormat::CU8: @@ -60,6 +66,7 @@ CU8Converter::CU8Converter() : SampleConverter(SampleFormat::CU8) { u.iq[1] = q; lookup_phase_[u.iq16] = scaled_atan2(d_q, d_i); lookup_magsq_[u.iq16] = magsq(d_i, d_q); + lookup_mag16_[u.iq16] = scaled_mag16(lookup_magsq_[u.iq16]); } } } @@ -110,6 +117,29 @@ void CU8Converter::ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iterat } } +void CU8Converter::ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) { + const cu8_alias *in_iq = reinterpret_cast(&*begin); + + // unroll the loop + const auto n = std::distance(begin, end) / 2; + const auto n8 = n / 8; + const auto n7 = n & 7; + + for (auto i = 0; i < n8; ++i, in_iq += 8) { + *out++ = lookup_mag16_[in_iq[0].iq16]; + *out++ = lookup_mag16_[in_iq[1].iq16]; + *out++ = lookup_mag16_[in_iq[2].iq16]; + *out++ = lookup_mag16_[in_iq[3].iq16]; + *out++ = lookup_mag16_[in_iq[4].iq16]; + *out++ = lookup_mag16_[in_iq[5].iq16]; + *out++ = lookup_mag16_[in_iq[6].iq16]; + *out++ = lookup_mag16_[in_iq[7].iq16]; + } + for (auto i = 0; i < n7; ++i, ++in_iq) { + *out++ = lookup_mag16_[in_iq[0].iq16]; + } +} + CS8Converter::CS8Converter() : SampleConverter(SampleFormat::CS8_) { cs8_alias u; @@ -122,6 +152,7 @@ CS8Converter::CS8Converter() : SampleConverter(SampleFormat::CS8_) { u.iq[1] = q; lookup_phase_[u.iq16] = scaled_atan2(d_q, d_i); lookup_magsq_[u.iq16] = magsq(d_i, d_q); + lookup_mag16_[u.iq16] = scaled_mag16(lookup_magsq_[u.iq16]); } } } @@ -172,6 +203,29 @@ void CS8Converter::ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iterat } } +void CS8Converter::ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) { + auto in_iq = reinterpret_cast(&*begin); + + // unroll the loop + const auto n = std::distance(begin, end) / 2; + const auto n8 = n / 8; + const auto n7 = n & 7; + + for (auto i = 0; i < n8; ++i, in_iq += 8) { + *out++ = lookup_mag16_[in_iq[0].iq16]; + *out++ = lookup_mag16_[in_iq[1].iq16]; + *out++ = lookup_mag16_[in_iq[2].iq16]; + *out++ = lookup_mag16_[in_iq[3].iq16]; + *out++ = lookup_mag16_[in_iq[4].iq16]; + *out++ = lookup_mag16_[in_iq[5].iq16]; + *out++ = lookup_mag16_[in_iq[6].iq16]; + *out++ = lookup_mag16_[in_iq[7].iq16]; + } + for (auto i = 0; i < n7; ++i, ++in_iq) { + *out++ = lookup_mag16_[in_iq[0].iq16]; + } +} + CS16HConverter::CS16HConverter() : SampleConverter(SampleFormat::CS16H) { // atan lookup, positive values only, 8-bit fixed point covering 0.0 .. 256.0 for (std::size_t i = 0; i < lookup_atan_.size(); ++i) { @@ -269,6 +323,15 @@ void CS16HConverter::ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iter } } +void CS16HConverter::ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) { + auto in_iq = reinterpret_cast(&*begin); + + const auto n = std::distance(begin, end) / 4; + for (auto i = 0; i < n; ++i, in_iq += 2) { + *out++ = scaled_mag16(magsq(in_iq[1], in_iq[0]) / 32768.0 / 32768.0); + } +} + void CF32HConverter::ConvertPhase(Bytes::const_iterator begin, Bytes::const_iterator end, PhaseBuffer::iterator out) { auto in_iq = reinterpret_cast(&*begin); @@ -314,3 +377,12 @@ void CF32HConverter::ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iter *out++ = magsq(in_iq[1], in_iq[0]); } } + +void CF32HConverter::ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) { + auto in_iq = reinterpret_cast(&*begin); + + const auto n = std::distance(begin, end) / 8; + for (auto i = 0; i < n; ++i, in_iq += 2) { + *out++ = scaled_mag16(magsq(in_iq[1], in_iq[0])); + } +} diff --git a/convert.h b/convert.h index 1b8aaa3..5e2a107 100644 --- a/convert.h +++ b/convert.h @@ -57,6 +57,12 @@ namespace flightaware { // samples (trailing partial samples are ignored, not buffered). virtual void ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) = 0; + // Read samples from `begin` .. `end` and write one uint16 magnitude value + // per sample to `out` (65536 = full scale, clamped to 65535). The input + // buffer should contain an integral number of samples (trailing partial + // samples are ignored, not buffered). + virtual void ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) = 0; + SampleFormat Format() const { return format_; } unsigned BytesPerSample() const { return bytes_per_sample_; } @@ -74,6 +80,7 @@ namespace flightaware { void ConvertPhase(Bytes::const_iterator begin, Bytes::const_iterator end, PhaseBuffer::iterator out) override; void ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; + void ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; private: union cu8_alias { @@ -83,6 +90,7 @@ namespace flightaware { std::array lookup_phase_; std::array lookup_magsq_; + std::array lookup_mag16_; }; class CS8Converter : public SampleConverter { @@ -91,6 +99,7 @@ namespace flightaware { void ConvertPhase(Bytes::const_iterator begin, Bytes::const_iterator end, PhaseBuffer::iterator out) override; void ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; + void ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; private: union cs8_alias { @@ -100,6 +109,7 @@ namespace flightaware { std::array lookup_phase_; std::array lookup_magsq_; + std::array lookup_mag16_; }; class CS16HConverter : public SampleConverter { @@ -107,6 +117,7 @@ namespace flightaware { CS16HConverter(); void ConvertPhase(Bytes::const_iterator begin, Bytes::const_iterator end, PhaseBuffer::iterator out) override; void ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; + void ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; private: std::uint16_t TableAtan(std::uint32_t r); @@ -119,6 +130,7 @@ namespace flightaware { CF32HConverter() : SampleConverter(SampleFormat::CF32H) {} void ConvertPhase(Bytes::const_iterator begin, Bytes::const_iterator end, PhaseBuffer::iterator out) override; void ConvertMagSq(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; + void ConvertMag16(Bytes::const_iterator begin, Bytes::const_iterator end, std::vector::iterator out) override; }; }; // namespace flightaware::uat }; // namespace flightaware diff --git a/demodulator.cc b/demodulator.cc index 6d4c6ba..bca801f 100644 --- a/demodulator.cc +++ b/demodulator.cc @@ -11,7 +11,7 @@ using namespace flightaware::uat; -SingleThreadReceiver::SingleThreadReceiver(SampleFormat format) : converter_(SampleConverter::Create(format)), demodulator_(new TwoMegDemodulator()) {} +SingleThreadReceiver::SingleThreadReceiver(SampleFormat format, AdaptiveGainControl::Pointer adaptive) : converter_(SampleConverter::Create(format)), demodulator_(new TwoMegDemodulator()), adaptive_(adaptive) {} // Handle samples in 'buffer' by: // converting them to a phase buffer @@ -44,6 +44,14 @@ void SingleThreadReceiver::HandleSamples(std::uint64_t timestamp, Bytes::const_i converter_->ConvertPhase(samples_.begin(), samples_.begin() + total_bytes, phase_.begin()); auto messages = demodulator_->Demodulate(phase_.begin(), phase_.begin() + total_samples); + // decoded message extents & mean power, for the adaptive gain controller + struct DecodedSpan { + std::size_t begin; + std::size_t end; + double power; + }; + std::vector decoded_spans; + if (!messages.empty()) { SharedMessageVector dispatch = std::make_shared(); dispatch->reserve(messages.size()); @@ -61,15 +69,71 @@ void SingleThreadReceiver::HandleSamples(std::uint64_t timestamp, Bytes::const_i total_power += m; } - auto rssi = (total_power == 0 ? -1000 : 10 * std::log10(total_power / magsq.size())); + auto mean_power = total_power / magsq.size(); + auto rssi = (total_power == 0 ? -1000 : 10 * std::log10(mean_power)); std::uint64_t message_timestamp = timestamp - (1000 * previous_samples / 2083333) + (1000 * std::distance(phase_.cbegin(), message.begin) / 2083333); dispatch->emplace_back(std::move(message.payload), message_timestamp, message.corrected_errors, rssi); + + if (adaptive_) { + decoded_spans.push_back(DecodedSpan{(std::size_t)std::distance(phase_.cbegin(), message.begin), (std::size_t)std::distance(phase_.cbegin(), message.end), mean_power}); + } } DispatchMessages(dispatch); } + // Feed sample magnitudes to the adaptive gain controller. The last tail_size + // samples are re-presented on the next call, so to feed every sample exactly + // once we stop short of them here and pick them up next time; by then, any + // message starting in that region has been demodulated, so decoded spans are + // known when their samples are fed. Decoded message samples are fed with the + // message's mean power; everything else is fed as undecoded samples. A + // decoded message clipped at the feed boundary has its remaining samples + // carried over and fed as decoded samples on the next call, so that burst + // detection never sees part of a decoded message as an undecoded burst. + if (adaptive_ && adaptive_->Enabled()) { + const std::size_t tail = demodulator_->NumTrailingSamples(); + const std::size_t feed_end = ((std::size_t)total_samples > tail) ? (std::size_t)total_samples - tail : 0; + if (feed_end > 0) { + if (mag16_.size() < feed_end) { + mag16_.resize(feed_end); + } + converter_->ConvertMag16(samples_.begin(), samples_.begin() + feed_end * converter_->BytesPerSample(), mag16_.begin()); + + std::size_t cursor = 0; + if (adaptive_carry_samples_ > 0) { + // remainder of a decoded message clipped at the previous feed + // boundary; it now sits at the start of the buffer + auto carry_end = std::min(adaptive_carry_samples_, feed_end); + adaptive_->Update(&mag16_[0], carry_end, &adaptive_carry_power_); + cursor = carry_end; + adaptive_carry_samples_ -= carry_end; + } + for (const auto &span : decoded_spans) { + if (span.begin >= feed_end) { + break; + } + auto span_begin = std::max(span.begin, cursor); + if (span_begin > cursor) { + adaptive_->Update(&mag16_[cursor], span_begin - cursor, nullptr); + } + auto span_end = std::min(span.end, feed_end); + if (span_end > span_begin) { + adaptive_->Update(&mag16_[span_begin], span_end - span_begin, &span.power); + cursor = span_end; + } + if (span.end > feed_end) { + adaptive_carry_samples_ = span.end - feed_end; + adaptive_carry_power_ = span.power; + } + } + if (cursor < feed_end) { + adaptive_->Update(&mag16_[cursor], feed_end - cursor, nullptr); + } + } + } + // preserve the tail of the sample buffer for next time const auto tail_size = demodulator_->NumTrailingSamples(); if (total_samples > tail_size) { diff --git a/demodulator.h b/demodulator.h index 47f0020..9ce136d 100644 --- a/demodulator.h +++ b/demodulator.h @@ -10,6 +10,7 @@ #include #include +#include "adaptive.h" #include "common.h" #include "convert.h" #include "fec.h" @@ -57,18 +58,25 @@ namespace flightaware { class SingleThreadReceiver : public Receiver { public: - SingleThreadReceiver(SampleFormat format); + SingleThreadReceiver(SampleFormat format, AdaptiveGainControl::Pointer adaptive = nullptr); void HandleSamples(std::uint64_t timestamp, Bytes::const_iterator begin, Bytes::const_iterator end) override; private: SampleConverter::Pointer converter_; std::unique_ptr demodulator_; + AdaptiveGainControl::Pointer adaptive_; Bytes samples_; std::size_t saved_samples_ = 0; PhaseBuffer phase_; + std::vector mag16_; + + // remainder of a decoded message clipped at the adaptive feed + // boundary, to be fed as decoded samples on the next call + std::size_t adaptive_carry_samples_ = 0; + double adaptive_carry_power_ = 0; }; }; // namespace flightaware::uat }; // namespace flightaware diff --git a/dump978_main.cc b/dump978_main.cc index 5e08dfd..b9e5db8 100644 --- a/dump978_main.cc +++ b/dump978_main.cc @@ -12,6 +12,7 @@ #include #include +#include "adaptive.h" #include "convert.h" #include "demodulator.h" #include "exception.h" @@ -98,6 +99,23 @@ static int realmain(int argc, char **argv) { ("sdr-antenna", po::value(), "set SDR antenna name") ("sdr-stream-settings", po::value(), "set SDR stream key-value settings") ("sdr-device-settings", po::value(), "set SDR device key-value settings") + ("adaptive-range", "enable adaptive gain management targeting a dynamic range") + ("adaptive-burst", "enable adaptive gain management based on detecting loud undecoded bursts") + ("adaptive-min-gain", po::value(), "adaptive gain: never use gains below this, in dB") + ("adaptive-max-gain", po::value(), "adaptive gain: never use gains above this, in dB") + ("adaptive-duty-cycle", po::value(), "adaptive gain: fraction of samples to inspect, in percent (default: 50)") + ("adaptive-range-target", po::value(), "adaptive gain: target dynamic range, in dB (default: auto)") + ("adaptive-range-alpha", po::value(), "adaptive gain: noise floor smoothing factor, 0..1") + ("adaptive-range-percentile", po::value(), "adaptive gain: percentile of sample magnitude to use as the noise floor") + ("adaptive-range-change-delay", po::value(), "adaptive gain: seconds to wait after a gain change before acting on noise measurements") + ("adaptive-range-scan-delay", po::value(), "adaptive gain: seconds to wait before rescanning after an abortive dynamic range scan") + ("adaptive-range-rescan-delay", po::value(), "adaptive gain: seconds to wait between periodic scans for increased gain") + ("adaptive-burst-alpha", po::value(), "adaptive gain: burst rate smoothing factor, 0..1") + ("adaptive-burst-change-delay", po::value(), "adaptive gain: seconds to wait after a gain change before acting on burst measurements") + ("adaptive-burst-loud-rate", po::value(), "adaptive gain: decrease gain if the loud undecoded burst rate exceeds this rate (per second)") + ("adaptive-burst-loud-runlength", po::value(), "adaptive gain: decrease gain if the loud undecoded burst rate is exceeded for this many consecutive seconds") + ("adaptive-burst-quiet-rate", po::value(), "adaptive gain: allow gain increases only if the loud decoded message rate is below this rate (per second)") + ("adaptive-burst-quiet-runlength", po::value(), "adaptive gain: allow gain increases only after the loud decoded message rate has been low for this many consecutive seconds") ("stratuxv3", po::value(), "read messages from Stratux v3 UAT dongle on given serial port") ("raw-port", po::value>(), "listen for connections on [host:]port and provide raw messages") ("raw-legacy-port", po::value>(), "listen for connections on [host:]port and provide raw messages, with no initial metadata header") @@ -128,6 +146,7 @@ static int realmain(int argc, char **argv) { MessageDispatch dispatch; SampleSource::Pointer sample_source; + std::shared_ptr soapy_source; MessageSource::Pointer message_source; tcp::resolver resolver(io_context); @@ -137,6 +156,63 @@ static int realmain(int argc, char **argv) { return EXIT_NO_RESTART; } + // clang-format off + static const char *adaptive_option_names[] = { + "adaptive-range", "adaptive-burst", "adaptive-min-gain", "adaptive-max-gain", "adaptive-duty-cycle", + "adaptive-range-target", "adaptive-range-alpha", "adaptive-range-percentile", + "adaptive-range-change-delay", "adaptive-range-scan-delay", "adaptive-range-rescan-delay", + "adaptive-burst-alpha", "adaptive-burst-change-delay", "adaptive-burst-loud-rate", + "adaptive-burst-loud-runlength", "adaptive-burst-quiet-rate", "adaptive-burst-quiet-runlength" + }; + // clang-format on + + for (auto option : adaptive_option_names) { + if (opts.count(option) && !opts.count("sdr")) { + std::cerr << "--" << option << " can only be used with --sdr" << std::endl; + return EXIT_NO_RESTART; + } + } + + AdaptiveParams adaptive_params; + adaptive_params.range_control = (opts.count("adaptive-range") > 0); + adaptive_params.burst_control = (opts.count("adaptive-burst") > 0); + + if ((adaptive_params.range_control || adaptive_params.burst_control) && opts.count("sdr-auto-gain")) { + std::cerr << "--sdr-auto-gain cannot be combined with adaptive gain control" << std::endl; + return EXIT_NO_RESTART; + } + + if (opts.count("adaptive-min-gain")) + adaptive_params.min_gain_db = opts["adaptive-min-gain"].as(); + if (opts.count("adaptive-max-gain")) + adaptive_params.max_gain_db = opts["adaptive-max-gain"].as(); + if (opts.count("adaptive-duty-cycle")) + adaptive_params.duty_cycle = opts["adaptive-duty-cycle"].as() / 100.0; + if (opts.count("adaptive-range-target")) + adaptive_params.range_target = opts["adaptive-range-target"].as(); + if (opts.count("adaptive-range-alpha")) + adaptive_params.range_alpha = opts["adaptive-range-alpha"].as(); + if (opts.count("adaptive-range-percentile")) + adaptive_params.range_percentile = opts["adaptive-range-percentile"].as(); + if (opts.count("adaptive-range-change-delay")) + adaptive_params.range_change_delay = opts["adaptive-range-change-delay"].as(); + if (opts.count("adaptive-range-scan-delay")) + adaptive_params.range_scan_delay = opts["adaptive-range-scan-delay"].as(); + if (opts.count("adaptive-range-rescan-delay")) + adaptive_params.range_rescan_delay = opts["adaptive-range-rescan-delay"].as(); + if (opts.count("adaptive-burst-alpha")) + adaptive_params.burst_alpha = opts["adaptive-burst-alpha"].as(); + if (opts.count("adaptive-burst-change-delay")) + adaptive_params.burst_change_delay = opts["adaptive-burst-change-delay"].as(); + if (opts.count("adaptive-burst-loud-rate")) + adaptive_params.burst_loud_rate = opts["adaptive-burst-loud-rate"].as(); + if (opts.count("adaptive-burst-loud-runlength")) + adaptive_params.burst_loud_runlength = opts["adaptive-burst-loud-runlength"].as(); + if (opts.count("adaptive-burst-quiet-rate")) + adaptive_params.burst_quiet_rate = opts["adaptive-burst-quiet-rate"].as(); + if (opts.count("adaptive-burst-quiet-runlength")) + adaptive_params.burst_quiet_runlength = opts["adaptive-burst-quiet-runlength"].as(); + if (opts.count("stdin")) { sample_source = StdinSampleSource::Create(io_context, opts); } else if (opts.count("file")) { @@ -144,7 +220,8 @@ static int realmain(int argc, char **argv) { sample_source = FileSampleSource::Create(io_context, path, opts); } else if (opts.count("sdr")) { auto device = opts["sdr"].as(); - sample_source = SoapySampleSource::Create(io_context, device, opts); + soapy_source = SoapySampleSource::Create(io_context, device, opts); + sample_source = soapy_source; } else if (opts.count("stratuxv3")) { auto path = opts["stratuxv3"].as(); message_source = StratuxSerial::Create(io_context, path); @@ -239,7 +316,15 @@ static int realmain(int argc, char **argv) { sample_source->Init(); auto format = sample_source->Format(); - auto receiver = std::make_shared(format); + AdaptiveGainControl::Pointer adaptive; + if (soapy_source && (adaptive_params.range_control || adaptive_params.burst_control)) { + adaptive = AdaptiveGainControl::Create(adaptive_params, soapy_source, 2083333); + if (!adaptive->Enabled()) { + adaptive.reset(); + } + } + + auto receiver = std::make_shared(format, adaptive); sample_source->SetConsumer([receiver](std::uint64_t timestamp, const Bytes &buffer) { receiver->HandleSamples(timestamp, buffer.begin(), buffer.end()); }); diff --git a/soapy_source.cc b/soapy_source.cc index 7b5221e..622ac1d 100644 --- a/soapy_source.cc +++ b/soapy_source.cc @@ -150,6 +150,20 @@ SoapySampleSource::SoapySampleSource(boost::asio::io_context &service, const std SoapySampleSource::~SoapySampleSource() { Stop(); } +double SoapySampleSource::GetGainDb() { return device_->getGain(SOAPY_SDR_RX, 0); } + +double SoapySampleSource::SetGainDb(double db) { + device_->setGain(SOAPY_SDR_RX, 0, db); + return device_->getGain(SOAPY_SDR_RX, 0); +} + +void SoapySampleSource::GetGainLimits(double &min, double &max, double &step) { + auto range = device_->getGainRange(SOAPY_SDR_RX, 0); + min = range.minimum(); + max = range.maximum(); + step = range.step(); +} + void SoapySampleSource::Init() { try { void (*unmake)(SoapySDR::Device *) = &SoapySDR::Device::unmake; // select the right overload diff --git a/soapy_source.h b/soapy_source.h index ee0278f..6e36049 100644 --- a/soapy_source.h +++ b/soapy_source.h @@ -13,13 +13,14 @@ #include +#include "adaptive.h" #include "sample_source.h" namespace flightaware { namespace uat { - class SoapySampleSource : public SampleSource { + class SoapySampleSource : public SampleSource, public SdrGainControl { public: - static SampleSource::Pointer Create(boost::asio::io_context &service, const std::string &device_name, const boost::program_options::variables_map &options) { return Pointer(new SoapySampleSource(service, device_name, options)); } + static std::shared_ptr Create(boost::asio::io_context &service, const std::string &device_name, const boost::program_options::variables_map &options) { return std::shared_ptr(new SoapySampleSource(service, device_name, options)); } virtual ~SoapySampleSource(); @@ -28,6 +29,11 @@ namespace flightaware { void Stop() override; SampleFormat Format() override { return format_; } + // SdrGainControl; valid only after Init() + double GetGainDb() override; + double SetGainDb(double db) override; + void GetGainLimits(double &min, double &max, double &step) override; + private: SoapySampleSource(boost::asio::io_context &service, const std::string &device_name, const boost::program_options::variables_map &options);