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