function dtx=rddtx(fil,numdec);
% RDDTX Reads .DT4 Files generated by a Biosonics DTX sounder
%   DTX=RDDTX(fils) reads sequentially the filenames given in the cell array fils (which can
%                   be a single file). This allows 'chaining' of long data segments which
%                   were originally split to keep file sizes reasonable.
%
%   DTX=RDDTX(fils,decnum) will decimate the files by taking the median-in-time
%                   over decnum samples (default 3). This is very useful in 
%                   removing interference from other sounders.
%
%  DTX is a structure whose contents are generally self-explanatory (but see notes below)
%  
%  Notes: GPS data is read but there are two ways to get this data. Internally
%        GPS location (but not times) are decoded and stored if you feed a 
%        NMEA datastream to the sounder. In addition, NMEA serial
%        strings are also stored "raw". Since these also contain TIMES they
%        are generally more useful for synchronizing with other datasets. I 
%        decode GGA and RMC strings, but if you only have (e.g.) GLL then you 
%        might want to add a bit more code. Note that if you use strings without the DATE
%        (like GGA) the date in the decoded files will be guessed (possibly incorrectly)
%        based on the clock of the laptop being used to log data.
%
%        Also, I remove duplicate GPS positions which may cause problems with simulated
%       data.
%
%        Two different timestamps are returned:
%             .imtime - internal sounder clock (resolution 1 ms)
%             .mtime  - time stamps from NMEA strings, interpolated to ping times.
%
%        At a pulse rate of > 1Hz you will see that imtime is very regular compared
%        with mtime (i.e. good short-term relative accuracy), but can have offsets and long-term
%        drift compared with correct time. mtime is better in absolute terms (i.e. correct 
%        time overall with no drift, but more "jitter" due to 1 sec resolution in input strings).
%        If you care about time stamps to that level of precision you should correct imtime with
%        (e.g.) a linear or quadratic drift to make it match mtime.
%

% Rich Pawlowicz (rich@eos.ubc.ca)  Sep/2003
%      25/5/2006 - added RMC decoding.

if ~iscell(fil), fil={fil}; end;
 


nfil=length(fil);
ifil=1;

fprintf('\nOpening file %s\n',fil{ifil});
fd=fopen(fil{ifil},'r','ieee-le');


if nargin<2, numdec=3; end;

% Removes duplicates in GPS input strings.
NODUP=1;


idec=0;
numrec=0;
numgps=0;

eof=0;

