Mhoro Habr. Pamwe munhu wese akambosangana kana kuona hama kana shamwari dziri mundege ashandisa mahara Flightradar24 service. Iyi inzira iri nyore yekutevera nzvimbo yendege munguva chaiyo.
В Nheyo yekushanda yebasa rakadaro repaIndaneti yakatsanangurwa. Isu tichaenda kumberi toona kuti ndeipi data iri kutumirwa nekugamuchirwa kubva kundege kuenda kunzvimbo yekugamuchira uye nekuzvinyora isu pachedu tichishandisa Python.
История
Zviripachena, data rendege harifambiswe kuti vashandisi vaone pamafoni avo. Iyo sisitimu inodaidzwa kuti ADS-B (Automatic dependent surveillance-broadcast), uye inoshandiswa kuendesa otomatiki ruzivo nezve ndege kunzvimbo yekudzora - identifier yayo, marongero, kutungamira, kumhanya, urefu uye imwe data inofambiswa. Kare, kusati kwasvika masisitimu akadaro, iyo dispatcher yaingoona poindi pane radar. Izvi zvakanga zvisisakwani kana ndege dzakawandisa.
Nehunyanzvi, ADS-B ine inotapurirana mundege iyo nguva nenguva inotumira mapaketi eruzivo pane yakakwira frequency ye1090 MHz (kune mamwe modhi, asi isu hatifarire mazviri, sezvo marongero anofambiswa pano chete). Ehe, mukuwedzera kune inotapurirana, kune zvakare inogamuchira kumwe kunhu panhandare yendege, asi kwatiri isu sevashandisi, yedu inogamuchira inonakidza.
Sezvineiwo, kuenzanisa, yekutanga system yakadaro, Airnav Radarbox, yakagadzirirwa vashandisiwo zvavo, yakaonekwa muna 2007, uye inodhura madhora mazana mapfumbamwe; kunyorera kune network masevhisi kunodhura imwe $900 pagore.
Ongororo yeavo varidzi vekutanga veRussia inogona kuverengerwa paforamu . Zvino izvo RTL-SDR vanogamuchira vave kuwanikwa zvakanyanya, mudziyo wakafanana unogona kuunganidzwa nemadhora makumi matatu; zvimwe nezve izvi zvaive mukati. . Ngatienderere mberi kune protocol pachayo - ngatione kuti inoshanda sei.
Kugamuchira zviratidzo
Kutanga, chiratidzo chinoda kurekodhwa. Chiratidzo chese chine nguva ye120 microseconds chete, saka kuti iparadze zvakanaka zvikamu zvayo, SDR inogamuchira ine sampling frequency inokwana 5 MHz inodiwa.
Mushure mekurekodha, tinogamuchira WAV faira ine sampling rate ye5000000 samples / sec; 30 masekonzi ekurekodha akadaro "kuyera" nezve 500MB. Kuiteerera neine media player, hongu, hazvina maturo - iyo faira haina ruzha, asi yakananga digitized redhiyo chiratidzo - izvi ndizvo chaizvo mashandiro anoita Software Defined Radio.
Tichavhura nekugadzirisa faira tichishandisa Python. Avo vanoda kuedza ivo pachavo vanogona kudhawunirodha muenzaniso wekurekodha .
Ngatitorei faira uye tione zviri mukati.
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()
Mhedzisiro: tinoona "pulses" ari pachena achipesana neruzha rwekumashure.
Imwe neimwe "pulse" chiratidzo, chimiro chayo chinonyatso kuoneka kana iwe ukawedzera kugadzirisa pane girafu.
Sezvauri kuona, mufananidzo wacho unonyatsoenderana nezvinopihwa murondedzero iri pamusoro. Unogona kutanga kugadzirisa data.
Decoding
Kutanga, iwe unofanirwa kuwana zvishoma rukova. Chiratidzo pachacho chakavharidzirwa uchishandisa Manchester encoding:
Kubva pamusiyano wenhanho mune nibbles zviri nyore kuwana chaiyo "0" uye "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"
Chimiro chechiratidzo pachacho ndeichi:
Ngatitarisei minda mune zvakadzama.
DF (Downlink Format, 5 bits) - inosarudza rudzi rwemashoko. Kune marudzi akawanda:
()
Isu tinongofarira mhando DF17, nekuti... Izvi ndizvo zvine zvinorongeka zvendege.
ICAO (24 bits) - kodhi yepasirese yakasarudzika yendege. Unogona kutarisa ndege nekodhi yayo (zvinosuruvarisa, munyori akamira kugadzirisa dhatabhesi, asi ichiri kukosha). Semuenzaniso, yekodhi 3c5ee2 isu tine inotevera ruzivo:
Rongedza: mu Rondedzero yeiyo ICAO kodhi inopihwa zvakadzama; Ini ndinokurudzira kuti avo vanofarira vaverenge.
dhata (56 kana 112 bits) - iyo chaiyo data yatichasarudza. Yekutanga 5 bits yedata ndiyo munda Type Code, ine subtype yedata iri kuchengetwa (kwete kuti ivhiringidzwe neDF). Kune akati wandei emhando idzi:
()
Ngatitarisei mienzaniso mishoma yemapakeji.
Kuzivikanwa kwendege
Muenzaniso mune binary fomu:
00100 011 000101 010111 000111 110111 110001 111000
Data fields:
+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+
TC = 00100b = 4, hunhu hwega hwega C1-C8 ine macode anoenderana neindices mumutsetse:
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_#################0123456789#######
Nekukodha tambo, zviri nyore kuwana iyo ndege kodhi: 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('#', ''))
Airborne position
Kana zita racho riri nyore, saka zvikorekedzo zvakanyanya kuoma. Iwo anotapurirwa muchimiro che2, kunyange uye asinganzwisisike mafuremu. Munda kodhi TC = 01011b = 11.
Muenzaniso wemapakiti akaenzana uye asinganzwisisike:
01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000
Iko kuverengerwa kwemacoordinates pachayo kunoitika zvinoenderana neyakaoma fomula:
()
Ini handisi nyanzvi yeGIS, saka handizive kuti inobva kupi. Ndiani anoziva, nyora mumashoko.
Hurefu hunoonekwa sehukareruka - zvichienderana nechaiyo bhiti, inogona kumiririrwa senge yakapetwa yemakumi maviri nemashanu kana zana tsoka.
Airborne Velocity
Pakeji ine TC = 19. Chinhu chinonakidza apa ndechekuti kumhanya kunogona kunge kwakaringana, kunoenderana nepasi (Ground Speed), kana kubhururuka nendege, kuyerwa nesensor yendege (Airspeed). Nzvimbo dzakawanda dzakasiyana dzinofambiswa zvakare:
()
mhedziso
Sezvauri kuona, tekinoroji yeADS-B yave inonakidza symbiosis, apo chiyero chinobatsira kwete kune nyanzvi chete, asiwo kune vese vashandisi. Asi hongu, basa rakakosha mune izvi rakaitwa neyakachipa tekinoroji yedhijitari SDR inogamuchira, iyo inobvumira mudziyo kuti ugamuchire masaini nema frequency pamusoro pegigahertz "yemapeni."
Mune mwero pachawo, hongu, kune zvimwe zvakawanda. Avo vanofarira vanogona kuona PDF iri papeji kana kuti shanyira yatotaurwa pamusoro apa .
Hazvigoneke kuti zvese zviri pamusoro zvichabatsira kune vakawanda, asi zvirinani pfungwa yakajairika yekuti inoshanda sei, ndinovimba, inoramba iripo.
Nenzira, yakagadzirira-yakagadzirwa decoder muPython yatovepo, unogona kuidzidza . Uye varidzi vevanogamuchira SDR vanogona kuungana uye kuvhura yakagadzirira-yakagadzirwa ADS-B decoder , izvi zvakakurukurwa zvakadzama mu .
Iko kunobva kodhi yeparser inotsanangurwa muchinyorwa inopiwa pazasi pekucheka. Uyu muenzaniso wekuyedza usingaite kunge kugadzirwa, asi zvimwe zvinhu zvinoshanda mauri, uye unogona kushandiswa kuburitsa faira rakanyorwa pamusoro.
Kunobva kodhi (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()
Ndinovimba mumwe munhu anga achifarira, ndatenda nekutarisisa kwako.
Source: www.habr.com