function varargout = read_tekisf(infile)
%
% function [ x, y, h ] = read_tekisf(infile)
%
% infile:  filename of the "isf" file downloaded from the oscilloscope
% x: time-data;   y: voltage[signal] data;   h: header
% 
% Read the waveform files in the "internal" data format recorded by 
% Tektronix oscilloscopes (e.g. tested on TDS3052 and DPO4102B) and 
% doubtless plenty of others).
%
% Giving no output arguments causes simple statistics of the data to 
% be displayed.   One output argument gives the measured ('y') data;
% two give the x then y data; three give also an additional header-structure.
%
% Nathaniel Taylor, partly adapted from a function by John Lipp (see comments).
% 2010-02-16 (specifically for Tek TDS3052)
% 2015-07-24 (adapted to Tek DPO/MSO4000B, but should still work for TDS3052)
% 

% 
fh = fopen(infile, 'r');
if fh<0,
	error('failed to open input file "%s" for reading', infile);
end
% read file as separate bytes
% [change to "uint8=>uint8" if wanting to be more space-efficient!]
bytes = fread(fh, Inf, 'uint8');

% convert to floating-point vectors
[x, y, h] = tekisf2num(bytes);

% make output variables
if nargout==0,
    header
    fprintf( 'x range: %g -- %g, [min %g, max %g]\n', ...
        x(1), x(end), min(x), max(x) );
    fprintf( 'y range: %g -- %g, [min %g, max %g]\n', ...
        y(1), y(end), min(y), max(y) );
    varargout = {};
elseif nargout==1,
    varargout = { y };
elseif nargout==2,
    varargout = { x, y };
elseif nargout==3,
    varargout = { x, y, h };
else
    help read_tekisf;
    error('expected from zero to three output args: y ; x,y ; x,y,header');
end


%----------------------------------------------------------------------------
function [x, y, h] = tekisf2num(bytes)
%
% function [x, y, h] = tekisf2num(bytes)
%
% tekisf2num: read data as a stream of bytes in the `internal' data format
% of Tektronix TDS3000-series oscilloscopes, converting to double-arrays
% of the x and y data.
%
% The original file 'isfread.m' from
%       http://www.mathworks.com/matlabcentral/fileexchange/6247
% was written 17/8/2004 by John Lipp - CCLRC Rutherford Appleton Laboratory.
% It was downloaded on 2009-05-04 and adapted by Nathaniel Taylor (KTH) for:
%       separate x,y output
%       variable number of output-arguments, for convenience
%       check success of opening input file
%       explicit big-endian mode (seems right for our scopes)
%       accept variable header length based on string #520000
%         (not sure how important this is; no guarantee the string is always there)
% It was further adapted on 2010-02-17 to accept a stream of bytes (e.g.
% reading data directly from the instrument) instead of from a file.
% Then, a lot more on 2015-07-24, to remove any assumption about the 
% fields in the header: now each semicolon-separated item is parsed as 
% a name-value pair, and the name is used as a fieldname for the header 
% structure output.
%

if ~isnumeric(bytes),
    error('must have numeric input (bytes from oscilloscope)');
end
if length(bytes)<2000,  % a very arbitrary choice of number
    error('expected many thousands of bytes of input: problem?');
end
bytes = reshape(bytes, 1, []);

% header (ascii)
% assumption:  the header [first part of file] is ascii characters, with
% no hash ('#') until the end of the header; at the end there is the 
% familiar method of specifying the following binary data by "#Nnnnn" 
% where the "N" is an ascii digit indicating how many ascii digits follow,
% then these following digits indicate how many binary bytes follow them;
% e.g. for the Tek3052 which always had 10k 2-byte samples to transfer, the 
% end of the header was followed by "#520000" indicating 20000 bytes (the 
% number "20000" takes 5 characters...);  for later oscilloscopes the 
% data length can vary
% 
% search at least this many bytes to find the end-of-header:
header_maxlength = 1000;  % (we usually see only about 500 bytes in the header)
% find where the marker occures
marker_position = find( bytes(1:header_maxlength) == double('#') );
if numel(marker_position)==0,
    error('failed to find end-of-header marker "#" in first %d bytes', ...
        header_maxlength);
end
% ignore any later bytes (in the data) that happen to be the same as the marker
marker_position = marker_position(1);
% take following digit to say how many more to read!
n_digit = str2double(char(bytes(marker_position+1)));
% so, find how many *more* (binary bytes) we should read
n_databytes = str2double(char(bytes(marker_position+1+(1:n_digit))));
% total header-length:
header_length = marker_position + 1 + n_digit;
% extract instrument settings etc from the header
h = parseHead( char(bytes(1:header_length)) );

% data (binary)
if length(bytes) ~= header_length + 2*h.NR_PT,  % the three are \n\r\n
    warning('expected total input bytes %d (%d+2*%d), but got %d', ...
         header_length + 2*h.NR_PT, ...
         header_length, h.NR_PT, length(bytes) );
end
inData = bytes((header_length+1):2:(header_length+2*h.NR_PT))*256 ...
       + bytes((header_length+2):2:(header_length+2*h.NR_PT));
inData(inData>=2^15) = inData(inData>=2^15) - 2^16;
lowerXLimit = h.XZERO;
upperXLimit = ((h.NR_PT-1)* h.XINCR + h.XZERO);
x = [lowerXLimit:h.XINCR:upperXLimit].';
y = h.YMULT*(inData-h.YOFF).';


%----------------------------------------------------------------------------
function h = parseHead(header_string)
%
% example of "header_string" content (from a DPO4102B, 2015)
% ;:WFMPRE:BYT_NR 2;BIT_NR 16;ENCDG BINARY;BN_FMT RI;BYT_OR MSB;WFID "Ch1, DC coupling, 1.000V/div, 20.00ms/div, 5000000 points, Sample mode";NR_PT 5000000;PT_FMT Y;XUNIT "s";XINCR 40.0000E-9;XZERO -99.8712E-3;PT_OFF 0;YUNIT "V";YMULT 156.2500E-6;YOFF 1.6640E+3;YZERO 0.0E+0;VSCALE 1.0000;HSCALE 20.0000E-3;VPOS 260.0000E-3;VOFFSET 0.0E+0;HDELAY 128.8000E-6;DOMAIN TIME;WFMTYPE ANALOG;CENTERFREQUENCY 0.0E+0;SPAN 0.0E+0;REFLEVEL 0.0E+0;:CURVE #810000000

h = struct;
% keep looping as long this variable still has ";" somewhere in it
while ~isempty(strfind(header_string,';')),
    % chop the part up to the first ";", leaving only the *rest* in 'h'
    [part, header_string] = strtok(header_string,';');
    % of that first chopped part, separate the name-value pair
    [name, value] = strtok(part);
    % try to convert value to a number
    numvalue = str2double(value);
    % str2double gives numeric value 'NaN' if it can't convert
    if ~isnan(numvalue),
        value = numvalue;   % take numeric value if it worked
    else
        value = strtrim(value);  % else, take text after removing padding
    end
    name = regexprep(name,'[^\w]','');  % clean to make valid as a variable-name
    h.(name) = value;
end