while ~eof,

  % Read the record tag and check size
  
  [tag,siz]=rd_tag(fd);
  
  bits=fread(fd,siz+2,'uint8');
  if  (bits(end)*256+bits(end-1))-4 ~=length(bits),
    disp('Length of record in error');
  end;  

  %%fprintf('%s - %d\n',tag,siz);

  % Decode the tage
  switch tag,
    case 'FFFF', % Signature
	rd_sig(bits);
	
    case '1E00', % Environment
        env=rd_env(bits);
    case '1200', % Sounder settings
        snd=rd_snd(bits);
	
	if ~exist('den'),

	   den=struct('imtime',repmat(NaN,1,numdec),...
	              'pingnum',repmat(NaN,1,numdec),...
		      'vals',repmat(NaN,snd.sampperping,numdec));

	   dtx=struct('name','dtx','mtime',repmat(NaN,1,ceil(snd.npings/numdec*nfil)),...
	              'env',env,'snd',snd,...
	              'imtime',repmat(NaN,1,ceil(snd.npings/numdec*nfil)),...
	              'pingnum',repmat(NaN,1,ceil(snd.npings/numdec*nfil)),...
		      'vals',repmat(NaN,snd.sampperping,ceil(snd.npings/numdec*nfil)),...
		      'range',([1:snd.sampperping]'+snd.blank)*snd.sampperiod*env.sv/2e6);

	   pos=struct('name','gps','latitude',repmat(NaN,1,snd.npings),...
	              'longitude',repmat(NaN,1,snd.npings),...
		      'mtime',repmat(NaN,1,snd.npings),'ptime',repmat(NaN,1,snd.npings));
	end;
		   	   
    case '1300', % Pulse settings
        pulse=rd_pulse(bits);
	
    case '1000', % Marker
        mark=rd_mark(bits);
	
    case '0000', % NULL
    case '1500', % SIngle beam ping
        ping=rd_ping(bits);
	
	if ~exist('dtx'), 
	  fprintf('ping data before environmental!\n');
	else
	  
	  if numdec==1;
	    numrec=numrec+1;
	    dtx.pingnum(numrec)=ping.pingnum;
	    dtx.vals(1:length(ping.samps),numrec)=log10(max(ping.samps,1));
	    dtx.vals((length(ping.samps)+1):end,numrec)=NaN;
	    dtx.imtime(numrec)=chron.mtime;
      	    if rem(numrec,30)==0, fprintf('.'); end;
	    if rem(numrec,30*50)==0, fprintf('\n'); end;
	  else  
	    idec=idec+1;
	    den.pingnum(idec)=ping.pingnum;
	    den.vals(1:length(ping.samps),idec)=ping.samps;
	    den.vals((length(ping.samps)+1):end,idec)=NaN;
	    den.imtime(idec)=chron.mtime;
	    if idec==numdec,
	      idec=0;
	      numrec=numrec+1;
	      dtx.pingnum(numrec)=median(den.pingnum);
	      dtx.imtime(numrec)=median(den.imtime);
	      dtx.vals(:,numrec)=log10(max(median(den.vals,2),1));
      	      if rem(numrec,30)==0, fprintf('.'); end;
	      if rem(numrec,30*50)==0, fprintf('\n'); end;
	    end;
	  end;
	end;
	
    case '1C00', % dual beam ping
        fprintf('%s - %d\n',tag,siz);
       
    case '1D00', % Split beam ping
        fprintf('%s - %d\n',tag,siz);
       
    case '0F00', % Chron
         
         chron=rd_chron(bits);
	 
    case '1F00', % Motor
        fprintf('%s - %d\n',tag,siz);
       
    case '0E00', % GPS Position
        gps=rd_gps(bits);
	gps.mtime=chron.mtime; 
 	
    case '1100', % GPS string
        string=rd_gpsstring(bits);
	
	% Note that the GPS block doesn't contain times, so we have to parse
	% this from the NMEA strings which are copied over into the file.
	% 
	% There are various possible strings to use, some better than others.

        % GPRMC has complete time stamp and appears to have enough decimal
	% places in the lat/long string in recent GPS units.
	rmc=strfind(string,'$GPRMC');
	if any(rmc),
	  string=string(rmc(1):end);  %Strip out the RMC info
	  icom=find(string==',');
 	  numgps=numgps+1;
	  % Time from string
 	  pos.mtime(numgps)=datenum( string([icom(9)+1:icom(10)-1 icom(1)+[1:6]]) ,'ddmmyyHHMMSS'); 
 	  % Internal clock time from last chron block
	  pos.ptime(numgps)=chron.mtime;
	  % lat/long.
	  lat=sscanf(string(icom(3)+1:icom(4)),'%f');
	  pos.latitude(numgps) = fix(lat/100) + rem(lat,100)/60;
	  lat=sscanf(string(icom(5)+1:icom(6)),'%f');
	  pos.longitude(numgps)= -( fix(lat/100) + rem(lat,100)/60 );
	else
	  
	  % GPGGA has lat/long and hhmmss, but not date.
          gga=strfind(string,'$GPGGA');
	  if any(gga),
	    string=string(gga(1):end);  %Strip out the GGA
	    icom=find(string==',');
	    numgps=numgps+1;
	    
	    % Take the date from the chron block, and the hhmmss from the NMEA string. This may screw
	    % up if clocks are not synchronized correctly.
	    % Fix 25/Jul/06 - use *last* date in general, take only first unknown from chron block.
	  %%  pos.mtime(numgps)=fix(chron.mtime)+datenum([string(icom(1)+[1:2]) ':' string(icom(1)+[3:4]) ':' string(icom(1)+[5:6])]);
	    if numgps>1,
	      pos.mtime(numgps)=fix(pos.mtime(numgps-1))+datenum(['00010000' string(icom(1)+[1:6])],'ddmmyyyyHHMMSS');
	    else
  	      pos.mtime(numgps)=fix(chron.mtime)+datenum(['00010000' string(icom(1)+[1:6])],'ddmmyyyyHHMMSS');
	    end;  
	%% [datestr(fix(chron.mtime)) '  ' datestr(pos.mtime(numgps))]
 	    pos.ptime(numgps)=chron.mtime;
	    lat=sscanf(string(icom(2)+1:icom(3)),'%f');
	    pos.latitude(numgps) = fix(lat/100) + rem(lat,100)/60;
	    lat=sscanf(string(icom(4)+1:icom(5)),'%f');
	    pos.longitude(numgps)= -( fix(lat/100) + rem(lat,100)/60 );
	  end
	end;
	  
	 
    case 'FEFF', % EOF
      fclose(fd);
      ifil=ifil+1;
      if ifil<=length(fil), 
        fprintf('\nOpening file %s\n',fil{ifil});
        fd=fopen(fil{ifil},'r','ieee-le'); 
      else 
        fd=-1; 
      end;
      if fd<0, eof=1; fprintf('No more files found!!!\n'); end;
      
  end;
 
end;

% Get rid of blanks at end of record

kk=finite(dtx.imtime);
dtx.imtime=dtx.imtime(kk);
dtx.pingnum=dtx.pingnum(kk);
dtx.vals=dtx.vals(:,kk);

% Interpolate in lat/long info to actual ping times.
% What if there is no GPS? - then don't do anything!
% Note that I ignore the GPS blocks.

%subplot(211);
%plot(pos.ptime);
%subplot(212);
%plot(pos.mtime);
%%plot((pos.ptime-pos.mtime-(12-7)/24)*86400,'.-');
 
if numgps>2,

 % Filtering gets rid of some of the 'bit noise'...but perhaps this is better done
 % by the user if required.
 % kk=finite(pos.mtime);
 % pos.mtime(kk)=filtfilt([1 1]/2,1,pos.mtime(kk));

  if NODUP
    kk=find([1,diff(pos.longitude)~=0 & diff(pos.latitude)~=0 ] & finite(pos.ptime));
    ptim=[0,pos.ptime(kk),max(max(pos.ptime),max(dtx.mtime))+1];
  else
    kk=1:numgps;
  end;
  
  if length(kk)<numgps,
    fprintf('\n******** Duplicate GPS positions found and removed\n');
  end;  
    
  dtx.mtime    =interp1( ptim,pos.mtime(kk([1 1:end end])),dtx.imtime);
  dtx.longitude=interp1( ptim,pos.longitude(kk([1 1:end end])),dtx.imtime);
  dtx.latitude =interp1( ptim,pos.latitude(kk([1 1:end end])) ,dtx.imtime);
else
  fprintf('\n******** No GPS data found\n');
end;




return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [tag,siz]=rd_tag(fd);

siz=fread(fd,1,'uint16');
tag=dec2hex(fread(fd,2,'uchar'),2);

tag=[tag(1,:) tag(2,:)];
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rd_sig(bits);
   

  
return
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function env=rd_env(bits);
 
  env.absorb      =        ( bits(1) + bits(2)*256 )/10000; 
  env.sv          =        ( bits(3) + bits(4)*256 )*.0025+1400; 
  env.temperature = twoscvt( bits(5) + bits(6)*256 )*.01; 
  env.salinity    =        ( bits(7) + bits(8)*256 )*.01; 
  env.power       = twoscvt( bits(8) + bits(10)*256 )*.1; 
  env.nsdr        = bits(11) + bits(12)*256; 

return;
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function snd=rd_snd(bits);
 
  
  snd.address     =        ( bits(1) + bits(2)*256 ); 
  snd.npings      = twoscvt( bits(3) + bits(4)*256+  bits(5)*65536 + bits(6)*16777216,4); 
  snd.sampperping =        ( bits(7) + bits(8)*256 ); 
  snd.sampperiod  =        ( bits(9) + bits(10)*256 )/1000; 
  snd.pulselen    =        ( bits(13) + bits(14)*256 )/1000; 
  snd.pingperiod  =        ( bits(15) + bits(16)*256 )/1000; % I'm guessing the scale factor
  snd.blank       =        ( bits(17) + bits(18)*256 ); 
  snd.maxdata     =        ( bits(19) + bits(20)*256 ); 
  
  snd.threshold   = twoscvt( bits(21) + bits(22)*256 )*.01; 
  
  snd.rxee        = bits(23+[0:127]);
  snd.txee        = bits(151+[0:127]);
  snd.ccor        =        ( bits(11) + bits(12)*256 )*.001; 

return;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pulse=rd_pulse(bits);
 
  switch twoscvt( bits(1) + bits(2)*256 ),
    case 0
      pulse.type='raw';
    case 1
      pulse.type='integrated';
    case 2
      pulse.type='chirp';
  end;    
  pulse.address     =        ( bits(3) + bits(4)*256 ); 

return;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function mark=rd_mark(bits);
 
  switch twoscvt( bits(1) + bits(2)*256 ),
    case 0
      mark.type='event';
    case 1
      mark.type='start';
    case 2
      mark.type='end';
  end;    
  mark.address     =        ( bits(3) + bits(4)*256 ); 

return;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function ping=rd_ping(bits);
  
  ping.address     =        ( bits(1) + bits(2)*256 ); 
  ping.pingnum     =        ( bits(3) + bits(4)*256+  bits(5)*65536 + bits(6)*16777216); 
  ping.ptime       =        ( bits(7) + bits(8)*256+  bits(9)*65536 + bits(10)*16777216);
  
  nsamp= ( bits(11) + bits(12)*256 );

  
  rle=[0;find(bits(14:2:end)==255);nsamp];   % Find RLE bits
  nrle=[-1;bits(11+2*rle(2:end-1))+2;-1];
  if 0,%length(rle)>2,
    fprintf('%d ',ping.pingnum);
    for ll=2:length(rle)-1,
     fprintf('rle = %d,%d  ',[rle(ll)',nrle(ll)']);
    end;
    fprintf('  nsamp=%d\n',nsamp);
  end;
  
  exponent=                        bitshift( bits(14:2:end), -4);
  mantissa=   bits(13:2:end-1)+256*bitand(   bits(14:2:end),15);
  
  il=exponent==0;
    samps(il)=mantissa(il);
    samps(~il)=bitshift( mantissa(~il)+hex2dec('1000'), exponent(~il)-1);
    
  nused=0;  
  for k=2:length(rle),
    ping.samps((rle(k-1)+1:rle(k)-1)+nused)=samps(rle(k-1)+1:rle(k)-1);
    if nrle(k)>0,
      ping.samps((rle(k)+[0:nrle(k)-1])+nused)=0;
      nused=nused+nrle(k)-1;
    end;  
  end;
  
return;
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
function chron=rd_chron(bits);



  chron.source=bit_source(bits(5));
  
  chron.ptime       =        ( bits(7) + bits(8)*256+  bits(9)*65536 + bits(10)*16777216);

  % Currently we have UNIX time, seconds since 1970
  chron.time=                ( bits(1) + bits(2)*256+  bits(3)*65536 + bits(4)*16777216);
  
  chron.mtime=chron.time/86400+datenum(1970,1,1) + rem(chron.ptime,1000)/86400000;


return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function source=bit_source(bit);
  switch dec2hex(bit,2),
    case '02',
      source='calendar';
    case '06',
      source='radio';
    case '09',
      source='chronometer';
    case '11',
      source='GPS';
    case '12',
      source='LORAN-C';
    otherwise
      fprintf('Unrecognized clock source - %s\n',dec2hex(bits(5)));
  end;
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function gps=rd_gps(bits);
 
  
  gps.latitude    = twoscvt( bits(1) + bits(2)*256+  bits(3)*65536 + bits(4)*16777216,4)*1e-5/60; 
  gps.longitude   = twoscvt( bits(5) + bits(6)*256+  bits(7)*65536 + bits(8)*16777216,4)*1e-5/60; 
  gps.source=bit_source(bits(9));

return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function string=rd_gpsstring(bits);

  string=setstr(bits');
  
return;

  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
function x = twoscvt( x, N )
%
% x = twoscvt( x, N )
%
% Converts N-byte integers in 2's complement form to numbers in the
% range [-(256^N)/2 : 256^N)/2-1].  N must be even, but is optional
% (the default value is 2 -i.e. 16-bit integers).
%
% For example:  [FFFF 8000 7FFF 0000] ---> [ -1  -32768  32767  0 ]


if( nargin < 2 ),  N = 2;  end         % use default if N is not specified

if( rem(N,2) )
   error('the number of bytes in the integers must be even')
end

M = 256^N;                             % FF....FF

i = find( x > M/2-1 );                 % x > 32767 for N = 2

x(i) = x(i) - M;