%% Working with the Radar HDF5 Files
% These HDF5 files were generated from wireshark network captures of the radar data, and this file will show an example analysis of a single file, exploring the data contained within.

close all; clear;

% Table of Contents
% Goals
% HDF5 Radar Files
%     Loading parts of the HDF5 File
%     AHRS Data
%     Radar Reports
% Scan Display
%     Next Steps
%% Goals
% I have developed a program in C that converts the network capture files into archival, self describing data files using the HDF5 format. These files are heirarchal, and is a format designed for large dataset analyses like the ones we want to acheive with the radar data. HDF5 files can be ingested by any software toolset that has supported libraries, most notably MATLAB, Python, C and Fortran. If you are arriving at this document without being at the repository, you can find the conversion tools, MATLAB and python scripts (including this one) here:
% https://bitbucket.org/mbari/wec-pcap-utils/src/master/
%% HDF5 Radar Files
% Assuming you are in the repo, let's select a file relative to this script. You may have to change this for non-posix.

hdf_file = '../data/mbariAtSeaWecRadar_00315_20210928171518.h5';

% Let's use matlab's tools to print out the information in this file

%h5disp(hdf_file,'/','min')
h5disp(hdf_file,'/')

%% Loading parts of the HDF5 File
% You can learn more about the HDF file format at https://portal.hdfgroup.org/display/HDF5/HDF5
% One thing that is great is that you only need to load the items into memory that you want to work with, so we'll step into loading something familiar, the AHRS data.
%% AHRS Data
% Load the ahrs data into memory from the hdf5 file.

ahrs.t = h5read(hdf_file, '/ahrs/timestamp');
ahrs.r = h5read(hdf_file, '/ahrs/rollAngle');
ahrs.p = h5read(hdf_file, '/ahrs/pitchAngle');
ahrs.y = h5read(hdf_file, '/ahrs/yawAngleTrue');

% Let's convert the timestamp into something more plottable.
ahrs.dt = datetime(ahrs.t, 'ConvertFrom','epochtime');

% Plot this 5 minute window of time.

figure();
subplot(311);
plot(ahrs.dt,ahrs.r);
title(hdf_file,'Interpreter','none');
ylabel('roll');
subplot(312);
plot(ahrs.dt,ahrs.p);
ylabel('pitch');
subplot(313);
plot(ahrs.dt,ahrs.y);
ylabel('yaw');

% There is obviously more data to explore in the AHRS records, but mostly we need the yawAngleTrue to offset the plots of the incoming radar data. We get our first indication in this file that heading is ±180˚ instead of the typical 0-360˚ range.
% Also, I have started loading functions for the three main types, let's test the ahrs load function:

ahrs = load_radar_ahrs(hdf_file);
head(ahrs)

%% Radar Reports
% The reports channel tells us at 1hz the operational state of the radar. For the most part these are changing very little. This is our first example of of a loading function to simplify the import process.

reports_a = load_radar_reports(hdf_file,'a');

% Let's summarize the data
head(reports_a)
headers = reports_a.Properties.VariableNames;

for i = headers
    rpt_mean = mean(double(reports_a{:,i}));
    rpt_std = std(double(reports_a{:,i}));
    disp(sprintf("Mean: %10.1f\tStd: %5.2f\t %s", rpt_mean, rpt_std, i{1}));
end

% I'm still working on decoding these, for now we will move forward.
%% Scan Display
% Load in scan data and plot one full rotation. Its worth describing the nomenclature in the variables. The radar collects data in a polar coordinate system, for each "scan", there is a primary angle and a total range that return intensity is sampled 1024 times along. Hence, each row of the datatype as compiled in the table below contains one scan. Each full revolution of the radar is 2048 scans, or spokes.

scans_a = load_radar_scans(hdf_file,'a');
head(scans_a)
% plot 1 revolution, each scan is 2048 lines
% idx = 1:2048;
% single_revolution = scans_a(idx,:);
% 
% rho = 1:1024; %create a range vector for plotting, units ignored for now.
% theta = ones(1,1024); % a ones vector we will change each loop through scan data.

%

%now build a plot one scan at a time, let's plot 1 scan
% figure();
% for i = idx
%     ang = deg2rad(single_revolution.angle(i))*theta;
%     col = single_revolution.scan(i,:);
%     polarscatter(ang,rho, 1,col,'filled');
%     hold on;
% end

% Now processing the whole file doesn't make a ton of sense, that plot above took 30 seconds to make, so it isn't computationally efficient whatsoever. This file has about 180 full revolutions, but it does illustrate how to get a first pass set of images together.
%% Next Steps
% The converter tool doesn't work on every file, having some issue when data from different scans overlap during packet reconstruction, I'm working on this.
% The converter right now is very efficient, a day of data can be converted in < 10min. Once I polish out the above issue, we can process the entire dataset into HDF5 files that can be worked on inside matlab/python, etc.
% Verify data conversion, interpretation with Dale.
% Validate the range value conversion in the reports.
% Start analying data with ahrs in mind and generate short animations to help verify we have data in regions Rob can verify.