% Skript vytvorený v rámci diplomovej práce Metody vyhodnocení kvality
% GNSS signálu za pomoci Android aplikací, pre účely vyhodnotenia SQM
% parametrov z nameraných GNSS dát z Android smartfónu
% 
% Autor: Bc. Simona Blašková, blasksim@cvut.cz

%% 
clear
clc

data1 = readlines("filename1.txt");    % namiesto filename 1 vložiť názov .txt súboru
data2 = readlines("filename2.nmea");   % namiesto filename 2 vložiť názov .nmea súboru

%% Spracovanie NMEA - parameter DOP

rawNMEAdata = erase(data2,'NMEA,'); % odstránenie NMEA zo začiatku riadku

% Vytvorenie objektov
pnmeaGGA=nmeaParser("MessageIDs","GGA");
pnmeaGSA=nmeaParser("MessageIDs","GSA");

c1=0; % pomocná premenná
for k=1:height(rawNMEAdata)
    [GGAisValid, ~] = extractNMEASentence(rawNMEAdata(k,1),'GGA'); % overenie validity GGA správy
    [GSAisValid, ~] = extractNMEASentence(rawNMEAdata(k,1),'GSA');
    if GGAisValid ==1
        c1=c1+1;
        [ggaData] = pnmeaGGA(rawNMEAdata(k,1));
        NMEA_table(c1,1) = {ggaData.MessageID};
        NMEA_table(c1,2) = {datetime(ggaData.UTCTime)};
        NMEA_table(c1,3) = {ggaData.NumSatellitesInUse};
        NMEA_table(c1,4) = {ggaData.Latitude};
        NMEA_table(c1,5) = {ggaData.Longitude};
        NMEA_table(c1,6) = {ggaData.Altitude};
        NMEA_table(c1,7) = {ggaData.QualityIndicator};
    elseif GSAisValid==1
        [gsaData] = pnmeaGSA(rawNMEAdata(k,1));
        NMEA_table(c1,1) = {ggaData.MessageID + gsaData.MessageID};
        NMEA_table(c1,8) = {gsaData.PDOP};
        NMEA_table(c1,9) = {gsaData.VDOP};
        NMEA_table(c1,10) = {gsaData.HDOP};
    else
    end
end

NMEA_table = cell2table(NMEA_table);
NMEA_table.Properties.VariableNames = ["MessageID","UTCTime","NumSatInUse","Latitude","Longitude","Altitude","QualityIndicator","PDOP","VDOP","HDOP"];
NMEA_table.UTCTime.Format = 'hh:mm:ss';

clear c1 filename2 k ggaData GGAisValid gsaData GSAisValid pnmeaGSA pnmeaGGA rawNMEAdata rawNMEAdata_full sz varNames varTypes

idx = NMEA_table.QualityIndicator == 0; % vymazanie riadkov s invalid FIX
NMEA_table(idx,:) = [];

% Vyfiltrovanie požadovaného časového úseku (iba UTC čas)
starts = datetime('30-Nov--0001 hh:mm:ss',"Format","default","TimeZone","UTC");     % vloženie počiatku časového úseku hh:mm:ss
ends = datetime('30-Nov--0001 hh:mm:ss',"Format","default","TimeZone","UTC");       % vloženie konca časového úseku hh:mm:ss
idx = isbetween(NMEA_table.UTCTime, starts, ends);
NMEA_table = NMEA_table(idx,:);
NMEA_table.UTCTime.Format = "hh:mm:ss";

clear starts ends idx

% Vykreslenie PDOP, HDOP, VDOP
figure(1)
plot(NMEA_table.UTCTime, NMEA_table.PDOP)
hold on
plot(NMEA_table.UTCTime, NMEA_table.HDOP)
plot(NMEA_table.UTCTime, NMEA_table.VDOP)
title("DOP")
xlabel("Čas (UTC)")
ylabel("DOP")
legend(["PDOP" "HDOP" "VDOP"],"Location","eastoutside")
xlim([NMEA_table.UTCTime(1) NMEA_table.UTCTime(end)])
hold off

% Štatistika PDOP
avgPDOP = mean(NMEA_table.PDOP);
stdPDOP = std(NMEA_table.PDOP);
medPDOP = median(NMEA_table.PDOP);

