dustcouncil/bandsaunter
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
bandsaunter/bandsaunter/adsb.py
The Dust Council eae60cb04d Write the aircraft down, and draw where they went 
ADS-B was a live table and nothing else: an aircraft was overhead for four
minutes and then gone, with nothing kept.  Now everything heard goes into
adsb_<time>.jsonl as it arrives -- one object per frame, the raw hex beside
what was read out of it, flushed per line because a listening session ends
with control-C -- with a readable report beside it.

flights.py asks who the aircraft are: adsbdb for the airframe and the
route, hexdb behind it, cached for a month.  What needs no website is
answered without one, because the ICAO address block says which country
registered the aircraft and the first three letters of an airline callsign
are its designator.  Nothing but the address and the callsign heard on the
air is ever sent.

  bandsaunter flights [LOG...] --out sky.gif

reads a log back and draws the evening as a map with the clock running.
Every frame is a moment: each aircraft is where it actually was then,
interpolated between the position reports either side of it and
dead-reckoned from its last speed and heading between them, and dropped
rather than guessed at once it has not been heard for --stale seconds.
The GIF is written here -- palette, LZW, frame differencing against a
transparent index -- so nothing but numpy is needed; ffmpeg writes an MP4
where it happens to be installed, and .png draws the whole evening at once.

The decoder needed 6.3 s to read a second of sky, so a live capture was
losing six frames in seven.  Reading the bits off a running total instead
of summing each window takes that to 0.6 s, with identical output.

--simulate flies six aircraft that are not there past a receiver that is
not there, through the real encoder, the real checksum and the real
decoder, so all of this can be tried without an aerial.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_016PsWPTweCT6pwxKngvVxcg
2026-09-03 22:22:49 -07:00

723 lines
28 KiB
Python
Raw Blame History

