diff --git a/aerpaw/results/plotRadioLogs.m b/aerpaw/results/plotRadioLogs.m index 7371d73..c53ddbe 100644 --- a/aerpaw/results/plotRadioLogs.m +++ b/aerpaw/results/plotRadioLogs.m @@ -43,12 +43,44 @@ function [f, fDist, R] = plotRadioLogs(resultsPath, G, tLim) metricNames = ["SNR", "Power", "Quality", "PathLoss", "NoiseFloor", "FreqOffset"]; yLabels = ["SNR (dB)", "Power (dB)", "Quality", "Path Loss (dB)", "Noise Floor (dB)", "Freq Offset (MHz)"]; + nMetrics = numel(metricNames); % --- Time-based figure --- f = figure; - tl = tiledlayout(numel(metricNames), 1, 'TileSpacing', 'compact', 'Padding', 'compact'); + tl = tiledlayout(nMetrics + 1, 1, 'TileSpacing', 'compact', 'Padding', 'compact'); - for mi = 1:numel(metricNames) + % Distance vs time tile (first) + ax = nexttile(tl); + hold(ax, 'on'); grid(ax, 'on'); + legendEntries = string.empty; + ci = 1; + if ~isempty(G) + for rxIdx = 1:nUAV + tbl = R{rxIdx}; + txIDs = unique(tbl.TxUAVID); + for ti = 1:numel(txIDs) + txID = txIDs(ti); + rows = tbl(tbl.TxUAVID == txID, :); + rows = rows(rows.Timestamp >= tLim(1) & rows.Timestamp <= tLim(2), :); + if isempty(rows), continue; end + [~, ia] = unique(rows.Timestamp); + [radioPt, dist] = pairDist(rows(ia, :), G); + if isempty(dist) || all(isnan(dist)), continue; end + valid = ~isnan(dist); + si = mod(ci - 1, numel(styles)) + 1; + plot(ax, datetime(radioPt(valid), 'ConvertFrom', 'posixtime'), dist(valid), ... + styles(si), 'Color', colors(ci, :), 'MarkerSize', 3, 'LineWidth', 0.5); + legendEntries(end+1) = sprintf("TX %d → RX %d", txID, rows.RxUAVID(1)); %#ok + ci = ci + 1; + end + end + end + ylabel(ax, 'Distance (m)'); + xlabel(ax, 'Time'); + legend(ax, legendEntries, 'Location', 'best'); + hold(ax, 'off'); + + for mi = 1:nMetrics ax = nexttile(tl); hold(ax, 'on'); grid(ax, 'on'); @@ -63,21 +95,30 @@ function [f, fDist, R] = plotRadioLogs(resultsPath, G, tLim) rows = tbl(tbl.TxUAVID == txID, :); rows = rows(rows.Timestamp >= tLim(1) & rows.Timestamp <= tLim(2), :); vals = rows.(metricNames(mi)); + valid = ~isnan(vals); + rows = rows(valid, :); + vals = vals(valid); - if isempty(rows) || all(isnan(vals)) + if isempty(rows) continue; end si = mod(ci - 1, numel(styles)) + 1; plot(ax, rows.Timestamp, vals, styles(si), ... - 'Color', colors(ci, :), 'MarkerSize', 3, 'LineWidth', 1); + 'Color', colors(ci, :), 'MarkerSize', 3, 'LineWidth', 0.5); legendEntries(end+1) = sprintf("TX %d → RX %d", txID, tbl.RxUAVID(1)); %#ok + + % Median per 1/3-second time bin + [t_med, v_med] = timeBinMedian(posixtime(rows.Timestamp), vals, 1/3); + plot(ax, datetime(t_med, 'ConvertFrom', 'posixtime'), v_med, '-', ... + 'Color', 'r', 'LineWidth', 2); + legendEntries(end+1) = sprintf("TX %d → RX %d (median)", txID, tbl.RxUAVID(1)); %#ok ci = ci + 1; end end ylabel(ax, yLabels(mi)); - if mi == numel(metricNames) + if mi == nMetrics xlabel(ax, 'Time'); end legend(ax, legendEntries, 'Location', 'best'); @@ -93,9 +134,38 @@ function [f, fDist, R] = plotRadioLogs(resultsPath, G, tLim) return; end - tl2 = tiledlayout(numel(metricNames), 1, 'TileSpacing', 'compact', 'Padding', 'compact'); + tl2 = tiledlayout(nMetrics + 1, 1, 'TileSpacing', 'compact', 'Padding', 'compact'); - for mi = 1:numel(metricNames) + % Distance vs time tile (first) + ax = nexttile(tl2); + hold(ax, 'on'); grid(ax, 'on'); + legendEntries = string.empty; + ci = 1; + for rxIdx = 1:nUAV + tbl = R{rxIdx}; + txIDs = unique(tbl.TxUAVID); + for ti = 1:numel(txIDs) + txID = txIDs(ti); + rows = tbl(tbl.TxUAVID == txID, :); + rows = rows(rows.Timestamp >= tLim(1) & rows.Timestamp <= tLim(2), :); + if isempty(rows), continue; end + [~, ia] = unique(rows.Timestamp); + [radioPt, dist] = pairDist(rows(ia, :), G); + if isempty(dist) || all(isnan(dist)), continue; end + valid = ~isnan(dist); + si = mod(ci - 1, numel(styles)) + 1; + plot(ax, datetime(radioPt(valid), 'ConvertFrom', 'posixtime'), dist(valid), ... + styles(si), 'Color', colors(ci, :), 'MarkerSize', 3, 'LineWidth', 0.5); + legendEntries(end+1) = sprintf("TX %d → RX %d", txID, rows.RxUAVID(1)); %#ok + ci = ci + 1; + end + end + ylabel(ax, 'Distance (m)'); + xlabel(ax, 'Time'); + legend(ax, legendEntries, 'Location', 'best'); + hold(ax, 'off'); + + for mi = 1:nMetrics ax = nexttile(tl2); hold(ax, 'on'); grid(ax, 'on'); @@ -119,61 +189,39 @@ function [f, fDist, R] = plotRadioLogs(resultsPath, G, tLim) end vals = rows.(metricNames(mi)); - if all(isnan(vals)) + valid = ~isnan(vals); + rows = rows(valid, :); + vals = vals(valid); + + if isempty(rows) continue; end - % Map 0-based UAV IDs to 1-based GPS cell indices - txGpsIdx = double(txID) + 1; - rxGpsIdx = double(rows.RxUAVID(1)) + 1; + [radioPt, dist] = pairDist(rows, G); + if isempty(dist) || all(isnan(dist)), continue; end - if txGpsIdx > numel(G) || rxGpsIdx > numel(G) - continue; - end - - Gtx = G{txGpsIdx}; - Grx = G{rxGpsIdx}; - - if ~ismember('East', Gtx.Properties.VariableNames) || ... - ~ismember('East', Grx.Properties.VariableNames) - continue; - end - - % Strip timezone before posixtime so radio and GPS timestamps - % are treated on the same scale (both are AERPAW wall-clock time) - txTs = Gtx.Timestamp; txTs.TimeZone = ''; - rxTs = Grx.Timestamp; rxTs.TimeZone = ''; - txPt = posixtime(txTs); - rxPt = posixtime(rxTs); - radioPt = posixtime(rows.Timestamp); - - % Interpolate GPS positions at radio measurement times. - % Exclude NaN ENU entries (outside algorithm flight range). - validTx = ~isnan(Gtx.East); - validRx = ~isnan(Grx.East); - - txE = interp1(txPt(validTx), Gtx.East(validTx), radioPt, 'linear', NaN); - txN = interp1(txPt(validTx), Gtx.North(validTx), radioPt, 'linear', NaN); - txU = interp1(txPt(validTx), Gtx.Up(validTx), radioPt, 'linear', NaN); - rxE = interp1(rxPt(validRx), Grx.East(validRx), radioPt, 'linear', NaN); - rxN = interp1(rxPt(validRx), Grx.North(validRx), radioPt, 'linear', NaN); - rxU = interp1(rxPt(validRx), Grx.Up(validRx), radioPt, 'linear', NaN); - - dist = vecnorm([txE - rxE, txN - rxN, txU - rxU], 2, 2); - - if all(isnan(dist)) - continue; - end + % Drop points where GPS interpolation returned NaN + validDist = ~isnan(dist); + rowTs = radioPt(validDist); + dist = dist(validDist); + vals = vals(validDist); si = mod(ci - 1, numel(styles)) + 1; scatter(ax, dist, vals, 9, colors(ci, :), strrep(styles(si), "-", ""), 'filled'); legendEntries(end+1) = sprintf("TX %d → RX %d", txID, rows.RxUAVID(1)); %#ok + + % Median per 1/3-second time bin, plotted against median distance + [~, dv_med] = timeBinMedian(rowTs, [dist, vals], 1/3); + [d_med, si_sort] = sort(dv_med(:, 1)); + v_med = dv_med(si_sort, 2); + plot(ax, d_med, v_med, '-', 'Color', 'r', 'LineWidth', 2); + legendEntries(end+1) = sprintf("TX %d → RX %d (median)", txID, rows.RxUAVID(1)); %#ok ci = ci + 1; end end ylabel(ax, yLabels(mi)); - if mi == numel(metricNames) + if mi == nMetrics xlabel(ax, 'Distance (m)'); end legend(ax, legendEntries, 'Location', 'best'); @@ -182,3 +230,30 @@ function [f, fDist, R] = plotRadioLogs(resultsPath, G, tLim) title(tl2, 'Radio Channel Metrics vs Distance'); end + + +function [radioPt, dist] = pairDist(rows, G) + % Interpolate GPS-based inter-UAV distance at each row's timestamp. + radioPt = []; dist = []; + txGpsIdx = double(rows.TxUAVID(1)) + 1; + rxGpsIdx = double(rows.RxUAVID(1)) + 1; + if txGpsIdx > numel(G) || rxGpsIdx > numel(G), return; end + Gtx = G{txGpsIdx}; + Grx = G{rxGpsIdx}; + if ~ismember('East', Gtx.Properties.VariableNames) || ... + ~ismember('East', Grx.Properties.VariableNames), return; end + txTs = Gtx.Timestamp; txTs.TimeZone = ''; + rxTs = Grx.Timestamp; rxTs.TimeZone = ''; + txPt = posixtime(txTs); + rxPt = posixtime(rxTs); + radioPt = posixtime(rows.Timestamp); + validTx = ~isnan(Gtx.East); + validRx = ~isnan(Grx.East); + txE = interp1(txPt(validTx), Gtx.East(validTx), radioPt, 'linear', NaN); + txN = interp1(txPt(validTx), Gtx.North(validTx), radioPt, 'linear', NaN); + txU = interp1(txPt(validTx), Gtx.Up(validTx), radioPt, 'linear', NaN); + rxE = interp1(rxPt(validRx), Grx.East(validRx), radioPt, 'linear', NaN); + rxN = interp1(rxPt(validRx), Grx.North(validRx), radioPt, 'linear', NaN); + rxU = interp1(rxPt(validRx), Grx.Up(validRx), radioPt, 'linear', NaN); + dist = vecnorm([txE - rxE, txN - rxN, txU - rxU], 2, 2); +end diff --git a/aerpaw/results/readRadioLogs.m b/aerpaw/results/readRadioLogs.m index a54079d..07c9007 100644 --- a/aerpaw/results/readRadioLogs.m +++ b/aerpaw/results/readRadioLogs.m @@ -70,6 +70,40 @@ function R = readRadioLogs(logPath) R = sortrows(R, "Timestamp"); + % Reconstruct per-measurement timestamps within GNURadio processing batches. + % The flowgraph accumulates one full PN sequence (4095 chips at samp_rate/sps) + % per measurement, but outputs the whole batch simultaneously with a single + % wall-clock timestamp. We reassign timestamps by counting backward from the + % batch processing time at the known PN period interval. + pn_period = 4095 / (2e6 / 16); % 32.76 ms per PN correlation period + + for txId = unique(R.TxUAVID)' + rows = find(R.TxUAVID == txId); + if numel(rows) < 2, continue; end + + dt = seconds(diff(R.Timestamp(rows))); + break_pos = [1; find(dt > 0.5) + 1]; + end_pos = [break_pos(2:end) - 1; numel(rows)]; + + for b = 1:numel(break_pos) + idx = rows(break_pos(b) : end_pos(b)); + batch_ts = posixtime(R.Timestamp(idx)); + t_ref = max(batch_ts); + + % Multiple metric files share the same processing timestamp for + % each PN period, so group by unique original timestamp rather + % than treating every row as a separate PN period. + [~, ~, group_id] = unique(batch_ts); + n_groups = max(group_id); + new_ts = t_ref - (n_groups - 1 : -1 : 0)' * pn_period; + + for g = 1:n_groups + R.Timestamp(idx(group_id == g)) = ... + datetime(new_ts(g), 'ConvertFrom', 'posixtime'); + end + end + end + % Remove rows during defined guard period between TDM shifts R(R.TxUAVID == -1, :) = []; diff --git a/aerpaw/results/timeBinMedian.m b/aerpaw/results/timeBinMedian.m new file mode 100644 index 0000000..f2972fa --- /dev/null +++ b/aerpaw/results/timeBinMedian.m @@ -0,0 +1,29 @@ +function [t_med, v_med] = timeBinMedian(t, v, binWidth) + % Compute median of each column of v within fixed-width time bins. + % + % t - (N,1) posixtime values + % v - (N,K) data matrix; one column per quantity + % binWidth - scalar bin width in seconds + % + % t_med - (B,1) median time of each non-empty bin + % v_med - (B,K) median of each column per non-empty bin + + edges = (floor(min(t) / binWidth) * binWidth) : binWidth : ... + (floor(max(t) / binWidth) * binWidth + binWidth); + bins = discretize(t, edges); + nBins = numel(edges) - 1; + K = size(v, 2); + + t_all = NaN(nBins, 1); + v_all = NaN(nBins, K); + for bi = 1:nBins + mask = bins == bi; + if ~any(mask), continue; end + t_all(bi) = median(t(mask)); + v_all(bi,:) = median(v(mask,:), 1); + end + + ok = ~isnan(t_all); + t_med = t_all(ok); + v_med = v_all(ok, :); +end diff --git a/test/test_miSim.m b/test/test_miSim.m index 0d0572c..0a517d3 100644 --- a/test/test_miSim.m +++ b/test/test_miSim.m @@ -9,8 +9,8 @@ classdef test_miSim < matlab.unittest.TestCase plotCommsGeometry = false; % disable plotting communications geometries % Sim - maxIter = 50; - timestep = 0.05; + maxIter = 250; + timestep = 0.1; % Domain domain = rectangularPrism; % domain geometry @@ -31,7 +31,7 @@ classdef test_miSim < matlab.unittest.TestCase % Agents initialStepSize = 0.2; % gradient ascent step size at the first iteration. Decreases linearly to 0 based on maxIter. - initialMaxAngleStepSize = 5; % angular step size (degrees) for tilt/azimuth gradient ascent per timestep. + initialMaxAngleStepSize = 0.1; % angular step size (degrees) for tilt/azimuth gradient ascent per timestep. minAgents = 3; % Minimum number of agents to be randomly generated maxAgents = 4; % Maximum number of agents to be randomly generated useDoubleIntegrator = false;