Flightradar24 - តើវាដំណើរការយ៉ាងដូចម្តេច? ផ្នែកទី 2 ពិធីការ ADS-B

សួស្តី Habr ។ ប្រហែល​ជា​អ្នក​រាល់​គ្នា​ដែល​ធ្លាប់​ជួប ឬ​ឃើញ​សាច់​ញាតិ ឬ​មិត្ត​ភក្តិ​នៅ​លើ​យន្ត​ហោះ​បាន​ប្រើ​សេវា Flightradar24 ដោយ​ឥត​គិត​ថ្លៃ។ នេះគឺជាមធ្យោបាយដ៏ងាយស្រួលបំផុតក្នុងការតាមដានទីតាំងរបស់យន្តហោះក្នុងពេលជាក់ស្តែង។

В ផ្នែកដំបូង។ គោលការណ៍ប្រតិបត្តិការនៃសេវាកម្មអនឡាញបែបនេះត្រូវបានពិពណ៌នា។ ឥឡូវនេះ យើងនឹងបន្តស្វែងយល់ថាតើទិន្នន័យអ្វីខ្លះដែលត្រូវបានផ្ញើ និងទទួលពីយន្តហោះទៅកាន់ស្ថានីយទទួល ហើយធ្វើការឌិកូដវាដោយខ្លួនឯងដោយប្រើ Python ។

История

ជាក់ស្តែង ទិន្នន័យយន្តហោះមិនត្រូវបានបញ្ជូនឱ្យអ្នកប្រើប្រាស់មើលឃើញនៅលើស្មាតហ្វូនរបស់ពួកគេនោះទេ។ ប្រព័ន្ធនេះត្រូវបានគេហៅថា ADS-B (ការឃ្លាំមើលដោយស្វ័យប្រវត្ត-ផ្សាយ) ហើយត្រូវបានប្រើដើម្បីបញ្ជូនព័ត៌មានអំពីយន្តហោះដោយស្វ័យប្រវត្តិទៅកាន់មជ្ឈមណ្ឌលបញ្ជា - ឧបករណ៍កំណត់អត្តសញ្ញាណ កូអរដោនេ ទិសដៅ ល្បឿន រយៈកម្ពស់ និងទិន្នន័យផ្សេងទៀតរបស់វាត្រូវបានបញ្ជូន។ ពីមុន មុនពេលការមកដល់នៃប្រព័ន្ធបែបនេះ អ្នកបញ្ជូនអាចមើលឃើញតែចំណុចមួយនៅលើរ៉ាដាប៉ុណ្ណោះ។ នេះមិនគ្រប់គ្រាន់ទៀតទេ នៅពេលដែលមានយន្តហោះច្រើនពេក។

តាមបច្ចេកទេស ADS-B មានឧបករណ៍បញ្ជូននៅលើយន្តហោះដែលបញ្ជូនកញ្ចប់ព័ត៌មានតាមកាលកំណត់នៅប្រេកង់ខ្ពស់ 1090 MHz (មានរបៀបផ្សេងទៀត ប៉ុន្តែយើងមិនសូវចាប់អារម្មណ៍នឹងពួកវាទេ ព្រោះកូអរដោនេត្រូវបានបញ្ជូនតែនៅទីនេះ)។ ជាការពិតណាស់ បន្ថែមពីលើឧបករណ៍បញ្ជូន ក៏មានអ្នកទទួលនៅកន្លែងណាមួយនៅអាកាសយានដ្ឋានផងដែរ ប៉ុន្តែសម្រាប់យើងជាអ្នកប្រើប្រាស់ អ្នកទទួលផ្ទាល់របស់យើងគឺគួរឱ្យចាប់អារម្មណ៍។

ដោយវិធីនេះ សម្រាប់ការប្រៀបធៀប ប្រព័ន្ធបែបនេះដំបូងគឺ Airnav Radarbox ដែលត្រូវបានរចនាឡើងសម្រាប់អ្នកប្រើប្រាស់ធម្មតាបានបង្ហាញខ្លួននៅក្នុងឆ្នាំ 2007 ហើយមានតម្លៃប្រហែល 900 ដុល្លារ ហើយការជាវសេវាបណ្តាញមានតម្លៃ 250 ដុល្លារទៀតក្នុងមួយឆ្នាំ។