%% Vyfiltrovanie riadkov #Raw
raw_num = contains(data1,"Raw");
data_raw = data1(raw_num,:);
data_raw_split = split(data_raw,',');
header = data_raw_split(1,:);
data_raw_split(1,:) = [];
raw = splitvars(table(data_raw_split),"data_raw_split","NewVariableNames",header);
raw = convertvars(raw,{'AccumulatedDeltaRangeMeters' 'AccumulatedDeltaRangeUncertaintyMeters' 'AgcDb' 'AccumulatedDeltaRangeState' 'BasebandCn0DbHz' 'BiasNanos' 'BiasUncertaintyNanos' 'CarrierCycles' 'CarrierFrequencyHz' 'CarrierPhase' 'CarrierPhaseUncertainty' 'ChipsetElapsedRealtimeNanos' 'Cn0DbHz' 'DriftNanosPerSecond' 'DriftUncertaintyNanosPerSecond' 'FullBiasNanos' 'FullInterSignalBiasNanos' 'FullInterSignalBiasUncertaintyNanos' 'HardwareClockDiscontinuityCount' 'LeapSecond' 'MultipathIndicator' 'PseudorangeRateMetersPerSecond' 'PseudorangeRateUncertaintyMetersPerSecond' 'ReceivedSvTimeNanos' 'ReceivedSvTimeUncertaintyNanos' 'SatelliteInterSignalBiasNanos' 'SatelliteInterSignalBiasUncertaintyNanos' 'SnrInDb' 'State' 'Svid' 'TimeNanos' 'TimeOffsetNanos' 'TimeUncertaintyNanos' 'utcTimeMillis'},"double");
raw.utcTimeMillis = datetime(1970,1,1,'Format','dd-MMM-yyyy HH:mm:ss') + milliseconds(raw.utcTimeMillis);
raw = renamevars(raw,"utcTimeMillis","UTCTime");
clear data_raw_split raw_num data_raw

% Vyfiltrovanie požadovaného časového úseku (čas zadávať v UTC)
starts = datetime(20yy,mm,dd,hh,mm,ss);     % vloženie počiatku časového úseku (year, month, day, hour, minutes, seconds)
ends = datetime(20yy,mm,dd,hh,mm,ss);       % vloženie konca časového úseku (year, month, day, hour, minutes, seconds)
idx = isbetween(raw.UTCTime, starts, ends);
raw = raw(idx,:);
clear idx header starts ends

%% Vyfiltrovanie riadkov #Fix
fix_num = contains(data1,"Fix");
data_fix = data1(fix_num,:);
data_fix(2,:) = []; % odstránenie 2. riadku - string "Fix" sa nachádza aj v riadkoch označených #Status
data_fix = split(data_fix,',');
header = data_fix(1,:);
data_fix(1,:) = [];
fix = splitvars(table(data_fix),"data_fix","NewVariableNames",header);
fix = convertvars(fix,{'LatitudeDegrees' 'LongitudeDegrees' 'AccuracyMeters' 'AltitudeMeters' 'BearingAccuracyDegrees' 'BearingDegrees' 'elapsedRealtimeNanos' 'SpeedAccuracyMps' 'SpeedMps' 'VerticalAccuracyMeters' 'UnixTimeMillis'},"double");
fix.UnixTimeMillis = datetime(1970,1,1,'Format','dd-MMM-yyyy HH:mm:ss') + milliseconds(fix.UnixTimeMillis);
fix = renamevars(fix,"UnixTimeMillis","UTCTime");
fix(:,1) = [];
clear header fix_num data_fix

% Vyfiltrovanie požadovaného časového úseku
starts = datetime(20yy,mm,dd,hh,mm,ss);     % vloženie počiatku časového úseku (year, month, day, hour, minutes, seconds)
ends = datetime(20yy,mm,dd,hh,mm,ss);       % vloženie konca časového úseku (year, month, day, hour, minutes, seconds)
idx = isbetween(fix.UTCTime, starts, ends);
fix = fix(idx,:);
clear idx header starts ends

%% Časový vektor
t1 = raw.UTCTime(1);
t2 = raw.UTCTime(end);
t = (t1:seconds(1):t2)';
clear t1 t2

%% Vykreslenie vypočítanej polohy z #Fix a výpočet chyby polohy (ak je známa referenčná poloha)

% Referenčné súradnice (ak sú k dispozícii)
lat_ref = ;      % súradnice zemepisnej šírky v °
lon_ref = ;      % súradnice zemepisnej výšky v °
h_ref = ;        % výška v metroch

% Namerané súradnice
lat = fix.LatitudeDegrees;
lon = fix.LongitudeDegrees;
h = fix.AltitudeMeters;

