Aloha Habr. Malia paha ua hoʻohana ka poʻe a pau i hālāwai a ʻike paha i nā ʻohana a i ʻole nā hoaaloha ma ka mokulele i ka lawelawe manuahi Flightradar24. He ala kūpono loa kēia e nānā i ke kūlana o ka mokulele i ka manawa maoli.
В Ua wehewehe ʻia ke kumu hana o ia lawelawe pūnaewele. E hele mākou i mua a noʻonoʻo i ka ʻikepili e hoʻouna ʻia a loaʻa mai ka mokulele i ke kahua hoʻokipa a hoʻololi iā mākou iho me ka hoʻohana ʻana iā Python.
История
ʻIke loa, ʻaʻole i lawe ʻia ka ʻikepili mokulele no nā mea hoʻohana e ʻike ma kā lākou kelepona. Kapa ʻia ka ʻōnaehana ʻo ADS-B (Automatic dependent surveillance—broadcast), a hoʻohana ʻia e hoʻouna ʻokoʻa i ka ʻike e pili ana i ka mokulele i ke kikowaena hoʻokele - ua hoʻouna ʻia kona ʻike, nā hoʻonohonoho, kuhikuhi, ka wikiwiki, ke kiʻekiʻe a me nā ʻikepili ʻē aʻe. Ma mua, ma mua o ka hiki ʻana mai o ia mau ʻōnaehana, hiki i ka mea hoʻouna ke ʻike i kahi kiko ma ka radar. ʻAʻole lawa kēia i ka wā i nui loa nā mokulele.
ʻO ka ʻenehana, aia ʻo ADS-B i kahi transmitter ma luna o kahi mokulele e hoʻouna pinepine i nā ʻeke o ka ʻike ma kahi kiʻekiʻe kiʻekiʻe o 1090 MHz (he mau ʻano ʻē aʻe, akā ʻaʻole mākou makemake nui iā lākou, no ka mea, ua hoʻouna ʻia nā hoʻonohonoho ma aneʻi). ʻOiaʻiʻo, ma waho aʻe o ka transmitter, aia kekahi mea hoʻokipa ma kahi o ke kahua mokulele, akā no mākou, ma ke ʻano he mea hoʻohana, hoihoi kā mākou mea hoʻokipa.
Ma ke ala, no ka hoʻohālikelike ʻana, ʻo ka ʻōnaehana like ʻole, ʻo Airnav Radarbox, i hoʻolālā ʻia no nā mea hoʻohana maʻamau, i ʻike ʻia i ka makahiki 2007, a ua kūʻai ʻia ma kahi o $900; ʻo ke kau inoa ʻana i nā lawelawe pūnaewele e kūʻai aku i $ 250 i ka makahiki.
Hiki ke heluhelu ʻia nā loiloi o kēlā mau mea nona ka Lūkini mua ma ka ʻaha kūkā . I kēia manawa ua loaʻa nui nā mea hoʻokipa RTL-SDR, hiki ke hōʻuluʻulu ʻia kahi mea like no $30; ʻoi aku ka nui o kēia ma . E neʻe kākou i ka protocol ponoʻī - e ʻike kākou pehea e hana ai.
Loaʻa nā hōʻailona
ʻO ka mea mua, pono e hoʻopaʻa ʻia ka hōʻailona. He 120 microseconds ka lōʻihi o ka hōʻailona holoʻokoʻa, no laila e ʻoluʻolu e wehe i kāna mau ʻāpana, makemake ʻia kahi mea hoʻokipa SDR me kahi alapine sampling o ka liʻiliʻi o 5 MHz.
Ma hope o ka hoʻopaʻa ʻana, loaʻa iā mākou kahi faila WAV me ka helu hoʻohālike o 5000000 samples / sec; 30 kekona o ia hoʻopaʻa ʻana "kaumaha" ma kahi o 500MB. ʻO ka hoʻolohe ʻana iā ia me kahi mea pāʻani media, ʻoiaʻiʻo, he mea ʻole - ʻaʻole i loko o ka faila ke kani, akā he hōʻailona lekiō i hoʻopaʻa ʻia pololei - ʻo ia ke ʻano o ka hana ʻana o Software Defined Radio.
E wehe mākou a hana i ka faila me ka Python. Hiki i ka poʻe makemake e hoʻokolohua iā lākou iho ke hoʻoiho i kahi hoʻopaʻa hoʻopaʻa hoʻohālike .
E hoʻoiho i ka faila a ʻike i ka mea i loko.
from scipy.io import wavfile
import matplotlib.pyplot as plt
import numpy as np
fs, data = wavfile.read("adsb_20190311_191728Z_1090000kHz_RF.wav")
data = data.astype(float)
I, Q = data[:, 0], data[:, 1]
A = np.sqrt(I*I + Q*Q)
plt.plot(A)
plt.show()
Ka hopena: ʻike mākou i nā "pulses" kūʻē i ka walaʻau hope.
ʻO kēlā me kēia "pulse" he hōʻailona, ʻike maopopoʻia keʻano o ia mea inā hoʻonuiʻoe i ka hoʻonā ma ka pakuhi.
E like me kāu e ʻike ai, ua kūlike ke kiʻi me ka mea i hāʻawi ʻia ma ka wehewehe ma luna. Hiki iā ʻoe ke hoʻomaka i ka hana ʻana i ka ʻikepili.
Wehewehewehe
ʻO ka mea mua, pono ʻoe e kiʻi i kahi kahawai liʻiliʻi. Hoʻopili ʻia ka hōʻailona ponoʻī me ka hoʻohana ʻana i Manchester encoding:
Mai ka ʻokoʻa pae i nā nibbles he maʻalahi ke kiʻi i ka "0" a me "1".
bits_str = ""
for p in range(8):
pos = start_data + bit_len*p
p1, p2 = A[pos: pos + bit_len/2], A[pos + bit_len/2: pos + bit_len]
avg1, avg2 = np.average(p1), np.average(p2)
if avg1 < avg2:
bits_str += "0"
elif avg1 > avg2:
bits_str += "1"
ʻO ke ʻano o ka hōʻailona ponoʻī penei:
E nānā hou aku kākou i nā kahua.
DF (Downlink Format, 5 bits) - hoʻoholo i ke ʻano o ka memo. Aia kekahi mau ʻano:
()
Makemake wale mākou i ke ʻano DF17, no ka mea... ʻO ia ka mea i loaʻa nā hoʻonohonoho o ka mokulele.
ICAO (24 bits) - code kū hoʻokahi honua o ka mokulele. Hiki iā ʻoe ke nānā i ka mokulele ma kāna code (akā, ua ho'ōki ka mea kākau i ka hoʻonui ʻana i ka waihona, akā pili nō naʻe). No ka laʻana, no ka code 3c5ee2 loaʻa iā mākou kēia ʻike:
Hoʻoponopono: in Hāʻawi ʻia ka wehewehe ʻana o ke code ICAO i nā kikoʻī hou aku; Paipai au i ka poʻe hoihoi e heluhelu.
ʻIkepili (56 a i ʻole 112 bits) - ka ʻikepili maoli a mākou e hoʻololi ai. ʻO nā 5 bits mua o ka ʻikepili ke kahua Kuhi helu, i loko o ka subtype o ka ʻikepili i mālama ʻia (ʻaʻole e huikau me DF). He kakaikahi o keia mau ano:
()
E nānā kākou i kekahi mau laʻana o nā pūʻolo.
ʻIke mokulele
Laʻana ma ke ʻano binary:
00100 011 000101 010111 000111 110111 110001 111000
Nā kahua ʻikepili:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, loaʻa i kēlā me kēia ʻano C1-C8 nā code e pili ana i nā helu helu ma ka laina:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ######_################0123456789######
Ma ka wehe ʻana i ke kaula, maʻalahi ka loaʻa ʻana o ke code mokulele: EWG7184
symbols = "#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######"
code_str = ""
for p in range(8):
c = int(bits_str[8 + 6*p:8 + 6*(p + 1)], 2)
code_str += symbols[c]
print("Aircraft Identification:", code_str.replace('#', ''))
Kūlana ea
Inā maʻalahi ka inoa, a laila ʻoi aku ka paʻakikī o nā hoʻonohonoho. Hoʻouna ʻia lākou ma ke ʻano o 2, ʻoiai a me nā papa ʻokoʻa. Ka helu kahua TC = 01011b = 11.
Ka laʻana o nā ʻeke like a me nā ʻeke ʻokoʻa:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
ʻO ka helu ʻana o nā hoʻonohonoho ponoʻī e like me kahi ʻano paʻakikī:
()
ʻAʻole wau he loea GIS, no laila ʻaʻole maopopo iaʻu mai hea mai ia. ʻO wai ka mea ʻike, e kākau i nā manaʻo.
ʻOi aku ka maʻalahi o ke kiʻekiʻe - ma muli o ke kiko kikoʻī, hiki ke hōʻike ʻia ma ke ʻano he nui o 25 a i ʻole 100 kapuaʻi.
Hololele Lewa
Pūʻolo me TC=19. ʻO ka mea hoihoi ma ʻaneʻi, ʻo ia ka hiki ke pololei ka māmā, pili i ka honua (Ground Speed), a i ʻole ea, i ana ʻia e ka mea ʻike mokulele (Airspeed). Hoʻouna ʻia nā ʻāpana like ʻole he nui:
()
hopena
E like me kāu e ʻike ai, ua lilo ka ʻenehana ADS-B i kahi symbiosis hoihoi, i ka wā e pono ai kahi maʻamau ʻaʻole wale i nā poʻe loea, akā i nā mea hoʻohana maʻamau. Akā ʻoiaʻiʻo, ua hoʻokani ʻia kahi hana koʻikoʻi e ka ʻenehana haʻahaʻa loa o nā mea hoʻokipa SDR digital, e hiki ai i ka hāmeʻa ke loaʻa maoli nā hōʻailona me nā alapine ma luna o kahi gigahertz "no nā peni."
Ma ka maʻamau pono'ī,ʻoiaʻiʻo, he nui aku. Hiki i ka poʻe makemake ke nānā i ka PDF ma ka ʻaoʻao a i ʻole e kipa i ka mea i ʻōlelo ʻia ma luna .
ʻAʻole hiki ke hoʻohana ʻia nā mea a pau i luna aʻe i nā mea he nui, akā ma ka liʻiliʻi o ka manaʻo nui o ke ʻano o ka hana ʻana, ke manaʻolana nei au.
Ma ke ala, aia kahi decoder mākaukau ma Python, hiki iā ʻoe ke aʻo . A hiki i nā mea hoʻokipa SDR ke hōʻuluʻulu a hoʻomaka i kahi decoder ADS-B i mākaukau , ua kūkākūkā nui ʻia kēia ma .
Hāʻawi ʻia ka code kumu o ka parser i wehewehe ʻia ma ka ʻatikala ma lalo o ka ʻoki. He laʻana hoʻāʻo kēia ʻaʻole i hoʻohālikelike ʻia i ka hana ʻana, akā hana kekahi mau mea i loko, a hiki ke hoʻohana ʻia e parse i ka faila i hoʻopaʻa ʻia ma luna.
Kumu kumu (Python)
from __future__ import print_function
from scipy.io import wavfile
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import math
import sys
def parse_message(data, start, bit_len):
max_len = bit_len*128
A = data[start:start + max_len]
A = signal.resample(A, 10*max_len)
bits = np.zeros(10*max_len)
bit_len *= 10
start_data = bit_len*8
# Parse first 8 bits
bits_str = ""
for p in range(8):
pos = start_data + bit_len*p
p1, p2 = A[pos: pos + bit_len/2], A[pos + bit_len/2: pos + bit_len]
avg1, avg2 = np.average(p1), np.average(p2)
if avg1 < avg2:
bits_str += "0"
elif avg1 > avg2:
bits_str += "1"
df = int(bits_str[0:5], 2)
# Aircraft address (db - https://junzis.com/adb/?q=3b1c5c )
bits_str = ""
for p in range(8, 32):
pos = start_data + bit_len * p
p1, p2 = A[pos: pos + bit_len / 2], A[pos + bit_len / 2: pos + bit_len]
avg1, avg2 = np.average(p1), np.average(p2)
if avg1 < avg2:
bits_str += "0"
elif avg1 > avg2:
bits_str += "1"
# print "Aircraft address:", bits_str, hex(int(bits_str, 2))
address = hex(int(bits_str, 2))
# Filter specific aircraft (optional)
# if address != "0x3c5ee2":
# return
if df == 16 or df == 17 or df == 18 or df == 19 or df == 20 or df == 21:
# print "Pos:", start, "DF:", msg_type
# Data (56bit)
bits_str = ""
for p in range(32, 88):
pos = start_data + bit_len*p
p1, p2 = A[pos: pos + bit_len/2], A[pos + bit_len/2: pos + bit_len]
avg1, avg2 = np.average(p1), np.average(p2)
if avg1 < avg2:
bits_str += "0"
# bits[pos + bit_len / 2] = 50
elif avg1 > avg2:
bits_str += "1"
# http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html
# print "Data:"
# print bits_str[:8], bits_str[8:20], bits_str[20:22], bits_str[22:22+17], bits_str[39:39+17]
# Type Code:
tc, ec = int(bits_str[:5], 2), int(bits_str[5:8], 2)
# print("DF:", df, "TC:", tc)
# 1 - 4 Aircraft identification
# 5 - 8 Surface position
# 9 - 18 Airborne position (w/ Baro Altitude)
# 19 Airborne velocities
if tc >= 1 and tc <= 4: # and (df == 17 or df == 18):
print("Aircraft address:", address)
print("Data:")
print(bits_str[:8], bits_str[8:14], bits_str[14:20], bits_str[20:26], bits_str[26:32], bits_str[32:38], bits_str[38:44])
symbols = "#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######"
code_str = ""
for p in range(8):
c = int(bits_str[8 + 6*p:8 + 6*(p + 1)], 2)
code_str += symbols[c]
print("Aircraft Identification:", code_str.replace('#', ''))
print()
if tc == 11:
print("Aircraft address:", address)
print("Data: (11)")
print(bits_str[:8], bits_str[8:20], bits_str[20:22], bits_str[22:22+17], bits_str[39:39+17])
# Bit 22 contains the F flag which indicates which CPR format is used (odd or even)
# First frame has F flag = 0 so is even and the second frame has F flag = 1 so odd
# f = bits_str[21:22]
# print("F:", int(f, 2))
# Altitude
alt1b = bits_str[8:20]
if alt1b[-5] == '1':
bits = alt1b[:-5] + alt1b[-4:]
n = int(bits, 2)
alt_ft = n*25 - 1000
print("Alt (ft)", alt_ft)
# lat_dec = int(bits_str[22:22+17], 2)
# lon_dec = int(bits_str[39:39+17], 2)
# print("Lat/Lon:", lat_dec, lon_dec)
# http://airmetar.main.jp/radio/ADS-B%20Decoding%20Guide.pdf
print()
if tc == 19:
print("Aircraft address:", address)
print("Data:")
# print(bits_str)
print(bits_str[:5], bits_str[5:8], bits_str[8:10], bits_str[10:13], bits_str[13] ,bits_str[14:24], bits_str[24], bits_str[25:35], bits_str[35:36], bits_str[36:65])
subtype = int(bits_str[5:8], 2)
# https://mode-s.org/decode/adsb/airborne-velocity.html
spd, hdg, rocd = -1, -1, -1
if subtype == 1 or subtype == 2:
print("Velocity Subtype 1: Ground speed")
v_ew_sign = int(bits_str[13], 2)
v_ew = int(bits_str[14:24], 2) - 1 # east-west velocity
v_ns_sign = int(bits_str[24], 2)
v_ns = int(bits_str[25:35], 2) - 1 # north-south velocity
v_we = -1*v_ew if v_ew_sign else v_ew
v_sn = -1*v_ns if v_ns_sign else v_ns
spd = math.sqrt(v_sn*v_sn + v_we*v_we) # unit in kts
hdg = math.atan2(v_we, v_sn)
hdg = math.degrees(hdg) # convert to degrees
hdg = hdg if hdg >= 0 else hdg + 360 # no negative val
if subtype == 3:
print("Subtype Subtype 3: Airspeed")
hdg = int(bits_str[14:24], 2)/1024.0*360.0
spd = int(bits_str[25:35], 2)
vr_sign = int(bits_str[36], 2)
vr = int(bits_str[36:45], 2)
rocd = -1*vr if vr_sign else vr # rate of climb/descend
print("Speed (kts):", spd, "Rate:", rocd, "Heading:", hdg)
print()
# print()
def calc_coordinates():
def _cprN(lat, is_odd):
nl = _cprNL(lat) - is_odd
return nl if nl > 1 else 1
def _cprNL(lat):
try:
nz = 15
a = 1 - math.cos(math.pi / (2 * nz))
b = math.cos(math.pi / 180.0 * abs(lat)) ** 2
nl = 2 * math.pi / (math.acos(1 - a/b))
return int(math.floor(nl))
except:
# happens when latitude is +/-90 degree
return 1
def floor_(x):
return int(math.floor(x))
lat1b, lon1b, alt1b = "10111000111010011", "10000110111111000", "000101111001"
lat2b, lon2b, alt2b = "10010011101011100", "10000011000011011", "000101110111"
lat1, lon1, alt1 = int(lat1b, 2), int(lon1b, 2), int(alt1b, 2)
lat2, lon2, alt2 = int(lat2b, 2), int(lon2b, 2), int(alt2b, 2)
# 131072 is 2^17, since CPR lat and lon are 17 bits each
cprlat_even, cprlon_even = lat1/131072.0, lon1/131072.0
cprlat_odd, cprlon_odd = lat2/131072.0, lon2/131072.0
print(cprlat_even, cprlon_even)
j = floor_(59*cprlat_even - 60*cprlat_odd)
print(j)
air_d_lat_even = 360.0 / 60
air_d_lat_odd = 360.0 / 59
# Lat
lat_even = float(air_d_lat_even * (j % 60 + cprlat_even))
lat_odd = float(air_d_lat_odd * (j % 59 + cprlat_odd))
if lat_even >= 270:
lat_even = lat_even - 360
if lat_odd >= 270:
lat_odd = lat_odd - 360
# Lon
ni = _cprN(lat_even, 0)
m = floor_(cprlon_even * (_cprNL(lat_even)-1) - cprlon_odd * _cprNL(lat_even) + 0.5)
lon = (360.0 / ni) * (m % ni + cprlon_even)
print("Lat", lat_even, "Lon", lon)
# Altitude
# Q-bit (bit 48) indicates whether the altitude is encoded in multiples of 25 or 100 ft (0: 100 ft, 1: 25 ft)
# The value can represent altitudes from -1000 to +50175 ft.
if alt1b[-5] == '1':
bits = alt1b[:-5] + alt1b[-4:]
n = int(bits, 2)
alt_ft = n*25 - 1000
print("Alt (ft)", alt_ft)
fs, data = wavfile.read("adsb_20190311_191728Z_1090000kHz_RF.wav")
T = 1/fs
print("Sample rate %f MS/s" % (fs / 1e6))
print("Cnt samples %d" % len(data))
print("Duration: %f s" % (T * len(data)))
data = data.astype(float)
cnt = data.shape[0]
# Processing only part on file (faster):
# cnt = 10000000
# data = data[:cnt]
print("Processing I/Q...")
I, Q = data[:, 0], data[:, 1]
A = np.sqrt(I*I + Q*Q)
bits = np.zeros(cnt)
# To see scope without any processing, uncomment
# plt.plot(A)
# plt.show()
# sys.exit(0)
print("Extracting signals...")
pos = 0
avg = 200
msg_start = 0
# Find beginning of each signal
while pos < cnt - 16*1024:
# P1 - message start
while pos < cnt - 16*1024:
if A[pos] < avg and A[pos+1] > avg and pos - msg_start > 1000:
msg_start = pos
bits[pos] = 100
pos += 4
break
pos += 1
start1, start2, start3, start4 = msg_start, 0, 0, 0
# P2
while pos < cnt - 16*1024:
if A[pos] < avg and A[pos+1] > avg:
start2 = pos
bits[pos] = 90
pos += 1
break
pos += 1
# P3
while pos < cnt - 16*1024:
if A[pos] < avg and A[pos+1] > avg:
start3 = pos
bits[pos] = 80
pos += 1
break
pos += 1
# P4
while pos < cnt - 16*1024:
if A[pos] < avg and A[pos+1] > avg:
start4 = pos
bits[pos] = 70
pos += 1
break
pos += 1
sig_diff = start4 - start1
if 20 < sig_diff < 25:
bits[msg_start] = 500
bit_len = int((start4 - start1) / 4.5)
# print(pos, start1, start4, ' - ', bit_len)
# start = start1 + 8*bit_len
parse_message(A, msg_start, bit_len)
pos += 450
# For debugging: check signal start
# plt.plot(A)
# plt.plot(bits)
# plt.show()
Manaʻo wau ua hoihoi kekahi, mahalo no kou nānā ʻana.
Source: www.habr.com