Writing a GSF-reader in Matlab

I recently ran across a couple issues with some GSF (Generic Sensor Format) files I have, so I decided to attempt to write a GSF reader in order to see what, exactly, I was dealing with. I have never written any kind of binary format reader before (though I did tinker around and add features to one written in Python), so I sort of jumped in the deep-end with this. I decided to use Matlab, since I am most comfortable with its scripting language at this point in time.

Here is a valuable lesson I just learned: Never take the binary specification file at face-value. Typos happen, and they certainly happened here. I have the GSF v. 3.01 specification file available from the SAIC website, and have been following it to write my reader. Everything has been going fine for the most part until I reached the attitude record.

According to the specification, the attitude record should look like this:
-->

Field Name	Description	Field Type	Count
NUM_MEASUREMENTS	Number of attitude measurements in this record.	I	2
ATTITUDE_TIME	Array of attitude measurement times	I	N*2*4
PITCH	Array of pitch measurements	I	N*4
ROLL	Array of roll measurements	I	N*4
HEAVE	Array of heave measurements	T	N*4
HEADING	Array of heading measurements	T	N*4
Attitude Record Size:			Variable

Should be pretty straightforward, so in my Matlab code I wrote the following function:

function [Num_meas, ATT_Time, ATT_Pitch, ATT_Roll, ATT_Heave, ATT_Heading] = readATTrecord(fid)
%% This function reads the Attitude record in a GSF file.

%% Element          Bytes   Type         Description
% Num_Measurements    2      int      number of attitude measurements (N)
% Attitude_Time       N*2*4  int      array of attitude meas. times
% Pitch               N*4    int      array of pitch measurements (hundreths of deg)
% Roll                N*4    int      array of roll measurements (hundreths of deg)
% Heave               N*4    char      array of heave measurements (hundreths of deg)
% Heading             N*4    char      array of heading measurements (hundreths of deg)

Num_meas = fread(fid,1,'int16');

for i = 1:Num_meas
ATT_Time(i,:) = fread(fid,2,'int32');
end
for i = 1:Num_meas
ATT_Pitch(i,:) = fread(fid,1,'int32');
end
for i = 1:Num_meas
ATT_Roll(i,:) = fread(fid,1,'int32');
end
for i = 1:Num_meas
ATT_Heave(i,:) = fread(fid,4,'char');
end
for i = 1:Num_meas
ATT_Heading(i,:) = fread(fid,4,'char');
end

This returned all sorts of funky data that made no sense, most notably the negative times and pitch angles of 120 degrees (that would be worse than the perfect storm!). I decided to look at the GSF library also available on the SAIC website in order to check out the source code. In the library directory, there is a file called GSF_dec.c, a C-code source file for decoding the GSF binary format. I looked up how the Attitude record was decoded in C, and saw this:

gsfDecodeAttitude(gsfAttitude *attitude, GSF_FILE_TABLE *ft, unsigned char *sptr)
{
unsigned char  *p = sptr;
gsfuLong        ltemp;
gsfuShort       stemp;
gsfsShort       signed_short;
int             i;
struct timespec basetime;
double          time_offset;

/* First four byte integer contains the observation time seconds */
memcpy(&ltemp, p, 4);
p += 4;
basetime.tv_sec = ntohl(ltemp);

/* Next four byte integer contains the observation time nanoseconds */
memcpy(&ltemp, p, 4);
p += 4;
basetime.tv_nsec = ntohl(ltemp);

/* Next two byte integer contains the number of measurements in the record */
memcpy(&stemp, p, 2);
p += 2;
attitude->num_measurements = ntohs(stemp);

>and a little farther down I saw:

/* Now loop to decode the attitude measurements */
for (i = 0; i <>num_measurements; i++)
{
/* Next two byte integer contains the time offset */
memcpy(&stemp, p, 2);
time_offset = ((double) ntohs (stemp)) / 1000.0;
p += 2;

LocalAddTimes (&basetime, time_offset, &attitude->attitude_time[i]);

/* Next two byte integer contains the pitch */
memcpy(&stemp, p, 2);
signed_short = (signed) ntohs(stemp);
attitude->pitch[i] = ((double) signed_short) / 100.0;
p += 2;

/* Next two byte integer contains the roll */
memcpy(&stemp, p, 2);
signed_short = (signed) ntohs(stemp);
attitude->roll[i] = ((double) signed_short) / 100.0;
p += 2;

/* Next two byte integer contains the heave */
memcpy(&stemp, p, 2);
signed_short = (signed) ntohs(stemp);
attitude->heave[i] = ((double) signed_short) / 100.0;
p += 2;

/* Next two byte integer contains the heading */
memcpy(&stemp, p, 2);
attitude->heading[i] = ((double) ntohs(stemp)) / 100.0;
p += 2;

The C-code shows that the specification table was completely wrong! There is only one reference time for the attitude reading, and the rest of the times are simply how many seconds past that reference time the rest of the measurements take place. Furthermore, the specifications stated that the Heave and Heading fields were text (something I thought was weird anyway) and in the C-code we clearly see they are integers. The specification also stated everything was 4-bytes, but it is actually only 2 (16 bits). This is very frustrating that the specification and the C-code do not match, especially since they both come from the same source.

At any rate, I am now able to write a working function that can be called in my main code:

function [ATT_Time, Num_meas, ATT_Offset, ATT_Pitch, ATT_Roll, ATT_Heave, ATT_Heading] = readATTrecord(fid)
%% This function reads the Attitude record in a GSF file.

%% Element          Bytes   Type         Description
% ATT_Time            2*4    int      time of attitude meas.
% Num_Measurements    2      int      number of attitude measurements (N)
% Time_Offset         N*2    int      time offset in seconds
% Pitch               N*2    int      array of pitch measurements (hundreths of deg)
% Roll                N*2    int      array of roll measurements (hundreths of deg)
% Heave               N*2    int      array of heave measurements (hundreths of deg)
% Heading             N*2    unit     array of heading measurements (hundreths of deg)

ATT_Time = fread(fid,2,'int32');
Num_meas = fread(fid,1,'int16');

for i = 1:Num_meas
ATT_Offset(i,:) = fread(fid,1,'int16');
end
for i = 1:Num_meas
ATT_Pitch(i,:) = fread(fid,1,'int16');
end
for i = 1:Num_meas
ATT_Roll(i,:) = fread(fid,1,'int16');
end
for i = 1:Num_meas
ATT_Heave(i,:) = fread(fid,1,'int16');
end
for i = 1:Num_meas
ATT_Heading(i,:) = fread(fid,1,'uint16');
end