figure(2)
geoplot(lat,lon,".",'Color','blue','MarkerSize', 10,'DisplayName','Surové GNSS dáta zo smartfónu')
hold on
geoplot(lat_ref,lon_ref,"*",'MarkerSize',10,'Color','red','DisplayName','Referenčná poloha')
geolimits([lat_ref-0.00010,lat_ref+0.00010],[lon_ref-0.00010,lon_ref+0.00010])
geobasemap topographic
legend()
hold off

% Výpočet odchýlok
[E,N,Up] = geodetic2enu(lat, lon, h, lat_ref, lon_ref, h_ref, wgs84Ellipsoid);

% RMSE
% RMSE hodnoty pre každú observáciu
for i = 1:length(N)
    RMS_3D(i) = sqrt((N(i))^2+(E(i))^2+(Up(i))^2);
end
for i = 1:length(N)
    RMS_horiz(i) = sqrt((N(i))^2+(E(i))^2);
end
clear i
RMS_vert = abs(Up);

% RMSE hodnoty pre všetky observácie
RMS_3D_all = sqrt(sum((N).^2+(E).^2+(Up).^2)/length(N));
RMS_horiz_all = sqrt(sum((N).^2+(E).^2)/length(N));
RMS_vert_all = sqrt(sum((Up).^2)/length(Up));

% Grafy
% Chyba polohy v čase
figure(3)
plot(t, N,'DisplayName','N','Color',"#0072BD");
hold on
plot(t, E,'DisplayName','E','Color',"#D95319")
plot(t, Up,'DisplayName','Up','Color',"#EDB120")
title("Chyba určenia polohy v čase")
ylabel("Chyba polohy (m)")
xlabel("Čas (UTC)")
legend
hold off

% RMSE
figure(4)
subplot(3,1,1)
plot(t,RMS_3D,'.','Color',"#A2142F")
ylim([0 40])
title("3D RMSE")
subplot(3,1,2)
plot(t,RMS_horiz,'.','Color',"#77AC30")
ylim([0 40])
title("Horizontal RMSE")
subplot(3,1,3)
plot(t,RMS_vert,'.','Color',"#7E2F8E")
ylim([0 40])
xlabel("Time")
title("Vertical RMSE")

%% Vyfiltrovanie dát z jednotlivých GNSS konštelácií
% ConstellationType:
% GPS = 1
% SBAS = 2
% GLONASS = 3
% QZSS = 4
% BeiDou = 5
% Galileo = 6
% IRNSS = 7
% unknown = 0

% GPS
gps_rows = contains(raw.ConstellationType,"1");
gps = raw(gps_rows,:);
gps_L1 = gps(contains(gps.CodeType,"C"),:);              % L1
gps_L5 = gps(contains(gps.CodeType,"I"),:);              % L5

% Galileo
galileo_rows = contains(raw.ConstellationType,"6");
galileo = raw(galileo_rows,:);
galileo_E1 = galileo(contains(galileo.CodeType,"C"),:);  % E1
galileo_E5 = galileo(contains(galileo.CodeType,"I"),:);  % E5

% GLONASS
glonass_rows = contains(raw.ConstellationType,"3");
glonass = raw(glonass_rows,:);

% BeiDou
beidou_rows = contains(raw.ConstellationType,"5");
beidou = raw(beidou_rows,:);

clear gps_rows glonass_rows beidou_rows galileo_rows

%% Priemerné C/No v čase