ការពិនិត្យឡើងវិញនៃម្ចាស់ជនជាតិរុស្ស៊ីដំបូងគេអាចអានបាននៅលើវេទិកា ម៉ាស៊ីនស្កេនវិទ្យុ. ឥឡូវនេះអ្នកទទួល RTL-SDR បានក្លាយជាអាចប្រើបានយ៉ាងទូលំទូលាយ ឧបករណ៍ស្រដៀងគ្នានេះអាចត្រូវបានផ្គុំក្នុងតម្លៃ 30 ដុល្លារ។ បន្ថែមទៀតអំពីនេះគឺនៅក្នុង ផ្នែកដំបូង។. ចូរបន្តទៅពិធីការដោយខ្លួនឯង - តោះមើលរបៀបដែលវាដំណើរការ។

ការទទួលសញ្ញា

ទីមួយ សញ្ញាត្រូវកត់ត្រា។ សញ្ញាទាំងមូលមានថិរវេលាត្រឹមតែ 120 មីក្រូវិនាទីប៉ុណ្ណោះ ដូច្នេះដើម្បីងាយស្រួលបំបែកសមាសធាតុរបស់វា អ្នកទទួល SDR ដែលមានប្រេកង់គំរូយ៉ាងហោចណាស់ 5 MHz គឺជាការចង់បាន។

បន្ទាប់ពីការថត យើងទទួលបានឯកសារ WAV ដែលមានអត្រាគំរូ 5000000 គំរូ/វិនាទី; 30 វិនាទីនៃការថតបែបនេះ "ថ្លឹងថ្លែង" ប្រហែល 500MB ។ ការស្តាប់វាជាមួយកម្មវិធីចាក់មេឌៀ ពិតណាស់គឺគ្មានប្រយោជន៍ទេ - ឯកសារមិនមានសំឡេងទេ ប៉ុន្តែជាសញ្ញាវិទ្យុឌីជីថលផ្ទាល់ - នេះគឺជារបៀបដែលកម្មវិធីកំណត់វិទ្យុដំណើរការ។

យើងនឹងបើក និងដំណើរការឯកសារដោយប្រើ Python ។ អ្នកដែលចង់ពិសោធន៍ដោយខ្លួនឯងអាចទាញយកការថតឧទាហរណ៍ តំណភ្ជាប់.

តោះទាញយកឯកសារ ហើយមើលអ្វីដែលនៅខាងក្នុង។

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()

លទ្ធផល៖ យើងឃើញ "ជីពចរ" ជាក់ស្តែងប្រឆាំងនឹងសំលេងរំខានផ្ទៃខាងក្រោយ។

"ជីពចរ" នីមួយៗគឺជាសញ្ញាមួយដែលរចនាសម្ព័ន្ធអាចមើលឃើញយ៉ាងច្បាស់ប្រសិនបើអ្នកបង្កើនគុណភាពបង្ហាញនៅលើក្រាហ្វ។

ដូចដែលអ្នកអាចឃើញរូបភាពគឺពិតជាស្របទៅនឹងអ្វីដែលត្រូវបានផ្តល់ឱ្យនៅក្នុងការពិពណ៌នាខាងលើ។ អ្នកអាចចាប់ផ្តើមដំណើរការទិន្នន័យ។

ការឌិកូដ

ដំបូងអ្នកត្រូវទទួលបានស្ទ្រីមបន្តិច។ សញ្ញាខ្លួនវាត្រូវបានអ៊ិនកូដដោយប្រើការអ៊ិនកូដ Manchester៖

ពីភាពខុសគ្នានៃកម្រិតនៅក្នុង nibbles វាងាយស្រួលក្នុងការទទួលបាន "0" និង "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"

រចនាសម្ព័ន្ធនៃសញ្ញាខ្លួនឯងមានដូចខាងក្រោម:

សូមក្រឡេកមើលវាលឱ្យបានលំអិត។

DF (ទម្រង់តំណខាងក្រោម 5 ប៊ីត) - កំណត់ប្រភេទនៃសារ។ មានច្រើនប្រភេទ៖

