dustcouncil/bandsaunter
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
bandsaunter/bandsaunter/adsb.py
The Dust Council f9b21f7e35 Give the aerodromes a colour of their own, and say what each aircraft is 
Four things about the animated pictures, and the window kept in step with
them.

The aerodromes were amber, and so is an aeroplane at twelve thousand feet.
Sixteen units of CIELAB apart is not two colours, it is one: an aircraft low
over a field was drawn in the field's own colour and neither could be picked
out from the other.  They are magenta now, fifty-six units from the nearest
altitude colour, which is what the ramp leaves free once red, amber, green,
cyan and violet have gone on height -- and what an aeronautical chart marks
an aerodrome in anyway.  The test states that as the distance rather than as
the colour, so that changing the ramp cannot quietly walk an aircraft back
into the airports.

The height beside an aircraft was the flight level, which is shorter and is
what an aviator reads, but "376" is only a height to somebody who already
knows it is one.  It is feet with the unit on it now, rounded to the
twenty-five feet Mode S reports altitude in: a real reading is a multiple of
that and comes through untouched, while a moment interpolated between two
reports stops claiming to know the height to the foot.

The flag of the country of registration now comes off the address block
where no register answered.  Taking it from the register's answer alone left
the flag off exactly the aircraft that had nothing else beside them either.
Mexico was missing from the address table while we were in there, which a
receiver in the American southwest notices; two registers independently give
XA- registrations for that block.

And what sort of aircraft it is, which is two facts and not one.  What it is
comes off the air: every identification message carries three bits under its
type code saying whether it is light, large, heavy, a rotorcraft, a glider,
a drone or a van on the apron, and that is the only word about what an
aircraft *is* that needs no register.  They were being decoded and thrown
away.  They are inside the identification frame the log already writes down
in full, so every log this program has ever written has them, including the
ones written before anything here knew to look.

Whether it is military comes off no air at all -- a tanker calls itself
heavy exactly as an airliner does -- and is read from the address block
instead.  On one evening here AE07D3 broadcast "heavy" and sat in the United
States military block; the register, asked separately, came back with a
C-17A Globemaster III, tail 90-0534, United States Air Force.

Labels in an animation now behave as the boxes in the window do.  One keeps
its place for as long as that place still works and is moved only when
something takes it, which on a real log is about half as many moves as
deciding afresh every frame; the moves that are left are eased over half a
second of playback; and what the next label is laid out against is where a
moving one is going rather than where it has reached.  A still has no frame
before it and places its labels exactly as it always did.

The fading was only half done.  A label's name faded, because it is drawn in
the aircraft's own colour, but its rows and its flag were fixed colours and
stayed at full brightness -- so the brightest thing on that part of the
picture was the one aeroplane nothing had been heard from.  An indexed
picture cannot blend, so the row grey and all twelve flag colours now have
dimmed copies at each of the four fade steps, and the whole box goes
together.  A crowded frame, which drops back to the short label, drops the
class and the flag with it.

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

829 lines
33 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", "category_name",
"EMITTER_CATEGORIES"]
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)
# A position takes an even frame and an odd one, and the pair is only good
# for as long as the aircraft has not meaningfully moved between them. Ten
# seconds is what the standard allows; past that the two halves describe
# different places and the answer is not a position at all.
CPR_PAIR_SECONDS = 10.0
# How long a previous position stays worth checking a new one against.
CPR_TRUST_SECONDS = 300.0
# Faster than anything with a transponder on it, so that a real aircraft is
# never called an error -- Concorde cruised at 1150 kt.
MAX_GROUND_SPEED_KT = 2000.0
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)",
}
# The three bits under an identification message's type code say what sort
# of thing is transmitting. Which list they index depends on the type code
# -- the same three bits mean "heavy" under type 4 and "glider" under type 3
# -- which is why this is a table of tables and not a table.
#
# This is the only word about what an aircraft *is* that comes off the air.
# Everything else -- the model, the operator, the registration -- is a
# lookup in somebody's database against the address. There is no category
# for "military": that has to be read off the address block instead.
EMITTER_CATEGORIES = {
4: {1: "light", 2: "small", 3: "large", 4: "high vortex",
5: "heavy", 6: "high performance", 7: "rotorcraft"},
3: {1: "glider", 2: "airship", 3: "parachutist", 4: "ultralight",
6: "drone", 7: "spacecraft"},
2: {1: "emergency vehicle", 2: "service vehicle", 3: "obstacle",
4: "obstacle", 5: "obstacle"},
1: {}, # reserved, and nothing transmits it
}
def category_name(type_code: int, category: int) -> str:
"""What an identification message says the transmitter is, or "".
Zero means the aircraft declined to say, which is common and is not an
error: it is answered with nothing rather than with a guess.
"""
return EMITTER_CATEGORIES.get(type_code, {}).get(category, "")
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
category: int = 0 # the emitter category, with the type code
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)
# The three bits under the type code say what sort of thing this is
# -- heavy, rotorcraft, glider, ground vehicle. It comes off the
# air with the callsign and costs nothing to keep.
frame.category = int(me[5:8], 2)
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
category: str = "" # what it said it was: heavy, rotorcraft...
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)
_placed_at: float = field(default=0.0, 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.category and not seen.category:
seen.category = category_name(frame.type_code, frame.category)
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
self._place(seen)
return seen
def _place(self, seen: Aircraft) -> None:
"""Work out where an aircraft is, from the last even and odd frames.
The pair has to be recent, and the answer has to be reachable. Both
checks are the difference between a map and a scatter of nonsense:
an unpaired frame kept from ten minutes ago decodes against a fresh
one to a position on the wrong side of the world, because compact
position reporting sends a fraction of a zone and the two fractions
are then read as though the aircraft had not moved between them.
Measured against one night's recording, that produced positions up
to seven thousand miles out, on two aircraft in three.
"""
even, odd = seen._even, seen._odd
if even is None or odd is None:
return
if abs(even.received_at - odd.received_at) > CPR_PAIR_SECONDS:
return
# 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 = (even.received_at, even.at_sample) > \
(odd.received_at, odd.at_sample)
found = global_position(even, odd, even_first)
if found is None or not _on_earth(*found):
return
when = max(even.received_at, odd.received_at)
if not seen.located or self._reachable(seen, found, when):
seen.latitude, seen.longitude = found
seen._placed_at = when
@staticmethod
def _reachable(seen: Aircraft, found: tuple[float, float],
when: float) -> bool:
"""Whether an aircraft could have got there from where it was.
A position that would need eight hundred knots in the second since
the last one is not a position, whatever the checksum said about the
frames it came from. The bar is set well above anything that flies
so that a genuinely fast aircraft, or a gap in reception, is never
mistaken for an error.
"""
gap = when - seen._placed_at
if gap <= 0 or gap > CPR_TRUST_SECONDS:
return True # too long ago to argue with
miles = _distance_nm(seen.latitude, seen.longitude, *found)
return miles <= MAX_GROUND_SPEED_KT * gap / 3600.0
def _on_earth(lat: float, lon: float) -> bool:
"""Whether a decoded position is a place at all."""
return -90.0 <= lat <= 90.0 and -180.0 <= lon <= 180.0
def _distance_nm(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
"""Great-circle distance, in nautical miles."""
p1, p2 = math.radians(lat1), math.radians(lat2)
dp = p2 - p1
dl = math.radians(lon2 - lon1)
a = math.sin(dp / 2) ** 2 + \
math.cos(p1) * math.cos(p2) * math.sin(dl / 2) ** 2
return 2 * 3440.065 * math.asin(min(1.0, math.sqrt(a)))
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