"""ADS-B: aircraft on 1090 MHz saying where they are.
Every airliner overhead broadcasts its identity, position, altitude and speed
twice a second, unencrypted, to nobody in particular. The format is Mode S
extended squitter and it is the easiest useful thing a receiver can decode,
because every frame carries a 24-bit checksum that either comes out right or
does not -- there is no judgement anywhere in this module about whether a
decode is believable.
The signalling is pulse-position modulation at a megabit a second. Each bit
is one microsecond wide and split in half: energy in the first half is a one,
energy in the second half is a zero. A frame opens with a preamble of four
pulses at 0, 1, 3.5 and 4.5 microseconds, which is what is searched for.
That rate is why this does not ride on the ordinary scan path. A megabit a
second needs at least two megasamples a second of raw receiver output, and
the scanner's channels are twelve and a half kilohertz wide; ``bandsaunter
adsb`` parks the receiver on 1090 MHz at full rate instead.
"""
from __future__ import annotations
import math
import random
import time
from dataclasses import dataclass, field
import numpy as np
__all__ = ["decode_adsb", "Frame", "Aircraft", "AircraftRegistry", "crc24",
"ADSB_HZ", "SAMPLE_RATE", "PREAMBLE_US", "encode_identification",
"encode_position", "encode_velocity", "modulate", "SimulatedSky",
"VirtualAircraft", "default_sky"]
ADSB_HZ = 1_090_000_000.0
# The lowest rate this can work at: one megabit a second, sampled twice a bit.
SAMPLE_RATE = 2_000_000
PREAMBLE_US = (0.0, 1.0, 3.5, 4.5)
SHORT_BITS = 56
LONG_BITS = 112
# The characters a callsign can be built from, six bits each. The hashes are
# the code points the standard leaves unassigned.
CALLSIGN_CHARS = ("#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############"
"0123456789######")
TYPE_NAMES = {
(1, 4): "identification",
(5, 8): "surface position",
(9, 18): "airborne position",
(19, 19): "velocity",
(20, 22): "airborne position (GNSS height)",
}
def crc24(data: bytes) -> int:
"""The Mode S parity, polynomial 0xFFF409.
Run over the whole frame including its three parity bytes, a good frame
gives zero. That is the entire error check in this module and it is
enough: twenty-four bits of it means a frame passes by chance one time in
sixteen million.
"""
poly = 0xFFF409
crc = 0
for byte in data:
crc ^= byte << 16
for _ in range(8):
crc = (((crc << 1) ^ poly) & 0xFFFFFF if crc & 0x800000
else (crc << 1) & 0xFFFFFF)
return crc
# ---------------------------------------------------------------------------
# One frame
# ---------------------------------------------------------------------------
@dataclass
class Frame:
"""One Mode S frame that passed its checksum."""
bits: str = ""
data: bytes = b""
df: int = 0 # downlink format
icao: str = "" # the aircraft's permanent 24-bit address
type_code: int = 0
at_sample: int = 0
callsign: str = ""
altitude_ft: int = 0
latitude: float = 0.0
longitude: float = 0.0
cpr_odd: bool = False
cpr_lat: int = 0
cpr_lon: int = 0
received_at: float = 0.0 # when it arrived, in whatever clock the caller keeps
ground_speed_kt: float = 0.0
track_deg: float = 0.0
vertical_rate_fpm: int = 0
@property
def what(self) -> str:
for (low, high), name in TYPE_NAMES.items():
if low <= self.type_code <= high:
return name
return f"type {self.type_code}"
def describe(self) -> str:
bits = [self.icao]
if self.callsign:
bits.append(self.callsign)
if self.altitude_ft:
bits.append(f"{self.altitude_ft} ft")
if self.latitude or self.longitude:
bits.append(f"{self.latitude:.4f},{self.longitude:.4f}")
if self.ground_speed_kt:
bits.append(f"{self.ground_speed_kt:.0f} kt")
bits.append(f"{self.track_deg:.0f}°")
if self.vertical_rate_fpm:
bits.append(f"{self.vertical_rate_fpm:+d} fpm")
if len(bits) == 1:
bits.append(self.what)
return " ".join(bits)
# ---------------------------------------------------------------------------
# Finding frames in the samples
# ---------------------------------------------------------------------------
def _magnitude(iq: np.ndarray) -> np.ndarray:
if np.iscomplexobj(iq):
return np.abs(iq).astype(np.float32)
return np.abs(np.asarray(iq, dtype=np.float32))
def _preamble_score(mag: np.ndarray, per_us: float) -> np.ndarray:
"""How much each sample looks like the start of a preamble.
The four pulses minus the four gaps that have to be quiet between them.
A correlation with the pattern alone finds a steady carrier just as
happily; requiring the gaps is what makes it a preamble.
"""
pulses = [int(round(us * per_us)) for us in PREAMBLE_US]
quiet = [int(round(us * per_us)) for us in (2.0, 2.5, 3.0, 6.0, 6.5, 7.0)]
width = max(1, int(round(0.5 * per_us)))
need = int(round(8.0 * per_us))
if mag.size < need + width:
return np.zeros(0, dtype=np.float32)
def at(offsets):
total = np.zeros(mag.size - need - width, dtype=np.float32)
for offset in offsets:
total += mag[offset:offset + total.size]
return total / len(offsets)
return at(pulses) - at(quiet)
def _bits_at(running: np.ndarray, start: int, per_us: float,
count: int) -> str:
"""Read ``count`` pulse-position bits: loud first half is a one.
``running`` is a running total of the magnitudes, so the energy in any
stretch of samples is one subtraction rather than a sum. It matters:
every candidate preamble in a second of receiver output is tried at two
lengths, which is a quarter of a million half-microsecond windows a
second, and adding them up one at a time is the difference between
keeping up with the sky and hearing one frame in seven.
"""
half = 0.5 * per_us
firsts = start + 8.0 * per_us + np.arange(count) * per_us
a0 = np.rint(firsts).astype(np.int64)
a1 = np.rint(firsts + half).astype(np.int64)
b1 = np.rint(firsts + per_us).astype(np.int64)
size = running.size - 1
if b1[-1] > size:
keep = int(np.searchsorted(b1, size, side="right"))
a0, a1, b1 = a0[:keep], a1[:keep], b1[:keep]
early = running[a1] - running[a0]
late = running[b1] - running[a1]
return "".join(np.where(early > late, "1", "0"))
# The formats whose parity is the checksum itself. Everything else has the
# aircraft's address exclusive-ored into it, so a decoder without a list of
# the aircraft it expects to hear cannot check them at all.
CHECKABLE_FORMATS = frozenset({11, 17, 18})
def _plausible(data: bytes) -> bool:
"""Whether a candidate is a real frame rather than a run of silence.
The checksum does the work, with one exception that matters: a frame of
all zeros passes it, because zero divided by anything leaves nothing.
Silence between transmissions is exactly that, so it would otherwise
decode as an endless stream of aircraft 000000.
"""
if not any(data):
return False
if (data[0] >> 3) not in CHECKABLE_FORMATS:
return False
return crc24(data) == 0
def _bytes_of(bits: str) -> bytes:
whole = len(bits) - len(bits) % 8
return bytes(int(bits[i:i + 8], 2) for i in range(0, whole, 8))
def decode_frames(iq: np.ndarray, sample_rate: float) -> list[Frame]:
"""Every Mode S frame in a block of raw receiver output.
Only frames whose checksum comes out right are returned, so there is no
threshold to tune and nothing to disbelieve.
"""
if sample_rate < SAMPLE_RATE * 0.99:
return []
mag = _magnitude(iq)
per_us = sample_rate / 1e6
score = _preamble_score(mag, per_us)
if score.size == 0:
return []
# A preamble stands well above the noise around it. This is only a
# shortlist -- the checksum decides -- so it is set low enough to let
# weak frames through and high enough not to try every sample.
floor = float(np.median(mag)) * 2.0
candidates = np.flatnonzero(score > max(floor, float(np.std(score))))
# Once, for the whole block: every window a candidate asks about is then
# the difference between two of these.
running = np.concatenate(([0.0], np.cumsum(mag, dtype=np.float64)))
frames: list[Frame] = []
# The candidates come out in order, so the only accepted frame a new one
# can overlap is the last of them.
taken = -per_us * 2
for start in candidates:
if start - taken < per_us:
continue
for count in (LONG_BITS, SHORT_BITS):
bits = _bits_at(running, int(start), per_us, count)
if len(bits) < count:
continue
data = _bytes_of(bits)
if len(data) * 8 != count or not _plausible(data):
continue
frame = _read(bits, data)
frame.at_sample = int(start)
frames.append(frame)
taken = int(start)
break
return frames
# ---------------------------------------------------------------------------
# What a frame says
# ---------------------------------------------------------------------------
def _read(bits: str, data: bytes) -> Frame:
frame = Frame(bits=bits, data=data, df=data[0] >> 3)
frame.icao = f"{int.from_bytes(data[1:4], 'big'):06X}"
if frame.df not in (17, 18) or len(data) < 11:
return frame
me = bits[32:88]
frame.type_code = int(me[:5], 2)
if 1 <= frame.type_code <= 4:
frame.callsign = _callsign(me)
elif 9 <= frame.type_code <= 18 or 20 <= frame.type_code <= 22:
frame.altitude_ft = _altitude(me)
frame.cpr_odd = me[21] == "1"
frame.cpr_lat = int(me[22:39], 2)
frame.cpr_lon = int(me[39:56], 2)
elif frame.type_code == 19:
_velocity(frame, me)
return frame
def _callsign(me: str) -> str:
chars = [CALLSIGN_CHARS[int(me[8 + 6 * i:14 + 6 * i], 2)] for i in range(8)]
return "".join(chars).replace("#", "").strip("_ ").strip()
def _altitude(me: str) -> int:
"""The 12-bit altitude field, in feet.
The Q bit says which of two encodings is in use: 25-foot steps, which is
everything in normal service, or the older 100-foot Gillham code, which
is not decoded here and comes back as zero rather than as a wrong number.
"""
field = me[8:20]
if field == "0" * 12:
return 0
q_bit = field[7]
if q_bit != "1":
return 0
value = int(field[:7] + field[8:], 2)
return value * 25 - 1000
def _velocity(frame: Frame, me: str) -> None:
subtype = int(me[5:8], 2)
if subtype not in (1, 2):
return # airspeed rather than ground speed
sign_ew = -1 if me[13] == "1" else 1
ew = int(me[14:24], 2) - 1
sign_ns = -1 if me[24] == "1" else 1
ns = int(me[25:35], 2) - 1
if ew < 0 or ns < 0:
return
scale = 4.0 if subtype == 2 else 1.0 # supersonic
vx, vy = sign_ew * ew * scale, sign_ns * ns * scale
frame.ground_speed_kt = math.hypot(vx, vy)
frame.track_deg = math.degrees(math.atan2(vx, vy)) % 360.0
rate = int(me[37:46], 2)
if rate:
frame.vertical_rate_fpm = (rate - 1) * 64 * (-1 if me[36] == "1" else 1)
# ---------------------------------------------------------------------------
# Position, which takes two frames
# ---------------------------------------------------------------------------
def _nl(lat: float) -> int:
"""How many longitude zones there are at this latitude."""
if abs(lat) >= 87.0:
return 1
if lat == 0:
return 59
inner = 1 - (1 - math.cos(math.pi / (2 * 15))) / \
math.cos(math.radians(abs(lat))) ** 2
inner = max(-1.0, min(1.0, inner))
return int(math.floor(2 * math.pi / math.acos(inner)))
def global_position(even: Frame, odd: Frame,
even_first: bool = True) -> tuple[float, float] | None:
"""Where an aircraft is, from one even and one odd position frame.
Compact position reporting sends a latitude and longitude as fractions of
a zone, with the zones laid out differently in the two frame types. One
frame alone is ambiguous by hundreds of miles; the pair is not.
"""
lat_even, lon_even = even.cpr_lat / 131072.0, even.cpr_lon / 131072.0
lat_odd, lon_odd = odd.cpr_lat / 131072.0, odd.cpr_lon / 131072.0
j = math.floor(59 * lat_even - 60 * lat_odd + 0.5)
rlat_even = (360.0 / 60) * ((j % 60) + lat_even)
rlat_odd = (360.0 / 59) * ((j % 59) + lat_odd)
if rlat_even >= 270:
rlat_even -= 360
if rlat_odd >= 270:
rlat_odd -= 360
if _nl(rlat_even) != _nl(rlat_odd):
return None # the two straddle a zone boundary
nl = _nl(rlat_even)
if even_first:
lat = rlat_even
ni = max(nl, 1)
m = math.floor(lon_even * (nl - 1) - lon_odd * nl + 0.5)
lon = (360.0 / ni) * ((m % ni) + lon_even)
else:
lat = rlat_odd
ni = max(nl - 1, 1)
m = math.floor(lon_even * (nl - 1) - lon_odd * nl + 0.5)
lon = (360.0 / ni) * ((m % ni) + lon_odd)
if lon >= 180:
lon -= 360
return lat, lon
@dataclass
class Aircraft:
"""What has been heard from one aircraft."""
icao: str
callsign: str = ""
altitude_ft: int = 0
latitude: float = 0.0
longitude: float = 0.0
ground_speed_kt: float = 0.0
track_deg: float = 0.0
vertical_rate_fpm: int = 0
messages: int = 0
first_seen: float = 0.0
last_seen: float = 0.0
_even: Frame | None = field(default=None, repr=False)
_odd: Frame | None = field(default=None, repr=False)
@property
def located(self) -> bool:
return bool(self.latitude or self.longitude)
def describe(self) -> str:
bits = [self.icao]
if self.callsign:
bits.append(self.callsign)
if self.located:
bits.append(f"{self.latitude:.4f},{self.longitude:.4f}")
if self.altitude_ft:
bits.append(f"{self.altitude_ft} ft")
if self.ground_speed_kt:
bits.append(f"{self.ground_speed_kt:.0f} kt {self.track_deg:.0f}°")
return " ".join(bits)
class AircraftRegistry:
"""Everything heard, gathered by aircraft rather than by frame.
Position needs an even frame and an odd one, which arrive half a second
apart, so something has to remember the first while the second is on its
way. This is that, and it is also what turns four hundred frames into
the dozen aircraft they came from.
"""
def __init__(self):
self.aircraft: dict[str, Aircraft] = {}
def __len__(self) -> int:
return len(self.aircraft)
def add(self, frame: Frame, when: float = 0.0) -> Aircraft:
seen = self.aircraft.get(frame.icao)
if seen is None:
seen = Aircraft(icao=frame.icao, first_seen=when)
self.aircraft[frame.icao] = seen
seen.messages += 1
seen.last_seen = when
frame.received_at = when
if frame.callsign:
seen.callsign = frame.callsign
if frame.altitude_ft:
seen.altitude_ft = frame.altitude_ft
if frame.ground_speed_kt:
seen.ground_speed_kt = frame.ground_speed_kt
seen.track_deg = frame.track_deg
if frame.vertical_rate_fpm:
seen.vertical_rate_fpm = frame.vertical_rate_fpm
if frame.cpr_lat or frame.cpr_lon:
if frame.cpr_odd:
seen._odd = frame
else:
seen._even = frame
if seen._even is not None and seen._odd is not None:
# Whichever of the pair arrived later is the one the position
# is reported at. Compared by arrival rather than by sample
# offset: the offset restarts at zero every block, so a pair
# that straddles two blocks would otherwise be read backwards
# and put the aircraft in the wrong zone.
even_first = (seen._even.received_at, seen._even.at_sample) > \
(seen._odd.received_at, seen._odd.at_sample)
found = global_position(seen._even, seen._odd, even_first)
if found is not None:
seen.latitude, seen.longitude = found
return seen
def described(self) -> list[str]:
return [craft.describe() for craft in
sorted(self.aircraft.values(), key=lambda a: a.icao)]
def decode_adsb(iq: np.ndarray, sample_rate: float,
registry: AircraftRegistry | None = None
) -> tuple[list[Frame], AircraftRegistry]:
"""Decode every frame in a block, and fold them into aircraft."""
registry = registry if registry is not None else AircraftRegistry()
frames = decode_frames(iq, sample_rate)
for frame in frames:
registry.add(frame, when=frame.at_sample / sample_rate)
return frames, registry
# ---------------------------------------------------------------------------
# The other direction: making frames, for a receiver that has no aerial
# ---------------------------------------------------------------------------
#
# ADS-B is the one thing in this program that cannot be tried out indoors. A
# scanner can be pointed at a simulated transmitter, but 1090 MHz needs an
# aerial cut for it and an aeroplane in the sky, and a person deciding whether
# any of this is worth wiring up has neither. So the frames can be built as
# well as read, and a sky full of imaginary aircraft can be flown past an
# imaginary receiver: the same encoding, the same checksum, the same decoder.
def _with_parity(payload: bytes) -> bytes:
"""A frame with its 24 parity bits on the end, as a transmitter sends it."""
return payload + crc24(payload).to_bytes(3, "big")
def _squitter(icao: int, me: bytes) -> bytes:
"""DF17, capability 5, one aircraft address and 56 bits of message."""
return _with_parity(bytes([17 << 3 | 5]) + (icao & 0xFFFFFF).to_bytes(3, "big")
+ bytes(me))
def encode_identification(icao: int, callsign: str, category: int = 0) -> bytes:
"""The frame an aircraft sends to say what it is called."""
text = callsign.upper().ljust(8)[:8]
bits = ""
for char in text:
index = CALLSIGN_CHARS.find(char)
bits += format(index if index >= 0 else 32, "06b")
me = bytes([(4 << 3) | (category & 0x07)]) + int(bits, 2).to_bytes(6, "big")
return _squitter(icao, me)
def _cpr_encode(lat: float, lon: float, odd: bool) -> tuple[int, int]:
"""Compact position reporting, the transmitting side of :func:`global_position`."""
i = 1 if odd else 0
d_lat = 360.0 / (60 - i)
y = int(round(131072 * ((lat % d_lat) / d_lat)))
zones = _nl(lat) - i
d_lon = 360.0 / zones if zones > 0 else 360.0
x = int(round(131072 * ((lon % d_lon) / d_lon)))
return y & 0x1FFFF, x & 0x1FFFF
def encode_position(icao: int, lat: float, lon: float, altitude_ft: int,
odd: bool) -> bytes:
"""An airborne position frame: where the aircraft is and how high.
Half a position, strictly: it takes an even frame and an odd one to say
where anything is, which is the whole point of the encoding.
"""
steps = max(0, int(round((altitude_ft + 1000) / 25.0)))
field_bits = format(min(steps, 0x7FF), "011b")
altitude = field_bits[:7] + "1" + field_bits[7:] # the Q bit: 25 ft
y, x = _cpr_encode(lat, lon, odd)
me_bits = (format(11, "05b") + "000" + altitude + "0"
+ ("1" if odd else "0")
+ format(y, "017b") + format(x, "017b"))
return _squitter(icao, int(me_bits, 2).to_bytes(7, "big"))
def encode_velocity(icao: int, east_kt: float, north_kt: float,
vertical_fpm: int = 0) -> bytes:
"""A velocity frame: ground speed as two components, and climb rate."""
east, north = int(round(east_kt)), int(round(north_kt))
rate = min(511, abs(int(vertical_fpm)) // 64 + 1) if vertical_fpm else 0
me_bits = (format(19, "05b") + "001" + "00000"
+ ("1" if east < 0 else "0") + format(min(1023, abs(east) + 1), "010b")
+ ("1" if north < 0 else "0") + format(min(1023, abs(north) + 1), "010b")
+ "0" + ("1" if vertical_fpm < 0 else "0")
+ format(rate, "09b") + "0" * 10)
return _squitter(icao, int(me_bits[:56], 2).to_bytes(7, "big"))
def _burst(frame: bytes, per_us: float, amplitude: float = 1.0) -> np.ndarray:
"""One frame as the magnitude a receiver sees: preamble, then the bits."""
bits = "".join(format(byte, "08b") for byte in frame)
span = np.zeros(int(round((8 + len(bits) + 1) * per_us)), dtype=np.float32)
half = int(round(0.5 * per_us))
for at in PREAMBLE_US:
lo = int(round(at * per_us))
span[lo:lo + half] = amplitude
for i, bit in enumerate(bits):
base = (8 + i) * per_us
lo = int(round(base if bit == "1" else base + 0.5 * per_us))
span[lo:lo + half] = amplitude
return span
def modulate(frames, sample_rate: float = SAMPLE_RATE, gap_us: float = 60.0,
amplitude: float = 1.0, noise: float = 0.0,
seed: int = 0) -> np.ndarray:
"""Turn frames into what a receiver on 1090 MHz would have heard."""
per_us = sample_rate / 1e6
gap = np.zeros(int(round(gap_us * per_us)), dtype=np.float32)
parts = [gap]
for frame in frames:
parts.append(_burst(frame, per_us, amplitude))
parts.append(gap)
signal = np.concatenate(parts)
if noise:
rng = np.random.default_rng(seed)
signal = signal + noise * np.abs(rng.standard_normal(signal.size))
return signal.astype(np.complex64)
@dataclass
class VirtualAircraft:
"""An aeroplane that does not exist, flying in a straight line.
Enough of an aircraft to be worth drawing: it has an address, a callsign,
a place to be, a speed to get there at and a rate of climb. It reports
itself exactly as a real one does, so nothing downstream can tell the
difference -- which is the point, because everything downstream is being
tested.
"""
icao: int = 0
callsign: str = ""
latitude: float = 0.0
longitude: float = 0.0
altitude_ft: int = 30_000
speed_kt: float = 420.0
heading_deg: float = 90.0
climb_fpm: int = 0
strength: float = 1.0
def advance(self, seconds: float) -> None:
"""Fly on for a while, which is all this aircraft knows how to do."""
nm = self.speed_kt * seconds / 3600.0
theta = math.radians(self.heading_deg)
self.latitude += nm / 60.0 * math.cos(theta)
# A minute of longitude is a minute of latitude times the cosine, and
# at eighty degrees north that difference is most of the answer.
self.longitude += nm / 60.0 * math.sin(theta) / max(
0.05, math.cos(math.radians(self.latitude)))
self.altitude_ft = max(0, int(self.altitude_ft
+ self.climb_fpm * seconds / 60.0))
def frames(self, second: int) -> list[bytes]:
"""What it broadcasts in one second: position, velocity, sometimes a name."""
theta = math.radians(self.heading_deg)
out = [encode_position(self.icao, self.latitude, self.longitude,
self.altitude_ft, odd=False),
encode_position(self.icao, self.latitude, self.longitude,
self.altitude_ft, odd=True),
encode_velocity(self.icao, self.speed_kt * math.sin(theta),
self.speed_kt * math.cos(theta), self.climb_fpm)]
if second % 5 == 0 and self.callsign:
out.insert(0, encode_identification(self.icao, self.callsign))
return out
def default_sky(latitude: float = 47.55, longitude: float = -122.30,
seed: int = 7) -> list[VirtualAircraft]:
"""A handful of aircraft around a receiver, going about their business.
Airliners at height on their way past, one climbing out, one descending
towards the airport and a helicopter going nowhere in particular: enough
different heights and speeds that a map of them is worth looking at.
"""
rng = random.Random(seed)
def near(miles: float) -> tuple[float, float]:
bearing = rng.uniform(0, 360)
return (latitude + miles / 60.0 * math.cos(math.radians(bearing)),
longitude + miles / 60.0 * math.sin(math.radians(bearing))
/ max(0.05, math.cos(math.radians(latitude))))
plan = (("UAL1902", 0xA1B2C3, 36_000, 470.0, 78.0, 0, 40.0),
("ASA412", 0xA24C71, 12_500, 310.0, 155.0, -1800, 22.0),
("SWA2311", 0xA9E0F4, 4_200, 240.0, 342.0, 2200, 14.0),
("DAL88", 0xAB1D55, 39_000, 505.0, 265.0, 0, 55.0),
("N517HP", 0xA6F109, 1_200, 95.0, 20.0, 0, 6.0),
("BAW49", 0x4008F6, 33_000, 480.0, 300.0, 640, 48.0))
sky = []
for callsign, icao, altitude, speed, heading, climb, distance in plan:
lat, lon = near(distance)
sky.append(VirtualAircraft(icao=icao, callsign=callsign, latitude=lat,
longitude=lon, altitude_ft=altitude,
speed_kt=speed, heading_deg=heading,
climb_fpm=climb,
strength=rng.uniform(0.6, 1.0)))
return sky
class SimulatedSky:
"""A receiver-shaped source of aeroplanes that are not there.
It answers ``read_samples`` like the real device does and hands back the
same magnitudes a dongle would, so ``bandsaunter adsb --simulate`` runs
every line of the decoder, the log, the lookups and the map without an
aerial, an aircraft or a licence.
"""
def __init__(self, aircraft=None, sample_rate: float = SAMPLE_RATE,
noise: float = 0.02, seed: int = 0, realtime: bool = False):
self.aircraft = list(aircraft) if aircraft is not None else default_sky()
self.sample_rate = float(sample_rate)
self.noise = noise
self.rng = np.random.default_rng(seed)
self.second = 0
self.frequency = ADSB_HZ
# A block is a second of samples but takes longer than a second to
# decode, so an aircraft advanced by the length of the block falls
# behind the clock the frames are stamped with -- and a map drawn
# from the log would show a 480-knot airliner crawling. Listening
# for real, the sky moves by the time that actually passed; in a
# test, by the block, so the same seed gives the same sky twice.
self.realtime = realtime
self._last = None
# -- the shape of a device -------------------------------------------
def open(self):
return self
def close(self) -> None:
return None
def tune(self, hz: float, settle: bool = True) -> int:
self.frequency = hz
return int(hz)
def read_samples(self, count: int, flush: bool = False) -> np.ndarray:
"""One block of sky: everyone reports, everyone moves on.
The bursts are scattered through the block rather than lined up at
the front, because two aircraft transmitting at the same moment is a
thing that happens and a decoder that has never seen it is untested.
"""
seconds = count / self.sample_rate
if self.realtime:
now = time.monotonic()
if self._last is not None:
seconds = max(0.0, now - self._last)
self._last = now
block = np.zeros(count, dtype=np.float32)
per_us = self.sample_rate / 1e6
for craft in self.aircraft:
for frame in craft.frames(self.second):
burst = _burst(frame, per_us, craft.strength)
at = int(self.rng.integers(0, max(1, count - burst.size)))
block[at:at + burst.size] = np.maximum(
block[at:at + burst.size], burst[:count - at])
craft.advance(seconds)
self.second += 1
if self.noise:
block = block + self.noise * np.abs(
self.rng.standard_normal(count)).astype(np.float32)
return block.astype(np.complex64)
def read_seconds(self, seconds: float, flush: bool = False) -> np.ndarray:
return self.read_samples(int(self.sample_rate * seconds))
Reference in a new issue View git blame Copy permalink