(ប្រភពតារាង)

យើងចាប់អារម្មណ៍តែប្រភេទ DF17 ទេព្រោះ... វាគឺជានេះដែលមានកូអរដោនេនៃយន្តហោះ។

ICAO (24 ប៊ីត) - លេខកូដអន្តរជាតិនៃយន្តហោះ។ អ្នកអាចពិនិត្យមើលយន្តហោះដោយលេខកូដរបស់វា។ លើបណ្តាញ (ជាអកុសល អ្នកនិពន្ធបានឈប់ធ្វើបច្ចុប្បន្នភាពមូលដ្ឋានទិន្នន័យ ប៉ុន្តែវានៅតែពាក់ព័ន្ធ)។ ឧទាហរណ៍ សម្រាប់លេខកូដ 3c5ee2 យើងមានព័ត៌មានដូចខាងក្រោម៖

កែសម្រួល៖ ក្នុង យោបល់ទៅអត្ថបទ ការពិពណ៌នានៃកូដ ICAO ត្រូវបានផ្តល់ឱ្យយ៉ាងលម្អិត ខ្ញុំសូមណែនាំឱ្យអ្នកដែលចាប់អារម្មណ៍អានវា។

ទិន្នន័យ (56 ឬ 112 ប៊ីត) - ទិន្នន័យពិតប្រាកដដែលយើងនឹងឌិកូដ។ ទិន្នន័យ 5 ប៊ីតដំបូងគឺជាវាល ប្រភេទកូដដែលមានប្រភេទរងនៃទិន្នន័យដែលត្រូវបានរក្សាទុក (មិនត្រូវច្រឡំជាមួយ DF)។ មានប្រភេទទាំងនេះមួយចំនួន៖

(ប្រភពតារាង)

សូមក្រឡេកមើលឧទាហរណ៍មួយចំនួននៃកញ្ចប់។

ការកំណត់អត្តសញ្ញាណយន្តហោះ

ឧទាហរណ៍ក្នុងទម្រង់គោលពីរ៖

00100 011 000101 010111 000111 110111 110001 111000

វាលទិន្នន័យ៖

+------+------+------+------+------+------+------+------+------+------+
| TC,5 | EC,3 | C1,6 | C2,6 | C3,6 | C4,6 | C5,6 | C6,6 | C7,6 | C8,6 |
+------+------+------+------+------+------+------+------+------+------+

TC = 00100b = 4 តួអក្សរនីមួយៗ C1-C8 មានលេខកូដដែលត្រូវគ្នានឹងសន្ទស្សន៍ក្នុងបន្ទាត់៖
#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_##################0123456789######

ដោយការឌិកូដខ្សែអក្សរ វាងាយស្រួលក្នុងការទទួលបានលេខកូដយន្តហោះ៖ 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('#', ''))

ទីតាំងខ្យល់

ប្រសិនបើឈ្មោះសាមញ្ញ នោះកូអរដោណេកាន់តែស្មុគស្មាញ។ ពួកវាត្រូវបានបញ្ជូនក្នុងទម្រង់ជា 2 ស៊ុមគូ និងសេស។ លេខកូដវាល TC = 01011b = 11 ។

ឧទាហរណ៍នៃកញ្ចប់គូ និងសេស៖

01011 000 000101110110 00 10111000111001000 10000110101111001
01011 000 000110010000 01 10010011110000110 10000011110001000

ការ​គណនា​កូអរដោណេ​ដោយ​ខ្លួន​វា​កើត​ឡើង​ដោយ​យោង​តាម​រូបមន្ត​ដែល​ពិបាក​ជាង​នេះ៖

(ប្រភព)

ខ្ញុំមិនមែនជាអ្នកជំនាញ GIS ទេ ដូច្នេះខ្ញុំមិនដឹងថាវាមកពីណាទេ។ អ្នកណាដឹងសរសេរក្នុងមតិយោបល់។

កម្ពស់ត្រូវបានគេចាត់ទុកថាសាមញ្ញជាង - អាស្រ័យលើប៊ីតជាក់លាក់វាអាចត្រូវបានតំណាងថាជាពហុគុណនៃ 25 ឬ 100 ហ្វីត។

