Šiuo savo įrašu iš dalies noriu pagerbti savo prieš 20 (ar pan.) metų buvusius signalų apdorojimo dėstytojus Miniją Tamošiūnaitę ir a.a. Adolfą Laimutį Telksnį. Abu tuo metu buvo visiškai nerealūs ir jų dėstomi dalykai buvo labai įdomūs. Ne viską anuomet supratau ir mažai kas užsiliko, bet darbas su Matlab’u man iki šiol kelia daug smagių prisiminimų.
O dabar prie reikalo. Pradžiai, tai ši rašliava yra pažintinė. Aš su visais žemiau išpaistytais dalykais nesiruošiu sukurti jokio galutinio „produkto“, o tik patyrinėti signalų apdorojimą. Yra jau prikurta programinės įrangos, kuri viską daro kur kas geriau ir kokybiškiau – apie tai parašysiu vėliau. O šiam kartui noriu tiesiog pradėti nuo nulio ir pasigilinti į visokias matematines subtilybes.
Digital Selective Calling – trumpų skaitmeninių žinučių standartas, skirtų laivybos radijo sistemoms, veikiančioms vidutinių, trumpųjų ir ultratrumpųjų bangų diapazonuose. Tai visiškai atviras standartas. Deja, atviro kodo dekoderių nėra. Nu gal ne visai, bet nėra jokio, kuris būtu padarytas normaliai iki galo. Taip pat nėra nė vieno Linuxui. Kadangi standartas gan nesudėtingas, nusprendžiau bent jau kažkokį kreivą-šleivą dekoderį pasigamint. Yra windauzams labai geras nemokamas DSC dekoderis, pavadinimu YADD (Yet Another DSC Decoder). Bet aš nesuprantu, kodėl autorius nenorėjo jo kodo atverti. Aš šitos „tradicijos“ nesuvokiu: pilna Windauzams nemokamo softo, bet atviro kodo – minimaliai. Ir va, YADD dabar yra abandonware, nes prieš kelerius metus jo autorius tiesiog numirė 😦 Patapo silent key. Kol dar su naujesniais windauzais YADD veikia, gerai. Bet gali vieną dieną ir nebeveikti.
Truputį pezalų
DSC vidutinėmis ir trumposiomis bangomis siunčiamas 100 dvejetainių bodų (taigi gaunasi kad bitų) sparta. Moduliacija – FSK, t.y. dažnio moduliacijos (kaip FM). Nešėjas gali būti F1B (dažnio moduliacija nuo centrinio nešėjo) arba J2B (dažniausiai USB + 1700 Hz).
Mano karjeroje visiškai nėra buvę kokio nors darbo su signalais. Na, pažaidimai su mikrovaldikliais ir skaitmeniniais jutikliais nesiskaito – ten viskas paprasta, skaitai specifikaciją ir darai pagal ją. Tuo tarpu signalų apdorojimas yra visai kita opera ir variantų taip pat yra krūvos.
Tad kai sugalvojau biškį tais signalų dekodavimais užsiimti, iš karto prisiminiau savo dėstytojus. Kai mokiausi anuomet, tai Minija buvo žiauriai graži merga ir bent jau man seilė tekėdavo sėdint jos paskaitose 😀 Deja, ir mokytis teko…
Ir, deja, iš tų laikų nelabai liko praktinių prisiminimų, ką galėčiau panaudoti. Tad į signalų dorojimą teko nerti kaip visiškai žaliam. Na, gal ne visiškai, su Python, Numpy ir Scipy šiek tiek patirties turiu, programuot bendrąja prasme moku. Bet pati logika, signalų pavertimas tuo, ko reikia, aplaižymai visokie – konkretus kosmosas. Plius, visada turėjau problemų su kompleksiniais skaičiais, o čia supratau, kad be jų neišsiversiu – bent jau minimaliai. Ir kas įdomiausia, pradėjęs su jais krapštytis, vartyti literatūrą, pagaliau, drąsiai įžengęs į penktąją dešimtį, pradėjau juos suvokti. Aikit sau, žinokit, kaip faina pasidarė! Na, bet labai daug matematikos su jais neprireikė.
DSC yra dvejetainis FSK signalas. Tas pats principas, kaip RTTY ar SITOR/AMTOR/NAVTEX ar netgi STANAG 4481. Vienas dažnis – nuliukas, kitas – vienetukas. Tarp jų – 170 Hz tarpas. Ir viskas. Bitukų kiekis yra tarpelio tarp dažnių pasikeitimo ilgis. Nėra jokio bitstuffingo. Tad ir vienas bodas čia lygus vienam bitui. Sakykim, jei signalas būtų iš daugiau dažnių, pavyzdžiui, keturių, tai jau vienas bodas galėtų būti 4 reikšmių (2 bitų kitaip sakant).
Mano pirminis uždavinukas – nuskaityt WAV failą su DSC žinutėmis ir jį demoduliuoti. Priklausomai nuo to, kaip tas failas įrašytas, reikia nusifiltruoti visus nereikalingus dažnius ir palikti tik DSC šniūrelį. O tada tą šniūrelį analizuoti ir bandyti „atstatyti“ originalų skaitmeninį signalą. Tam reikia ištisinį garso signalą paversti dažnio reikšmėmis – demoduliuoti.
Demoduliavimas
FSK (arba FM) demoduliavimas gali būti daromas keliais būdais. Sakykim, kad mūsų WAV signalas įrašytas gražiai, jame garso bangelės svyruodamos kerta nuliuką. Galima „bėgti“ per tą failą ir skaičiuoti tarpelius tarp kirtimų. Dar kažkiek paapvalinti ir bandyti iš tų tarpelių paskaičiuoti dažnį. Matematiškai paprasta, galutinių duomenų kiekis gaunasi mažesnis už originalą (nes tarpelių ilgiai užima mažiau vietos), bet pametama laiko dedamoji, ją reikia paskui išsiskaičiuoti atgal. Plius, šitas būdas gerai tinka švariems diskrečių dažnių signalams. Jei signalas triukšmingesnis, arba jame bangelės yra kelių dažnių miksas – šūds gausis.
Kitas būdas – PLL algoritmas, kuriame programinis osciliatorius kabinamas ant to paties nuliuko kirtimo. Gan paprastai įgyvendinamas bet kokiomis programinėmis priemonėmis, kaip ir pirmasis, bet duoda tikslesnį rezultatą.
Dar galima remtis formule sin²x + cos²x = 1. Tiesiai šeriam sinusą ir kosinusą WAV failo skaičiukais po vieną, kas kelis padarom jų sumos kvadratų šaknį norimais dažniais ir atidedam. Gaunam du masyvus: vienas viršutinio dažnio, kitas – apatinio. Geras metodas, jei norim apdoroti kuo greičiau ir nekaupti buferyje garso amplitudės nuskaitymų. Bet visiškai nepasiteisina su triukšmingais signalais ir labai priklauso nuo pasirinkto „langelio“ šaknies traukimui.
Labai įdomus variantas dažnių analizei yra bangelės (wavelet’ai). Su jomis labai įdomias skaleogramas galima nupaišyti. Bet pabandęs supratau, kad dėl trumpučio signalo paleidinėt pusę minutės trunkančią analizę tiesiog neverta. Bet ten visą spektrą laiko atžvilgiu būtų galima analizuoti, labai įdomus variantas. Visgi, skirtas tiesiog ne tam.
Mano tyrinėti jau sukurti dekoderiai lupa iš signalo atskirus dažnius, kurie domina, ir juos išfiltruoja. Tada analizuoja lygiagrečius dažnių srautus. Tai, ko gero, tiksliausias metodas, bet aš jį nusprendžiau pasilikti realaus dekoderio gamybai. Tiksliau, copy/paste iš esamo kodo 😀
Paskutinis variantas, kuris man dar kažkiek suprantamas – analitinio signalo pasigaminimas (arba ėmimas tiesiai iš SDR) ir jo demoduliavimas. Pitono bibliotekos turi viską, kas reikalinga, galima demoduliaciją padaryti keliomis eilutėmis. O ko daugiau reikia tokiam signalų „ekspertui“, kaip aš? Priešpaskutinis dažnių išskyrimo metodas taip pat naudoja analitinį signalą.
Analitinis signalas
Kas tas analitinis signalas? Grubiai tariant, tai tokia matematinė koncepcija, kad su dviem sinuso ir kosinuso kreivėm galima sukurti bet kokią kitą sinusoidę anas visaip maišant. Aš nežinau, ar fiziškai siųstuvai naudoja šitą reikalą, ar ne, bet galim įsivaizduoti – kad taip. SDR tuo tarpu priimdamas signalą jį perleidžia per 90° keitiklį ir kompiuteriui duoda du skaičiukų srautus: IQ. I reiškia in phase, arba iš esmės realųjį signalą, o Q yra quadrature – kvadratūra, 90° ne fazėje esanti dedamoji. SDR’ai šitą bajerį atlieka su dviem ADC, vienas duoda I, kitas Q skaičiukus. Tai padeda išvengti Nyquisto problemos, kai nuskaitymo dažnis būna mažesnis už ateinantį dažnį, t.y. gali atsirasti visokie parazitiniai dažniai, signalų „vaiduokliai“ ir pan.
Analitinis signalas dar pasižymi tuo, kad yra su kompleksiniais skaičiukais. IQ srautas ir yra toks signalas. Tačiau WAV garso faile mes turime realų signalą. T.y. iš SDR IQ srauto pagamintą galutinį produktą. Kadangi mano tikslas yra būtent WAV failų skaitymas (o ne IQ), tai man būtų gerai iš jo kaip nors pasigaminti analitinį signalą. Nes turėtų būti įmanoma.
Ir yra įmanoma. Tokį signalą gali padaryti Hilberto transformacija.
Pirmas dalykas: gražus ir kokybiškas įrašytas DSC signalas, atsklidęs iš kažkur netoli, spėju, ir jo spektrogramos atvaizdas:
Tikrai labai gražus signalas, nors truputį yra jame patriukšmavimų, bet tokį galima praktiškai vizualiai „perskaityti“. Tiesa, FFT langelį reikėtų pamažinti – čia jis yra 1024 taškų 16000 skiriamosios gebos failui.
Taigi viskas prasidėjo nuo to, kad aš ilgokai ieškojau, kaip čia dabar FSK/FM demoduliavimą padaryti su Pitonu. Labai paprasti sprendimai neveikė, o labai sudėtingi – nu tiesiog sudėtingi, nesinorėjo man per daug su jais gaišti laiko į kiekvieną smulkmeną įsigilinant.
Pradėjau su PLL metodu. Jis man pasirodė visai neblogas, bet gautas demoduliuotas signalas buvo gan „banguotas“, su daug šiukšlių. Tada pradėjau kapstytis, kaip čia nufiltruoti viską, ko nereikia, o pasilikti tik 170 Hz pločio signalą.
Tam radau gerą resursą (nuoroda apačioje bus). Iš ten supratau vieną dalyką: skaitmeniniai filtrai kuriami +/- norimo pločio dažniams. T.y. jei mes turim DSC signalą, kuris yra 170 Hz pločio, jis nuo savo centrinio dažnio yra arba -85 arba +85 Hz. Filtras kuriamas būtent nuo „minus dažnio“ iki „plius dažnio“. Bet realiame gyvenime „minusinio“ dažnio nebūna. Tiksliau, būna, jei signalas yra analitinis arba IQ 🙂 Su kompleksiniais skaičiais yra įmanomas „neigiamas dažnis“. Bet kadangi aš noriu filtruoti realų signalą, tai procesas toks: sukuriam reikiamą filtrą ir jį „pastumiam“. Pastūmimo formulė paprasta:
Tik tiek, kad ji su menamu vienetuku. Gerai, kad Pitone ten viskas padaryta paprastai ir su menamais vienetukais nėr problemų.
Per filtrą perleistas signalas būna gerokai švaresnis ir lengviau apdorojamas.
Daug prirašiau, dabar prie praktikos. Kaip nuskaityti WAV failą, nerodysiu – ten viskas labai paprasta. Pirmiausia visi importai:
import sys
import wave
import struct
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.signal import hilbert, butter, lfilter
Aš nusiskaičiau visą 9 s ilgio WAV įrašą su DSC signalu į kintamąjį samples. Tada man patinka normalizuoti viską į (-1 .. 1) ribas ir, aišku, centruoti signalą aplink nulį – kad būtų galima rasti kirtimus rast, jei to reikia. Man to reikėjo pirmiems bandymams su PLL, kuriuos paskui mečiau.
Taigi normalizavimas. Pastaba: samples nėra paprastas Python masyvas, o numpy.array tipo. Dėl to su juo matematinės operacijos kode yra labai paprastos:
norm_samples = samples / max(samples)
norm_samples -= np.average(norm_samples)
Tada, kaip jau minėjau, filtro susikūrimas ir pastūmimas iki centrinio dažnio center_freq (mano įraše jis 1500 Hz), imant nuskaitymo gebą sample_rate, o period yra 1 / sample_rate:
dsc_filter = signal.firwin(101, 90, nyq = sample_rate / 2)
shift_time = np.arange(0.0, period * len(dsc_filter), period)
exponential = np.exp(2j*np.pi*(-1)*center_freq*shift_time)
dsc_filter = dsc_filter * exponential
101 čia yra filtro ilgis. Pats filtras, jei jį grafiškai nupaišytumėm, yra tiesiog bangelė, arba wavelet. Paveikslėlyje mėlynas yra originalus -90..+90 filtras, o oranžinis – pastumtas per 1500 Hz, t.y. 1410..1590 filtras:
Jis perleidžiamas per signalą konvoliucijos metodu ir gaunam apvalytą signalą:
norm_samples = np.convolve(norm_samples, dsc_filter, 'valid')
Kaip matote, kol kas viskas paprasta, kelios eilutės – ir pridaryta jau daug dalykų. Turime nufiltruotą signalą, nors vizualiai didelio skirtumo nuo originalaus lyg ir nėra:
O dabar jau laikas magijai. Reikia šitą nufiltruotą signalą paversti analitiniu ir pabandyti „išlupti“ iš jo momentinį dažnio pokytį. Momentinis dažnio pokytis turėtų parodyti būtent DSC signalo formą.
z_samples = hilbert(np.real(norm_samples))
z_inst_phase = np.unwrap(np.angle(z_samples))
z_inst_freq = np.diff(z_inst_phase) / (2 * np.pi) * sample_rate
Taigi iš normalizuoto signalo realiosios dalies darom Hilberto transformaciją. Kodėl realiosios? Nes kai filtrą pastumiam su kompleksiniais skaičiais, tų kompleksinių atsiranda ir nufiltruotame signale. Tada pasigaminam fazės pokyčio signalą, kuris kompleksinėje erdvėje skaičiuojamas kaip kampas. Na, ir galiausiai perleidžiam tą fazės pokytį per difereciatorių.
Rezultatas va toks:
Jau akivaizdu, kad tokį signalą galima analizuoti. Čia yra beveik atgamintas originalusis DSC signalas – galite palyginti su spektrograma pradžioje. Oranžinė linija – centrinis 1500 Hz dažnis. Tiesa, net ir filtruotame signale yra triukšmų, papildomai šitą signalą galima dar „aplaižyti“ su Butterworth filtru.
b, a = butter(N=5, Wn=85, fs=sample_rate, btype="lowpass", analog=False)
z_inst_freq_filt = lfilter(b, a, z_inst_freq)
Dabar jau „aplaižytą“ dažnio signalą užkloju ant pirmojo, kad būtų matyti skirtumas, pritraukta truputį iš arčiau:
Gali atrodyti, kad šitas „aplaižymas“ nelabai reikalingas. Bet kai kur aplink centrinį dažnį yra mažų pasvyravimų, kurie dekoduojant gali trukdyti.
Tad čia yra demoduliuotas signalas. Kitas žingsnis – jo dekodavimas. Ten jau taip paprasta nebus, kaip su plačiomis Pitono bibliotekomis. Laikas aiškintis DSC žinutės sandarą. O ji štai tokia:
Pradžioje siunčiamas dot pattern, kuris matosi ir mano pateiktuose grafikuose. Tai tiesiog vienetukų/nuliukų krūva. Imtuvas šią seką gali pasigauti ir susikalibruoti pagal jų svyravimus. Nes siųstuvų ir imtuvų laikrodžiai nebėga visiškai idealiai vienodai, plius dar gali kiti faktoriai turėti įtakos, tad čia yra labai gera signalo dalis sutikrinti imtuvo laikrodžiui. Arba analizuojant WAV failą pasigauti vieno bituko žingsnį. Toliau eina speciali fazavimo seka iš DX ir RX simbolių, tada žinutės tipas ir jau toliau visokie duomenys. Iš viso 62 simboliai po 10 bitų, t.y. 620 bitų. Dot pattern gali būti arba 20, arba 200 bitų ilgio, priklausomai nuo žinutės tipo.
Uždavinys iš dalies paprastas: skenuojam šitą signalą ir fiksuojam perėjimus per centrinį dažnį bei tarpus tarp perėjimų. Viena problema, kurią radau, kad signalas yra „smailesnis“ ties viršutiniu dažniu, t.y. tarpas matuojant viršutinę bangelę per centrinį dažnį yra trumpesnis, nei matuojant apatinę. Šitas efektas gali sukliudyti pagauti pirmą DX simbolį, kuris, kaip matom grafike, turi didesnį tarpą po dot pattern. Todėl analizuojant dot pattern reikia pasigauti vidutinį bito žingsnį, taip pat skirtumus nuo vidutinio viršutinei bangelei ir apatinei bei apsiskaičiuoti žingsnio pločio korekciją. Daugoka ekstra skaičiavimų, bet nelabai yr kas daryt. Galima būtų adaptyviai kilnoti aukštyn/žemyn centrinį dažnį, kol bėga dot pattern ir rasti tokį, ties kuriuo kirtimų tarpai vienodi. Bet bandžiau: su prastesnės kokybės signalu kartais taip priartėjama prie apatinio dažnio, kad jo triukšmai pradeda rodytis, kaip netikri kirtimai. Tad irgi netinka, pločio korekcija yra tiesiog patikimiau.
Čia paveiksliukas iš arti pritrauktos bangelės, kad būtų aišku, jog tas žingsnio skirtumas yra. Galvoju, gal čia dar ir dažnis nuo centrinio nukrypęs, kad taip akivaizdžiai skirtumas matyt:
Ir štai tiek. Aš tokią mėgėjišką pradžią pasidariau, bet signalo dekodavimas gavosi labai nestabilus ir priklausomas nuo centrinio dažnio parinkimo mažiau, nei dešimties hercų tikslumu. Tai nėra geras rezultatas ir klaidų dekodavime lieka daug. Aš tuo kodo gabalu dalintis nesiruošiu, nes jis visiškai „ant snarglių“. Kaip ir minėjau, šita rašliava buvo daugiau pažintinė, dekoderį aš darysiu Fldigi kodo pagrindu, galbūt netgi ir integruosiu į patį Fldigi pilnai. Pirmas etapas bus dekodavimas iš WAV failo, o antras – galima integracija. Tad kol kas tiek.
Literatūra:
- DSC Wikipedijoje.
- DSC žinučių specifikacija (PDF).
- I ir Q signalų paaiškinimas, vienas iš suprantamiausių, ką radau.
- Baisiai primityvus RTTY demoduliavimas su Pitonu. Dekodavimo nėra.
- PySDR: filtrų kūrimas, stumdymas ir naudojimas.
- PySDR: amplitudės ir fazės nustatymas IQ sraute (su kompleksiniais skaičiais)
- Sinusų/kosinusų kvadratas.
- Hilberto analizė amplitudės ir fazės nustatymui.
- FSK demoduliavimas su Pitonu.
- Signalo apglostymas su Butterworth filtru.
- Labai daug visko apie NAVTEX (ta pati moduliacija) demoduliavimą ir dekodavimą.
- Swling: DSC stebėjimas, dažniai ir paaiškinimai žmonių kalba.
Darau, blė 