% GPS
gps_table = table(gps.UTCTime, gps.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
gps_table = table2timetable(gps_table);
avg_CNo_gps = retime(gps_table,"secondly","mean");
avg_CNo_gps = avg_CNo_gps.CNo;

gps_L1_table = table(gps_L1.UTCTime, gps_L1.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
gps_L1_table = table2timetable(gps_L1_table);
avg_CNo_gps_L1 = retime(gps_L1_table,"secondly","mean");
avg_CNo_gps_L1 = avg_CNo_gps_L1.CNo;

gps_L5_table = table(gps_L5.UTCTime, gps_L5.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
gps_L5_table = table2timetable(gps_L5_table);
avg_CNo_gps_L5 = retime(gps_L5_table,"secondly","mean");
avg_CNo_gps_L5 = avg_CNo_gps_L5.CNo;
clear gps_table gps_L1_table gps_L5_table

% Galileo
galileo_table = table(galileo.UTCTime, galileo.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
galileo_table = table2timetable(galileo_table);
avg_CNo_galileo = retime(galileo_table,"secondly","mean");
avg_CNo_galileo = avg_CNo_galileo.CNo;

galileo_E1_table = table(galileo_E1.UTCTime, galileo_E1.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
galileo_E1_table = table2timetable(galileo_E1_table);
avg_CNo_galileo_E1 = retime(galileo_E1_table,"secondly","mean");
avg_CNo_galileo_E1 = avg_CNo_galileo_E1.CNo;

galileo_E5_table = table(galileo_E5.UTCTime, galileo_E5.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
galileo_E5_table = table2timetable(galileo_E5_table);
avg_CNo_galileo_E5 = retime(galileo_E5_table,"secondly","mean");
avg_CNo_galileo_E5 = avg_CNo_galileo_E5.CNo;
clear galileo_table galileo_L1_table galileo_L5_table

% GLONASS
glonass_table = table(glonass.UTCTime, glonass.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
glonass_table = table2timetable(glonass_table);
avg_CNo_glonass = retime(glonass_table,"secondly","mean");
avg_CNo_glonass = avg_CNo_glonass.CNo;
clear glonass_table

% BeiDou
beidou_table = table(beidou.UTCTime, beidou.Cn0DbHz,'VariableNames',{'UTCTime' 'CNo'});
beidou_table = table2timetable(beidou_table);
avg_CNo_beidou = retime(beidou_table,"secondly","mean");
avg_CNo_beidou = avg_CNo_beidou.CNo;
clear beidou_table

% Všetky GNSS konštelácie
G = findgroups(raw.UTCTime);
avg_CNo_all = splitapply(@(x) mean(x), raw.Cn0DbHz, G);

% Grafy
figure(5)
plot(t, avg_CNo_gps_L1,"Color","#0072BD")
hold on
plot(t, avg_CNo_galileo_E1,"Color","#77AC30")
plot(t, avg_CNo_glonass,"Color","#EDB120")
plot(t, avg_CNo_beidou,"Color","#7E2F8E")
plot(t, avg_CNo_all,"Color","#FF0000")
xlabel("Čas (UTC)")
ylabel("C/No (dB-Hz)")
xlim([t(1) t(end)])
title("Priemerné C/No v čase (L1/E1/G1/B1)")
legend({'GPS' 'Galileo' 'GLONASS' 'Beidou' 'Celkom (L1+L5)'}, "Location", "southoutside")
hold off

figure(6)
plot(t, avg_CNo_gps_L5,"Color","#0072BD")
hold on
plot(t, avg_CNo_galileo_E5,"Color","#77AC30")
xlabel("Čas (UTC)")
ylabel("C/No (dB-Hz)")
xlim([t(1) t(end)])
title("Priemerné C/No v čase (L5/E5)")
legend({'GPS' 'Galileo'}, "Location", "southoutside")
hold off

% Priemerné hodnoty C/No za celé meranie
CNo_all = mean(raw.Cn0DbHz);
CNo_gps_L1 = mean(gps_L1.Cn0DbHz);
CNo_gps_L5 = mean(gps_L5.Cn0DbHz);
CNo_galileo_E1 = mean(galileo_E1.Cn0DbHz);
CNo_galileo_E5 = mean(galileo_E5.Cn0DbHz);
CNo_glonass = mean(glonass.Cn0DbHz);
CNo_beidou = mean(beidou.Cn0DbHz);

%% Priemerné C/No jednotlivých satelitov za celú dobu merania
unique_gps_L1 = unique(gps_L1.Svid);
unique_gps_L5 = unique(gps_L5.Svid);
unique_galileo_E1 = unique(galileo_E1.Svid);
unique_galileo_E5 = unique(galileo_E5.Svid);
unique_glonass = unique(glonass.Svid);
unique_beidou = unique(beidou.Svid);

% GPS
j = 1;
prns_gps = zeros(length(unique_gps_L1),1);              % L1
for i = min(unique_gps_L1):max(unique_gps_L1)
    if unique_gps_L1 ~= i
        i = i+1;
    else
        ix = find(gps_L1.Svid == i);
        prn = gps_L1(ix,:);
        prns_gps_L1(j,1) = mean(prn.Cn0DbHz);
        j = j+1;
    end
end

prns_gps = zeros(length(unique_gps_L5),1);              % L5
j = 1;
for i = min(unique_gps_L5):max(unique_gps_L5)
    if unique_gps_L5 ~= i
        i = i+1;
    else
        ix = find(gps_L5.Svid == i);
        prn = gps_L5(ix,:);
        prns_gps_L5(j,1) = mean(prn.Cn0DbHz);
        j = j+1;
    end
end

% Galileo
j = 1;
prns_galileo = zeros(length(unique_galileo_E1),1);       % E1
for i = min(unique_galileo_E1):max(unique_galileo_E1)
    if unique_galileo_E1 ~= i
        i = i+1;
    else
        ix = find(galileo_E1.Svid == i);
        prn = galileo_E1(ix,:);
        prns_galileo_E1(j,1) = mean(prn.Cn0DbHz);
        j = j+1;
    end
end

j = 1;
prns_galileo = zeros(length(unique_galileo_E5),1);       % E5
for i = min(unique_galileo_E5):max(unique_galileo_E5)
    if unique_galileo_E5 ~= i
        i = i+1;
    else
        ix = find(galileo_E5.Svid == i);
        prn = galileo_E5(ix,:);
        prns_galileo_E5(j,1) = mean(prn.Cn0DbHz);
        j = j+1;
    end
end

% GLONASS
j = 1;
prns_glonass = zeros(length(unique_glonass),1);
for i = min(unique_glonass):max(unique_glonass)
    if unique_glonass ~= i
        i = i+1;
    else
        ix = find(glonass.Svid == i);
        prn = glonass(ix,:);
        prns_glonass(j,1) = mean(prn.Cn0DbHz);
        j = j+1;
    end
end

% BeiDou
j = 1;
prns_beidou = zeros(length(unique_beidou),1);
for i = min(unique_beidou):max(unique_beidou)
    if unique_beidou ~= i
        i = i+1;
    else
        ix = find(beidou.Svid == i);
        prn = beidou(ix,:);
        prns_beidou(j,1) = mean(prn.Cn0DbHz);
        j = j+1;
    end
end

figure(7)
subplot(3,2,1)
bar(prns_gps_L1,'FaceColor',"#0072BD");
labels = string(unique_gps_L1);
set(gca,'XTickLabel',labels)
title("GPS L1")
ylim([0 50])
xlabel("PRN")
ylabel("C/No (dB-Hz)")

subplot(3,2,2)
bar(prns_gps_L5,'FaceColor',"#0072BD");
labels = string(unique_gps_L5);
set(gca,'XTickLabel',labels)
ylim([0 50])
title("GPS L5")
xlabel("PRN")
ylabel("C/No (dB-Hz)")

subplot(3,2,3)
bar(prns_galileo_E1,'FaceColor',"#77AC30");
title("Galileo E1")
labels = string(unique_galileo_E1);
set(gca,'XTickLabel',labels)
ylim([0 50])
xlabel("PRN")
ylabel("C/No (dB-Hz)")

subplot(3,2,4)
bar(prns_galileo_E5,'FaceColor',"#77AC30");
title("Galileo E5")
labels = string(unique_galileo_E5);
set(gca,'XTickLabel',labels)
ylim([0 50])
xlabel("PRN")
ylabel("C/No (dB-Hz)")

subplot(3,2,5)
bar(prns_glonass,'FaceColor',"#EDB120");
title("GLONASS G1")
labels = string(unique_glonass);
set(gca,'XTickLabel',labels)
ylim([0 50])
xlabel("PRN")
ylabel("C/No (dB-Hz)")

subplot(3,2,6)
bar(prns_beidou,'FaceColor',"#7E2F8E");
title("BeiDou B1")
labels = string(unique_beidou);
set(gca,'XTick',1:unique_beidou(end),'XTickLabel',labels)
ylim([0 50])
xlabel("PRN")
ylabel("C/No (dB-Hz)")

clear unique_beidou unique_galileo_L1 unique_galileo_L5 unique_glonass unique_gps_L1 unique_gps_L5 labels i j ix prn prns_galileo prns_gps prns_beidou prns_glonass prns_galileo_L1 prns_galileo_L5 prns_gps_L1 prns_gps_L5

%% Počet sledovaných satelitov

% GPS
G = findgroups(gps.UTCTime);
num_sats_gps = splitapply(@(x) numel(unique(x)),gps.Svid,G);
G = findgroups(gps_L1.UTCTime);
num_sats_gps_L1 = splitapply(@(x) numel(unique(x)),gps_L1.Svid,G);
G = findgroups(gps_L5.UTCTime);
num_sats_gps_L5 = splitapply(@(x) numel(unique(x)),gps_L5.Svid,G);

% Galileo
G = findgroups(galileo.UTCTime);
num_sats_galileo = splitapply(@(x) numel(unique(x)),galileo.Svid,G);
G = findgroups(galileo_E1.UTCTime);
num_sats_galileo_E1 = splitapply(@(x) numel(unique(x)),galileo_E1.Svid,G);
G = findgroups(galileo_E5.UTCTime);
num_sats_galileo_E5 = splitapply(@(x) numel(unique(x)),galileo_E5.Svid,G);

% GLONASS
G = findgroups(glonass.UTCTime);
num_sats_glonass = splitapply(@(x) numel(unique(x)),glonass.Svid,G);

% BeiDou
G = findgroups(beidou.UTCTime);
num_sats_beidou = splitapply(@(x) numel(unique(x)),beidou.Svid,G);

% Všetky GNSS konštelácie
num_sats_all = num_sats_gps + num_sats_galileo + num_sats_glonass + num_sats_beidou;

% Grafy
figure(8)
plot(t,num_sats_all,"Color","#FF0000")
hold on
plot(t,num_sats_gps_L1,"Color","#0072BD")
plot(t,num_sats_galileo_E1,"Color","#77AC30")
plot(t,num_sats_glonass,"Color","#EDB120")
plot(t,num_sats_beidou,"Color","#7E2F8E")
ylabel('NSV')
xlabel('Čas (UTC)')
xlim([t(1) t(end)])
legend({'Celkom (L1+L5)' 'GPS' 'Galileo' 'GLONASS' 'BeiDou'},"Location","southoutside")
title("Počet sledovaných satelitov (L1/E1/G1/B1)")
hold off

figure(9)
plot(t,num_sats_gps_L5,"Color","#0072BD")
hold on
plot(t,num_sats_galileo_E5,"Color","#77AC30")
ylabel('NSV')
xlabel('Čas (UTC)')
xlim([t(1) t(end)])
legend({'GPS' 'Galileo'},"Location","southoutside")
title("Počet sledovaných satelitov (L5/E5)")
hold off
clear i t1 t2 idx

% Priemerné hodnoty NSV za celé meranie
avg_nsv_all = mean(num_sats_all);
avg_nsv_gps_L1 = mean(num_sats_gps_L1);
avg_nsv_gps_L5 = mean(num_sats_gps_L5);
avg_nsv_galileo_E1 = mean(num_sats_galileo_E1);
avg_nsv_galileo_E5 = mean(num_sats_galileo_E5);
avg_nsv_glonass = mean(num_sats_glonass);
avg_nsv_beidou = mean(num_sats_beidou);

%% AGC
agc = splitapply(@(x) mean(x), raw.AgcDb, findgroups(discretize(raw.UTCTime, t)));

figure(10)
plot(t(1:length(agc)),agc);
ylim([min(agc)-1 max(agc)+1])
xlim([min(t) max(t)])
title("AGC")
ylabel('AGC (dB)')
xlabel('Čas (UTC)')

% Priemerná hodnota AGC za celé meranie
avg_agc = mean(agc);

%% Skyplot

% Vyfiltrovanie riadkov #Status
status_num = contains(data1,"Status");
data_status = data1(status_num,:);
data_status(1,:) = [];
status_split = split(data_status,',');
header = status_split(1,:);
status_split(1,:) = [];
status = splitvars(table(status_split),"status_split","NewVariableNames",header);
status = convertvars(status,{'AzimuthDegrees' 'BasebandCn0DbHz' 'CarrierFrequencyHz' 'Cn0DbHz' 'ElevationDegrees' 'HasAlmanacData' 'HasEphemerisData' 'SignalCount' 'SignalIndex' 'Svid' 'UsedInFix' 'UnixTimeMillis'},"double");
idx = status.UnixTimeMillis == 0;
status(idx,:) = [];
status.UnixTimeMillis = datetime(1970,1,1,'Format','dd-MMM-yyyy HH:mm:ss') + milliseconds(status.UnixTimeMillis);
status = renamevars(status,"UnixTimeMillis","UTCTime");
clear idx status_split header status_num

% Vykreslenie skyplot grafu
status(status.ElevationDegrees < 0,:) = [];
azimuth = status.AzimuthDegrees;
elevation = status.ElevationDegrees;
groups = status.ConstellationType;
figure(11)
skyplot(azimuth, elevation, "GroupData",categorical(groups),"MarkerSizeData",5)
legend('GPS', 'SBAS', 'GLONASS', 'BeiDou', 'Galileo')
title('Skyplot')

