%{
LOSC Matlab Sample

Example script to load a LOSC data file in Matlab
Tested with MATLAB R2014a
For more documentation on the structure of the LOSC data files,
see the LOSC tutorials:
https://losc.ligo.org/tutorials/

Last updated June 2016 by the LOSC team
%}
close('all')
filename = 'L-L1_LOSC_4_V1-842682368-4096.hdf5';

%-- Print out the structure of the file
h5disp(filename)

% -- GET THE STRAIN VALUES 
strain = h5read(filename, '/strain/Strain');
ts = h5readatt(filename,'/strain/Strain','Xspacing');

% -- Compute the sample frequency
fs = 1/ts;

% -- Get some meta data
gpsstart = h5read(filename, '/meta/GPSstart');
duration = h5read(filename, '/meta/Duration');
ifo      = h5read(filename, '/meta/Detector');

% -- Get some DQ info
dqnames = h5read(filename, '/quality/simple/DQShortnames');
dqmask  = h5read(filename, '/quality/simple/DQmask');

% -- Convert the DQ mask to a bit array
% -- Each column of the arra is now a 1 Hz time-seres for a DQ flag 
dqbitarray = de2bi(dqmask);

% -- This will be a hash table where each entry is the 1 Hz 
% -- time series for a DQ flag 
dqdict = containers.Map;

% -- Plot each DQ 1 Hz time series
for n = 1:length(dqnames)
    name = dqnames(n);
    chname = deblank(name{1});
    namespl = strsplit(chname, '_');
    timeseries = dqbitarray(:,n);
    dqdict(chname) = timeseries;
    figure();
    plot(timeseries, 'LineWidth',2);
    ylabel(namespl);
    xlabel(['Time since ',num2str(gpsstart),' (s)'])
    axis([0, 4096, -1, 2]);
end

% -- Construct a segment list
sci = dqdict('DATA');
sci_change = diff(sci);
startList = find(sci_change == 1);
stopList  = find(sci_change == -1);

if sci(1) == 1
    startList = [1, startList];
end
if sci(end) == 1
    stopList = [stopList, 4096];
end

% -- Load data for first segment
start = startList(1);
stop  = stopList(1);
goodData = strain(start*fs:stop*fs);
disp(['Loaded some good data with ', num2str(length(goodData)), ' samples.'])

% -- Plot first 2 seconds of data segment
N = fs*2;
time = [start:ts:stop];
figure()
plot(time(1:N), goodData(1:N))
xlabel(['Time since ', num2str(gpsstart), ' (s)'])
ylabel('Strain')




    