ល្បឿនខ្យល់

កញ្ចប់ជាមួយ TC = 19 ។ អ្វី​ដែល​គួរ​ឱ្យ​ចាប់​អារម្មណ៍​នៅ​ទី​នេះ​គឺ​ថា ល្បឿន​អាច​មាន​ភាព​សុក្រឹត ទាក់ទង​ទៅ​នឹង​ដី (Ground Speed) ឬ​ខ្យល់​ដែល​វាស់​ដោយ​ឧបករណ៍​ចាប់​សញ្ញា​យន្តហោះ (Airspeed)។ វាលផ្សេងគ្នាជាច្រើនត្រូវបានបញ្ជូនផងដែរ:

(ប្រភព)

សេចក្តីសន្និដ្ឋាន

ដូចដែលអ្នកអាចឃើញ បច្ចេកវិទ្យា ADS-B បានក្លាយជានិមិត្តសញ្ញាដ៏គួរឱ្យចាប់អារម្មណ៍មួយ នៅពេលដែលស្តង់ដារមួយមានប្រយោជន៍មិនត្រឹមតែសម្រាប់អ្នកជំនាញប៉ុណ្ណោះទេ ប៉ុន្តែក៏សម្រាប់អ្នកប្រើប្រាស់ធម្មតាផងដែរ។ ប៉ុន្តែជាការពិតណាស់ តួនាទីសំខាន់ក្នុងរឿងនេះត្រូវបានលេងដោយបច្ចេកវិទ្យាថោកជាងនៃអ្នកទទួលឌីជីថល SDR ដែលអនុញ្ញាតឱ្យឧបករណ៍នេះទទួលបានសញ្ញាតាមព្យញ្ជនៈជាមួយនឹងប្រេកង់លើសពី gigahertz "សម្រាប់កាក់" ។

នៅក្នុងស្តង់ដារខ្លួនវាពិតណាស់មានច្រើនទៀត។ អ្នកដែលចាប់អារម្មណ៍អាចមើលឯកសារ PDF នៅលើទំព័រ ICAO ឬចូលទៅកាន់គេហទំព័រដែលបានរៀបរាប់ខាងលើ វេបសាយ.

វាមិនទំនងថាទាំងអស់ខាងលើនឹងមានប្រយោជន៍សម្រាប់មនុស្សជាច្រើននោះទេប៉ុន្តែយ៉ាងហោចណាស់គំនិតទូទៅនៃរបៀបដែលវាដំណើរការខ្ញុំសង្ឃឹមថានៅតែមាន។

ដោយវិធីនេះ ឧបករណ៍ឌិកូដដែលត្រៀមរួចជាស្រេចនៅក្នុង Python មានរួចហើយ អ្នកអាចសិក្សាវាបាន នៅទីនេះ. ហើយម្ចាស់អ្នកទទួល SDR អាចប្រមូលផ្តុំ និងបើកដំណើរការឧបករណ៍ឌិកូដ ADS-B ដែលត្រៀមរួចជាស្រេច ពីទំព័រនេះត្រូវបានពិភាក្សាលម្អិតបន្ថែមទៀតនៅក្នុង ផ្នែកដំបូង។.

កូដប្រភពនៃធាតុញែកដែលបានពិពណ៌នានៅក្នុងអត្ថបទត្រូវបានផ្តល់ឱ្យខាងក្រោមការកាត់។ នេះ​ជា​ឧទាហរណ៍​សាកល្បង​ដែល​មិន​ធ្វើ​ពុត​ជា​ផលិត​ទេ ប៉ុន្តែ​រឿង​ខ្លះ​ដំណើរការ​ក្នុង​វា ហើយ​វា​អាច​ប្រើ​ដើម្បី​ញែក​ឯកសារ​ដែល​បាន​កត់ត្រា​ខាងលើ។
កូដប្រភព (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()

ខ្ញុំសង្ឃឹមថាមាននរណាម្នាក់ចាប់អារម្មណ៍ សូមអរគុណចំពោះការយកចិត្តទុកដាក់របស់អ្នក។

ប្រភព: www.habr.com