Contents

classdef ts_channel
    % Time series data storage and manipulation.  Each channel can have its own start
    % time and sampling rate.  Multiple channels are combined in a ts_data
    % object.
    %
    % Data are uniformly spaced along the y-axis (time, frequency, etc)
    %
    % ts_channel properties:
    %   samprate      - Samples per second                    (double)
    %   starttime     - Time of first sample                  (ts_datetime)
    %   data          - Data storage                          (double vector)
    %   info          - Includes station, component ...       (ts_metadata)
    %   sourcefile    - Information about source data file    (ts_file)
    %   subplotinfo   -                                       (ts_subplotinfo)
    %   sampinterval  - READONLY: Interval between samples    (double)
    %                            (seconds)  ==1/samprate;
    %   nSamples      - READONLY: Number of samples in data   (double)
    %   nSeconds      - READONLY: Number of seconds in data   (double)
    %   endtime       - READONLY: time of last sample         (ts_datetime)
    %                             (NOT DRIFT CORRECTED)
    %   endtime_dc    - READONLY: time of last sample         (ts_datetime)
    %                             (DRIFT CORRECTED)
    %
    % ts_channel methods:
    %   Creator:
    %      ts_channel      - Creator
    %   Data modifiers:
    %      demean          - Remove mean of data
    %      detrend         - Remove mean plus trend
    %      filter          - Butterworth filter
    %      resample        - Resample data
    %      trueAmplitude   - Correct for instrument response
    %   Information:
    %      calcNumSamps    - Calculate the samples in a given number of seconds, hours, etc
    %      calcOffset      - Calculate time offset between two ts_channel objects
    %      isempty         - Is data vector is empty?
    %      isConsecutive   - Returns 1 if two data channels are consecutive in time
    %      isNotConsecutive - Returns 1 if two ts_channel objects are NOT consecutive
    %      overLapSeconds  - Return overlap between two channels, in seconds
    %      overLapSamples  - Gives overlap between two channels, in samples
    %   Transformation
    %      cat             - Concatenate two ts_channel objects
    %      cut             - Return a section of the data
    %      calcSpectrogram - Make a spectrogram
    %   Redefine Matlab methods:
    %      plot            - Plot a data channel
    %
    % EXAMPLE
    % sps=40; t=0:1/sps:100;
    % x=sin(t/5);
    % d=ts_channel(x,sps);
    % plot(d);
    % d
    % d.starttime='2007/01/01 00:00:00';
    % plot(d);

    properties
        % Samples per second
        samprate=[];
        % Exact time of first sample
        starttime=ts_datetime;
        % Data storage (double vector)
        data=[];
        % Includes station, component ...
        info=ts_metadata();
        % (ts_subplotinfo object)
        subplotinfo = ts_get('subplotinfo');
    end
    properties (Dependent=true)

        sourcefile     % Information about source data file
        sampinterval   % Interval between samples (seconds)  ==1/samprate;
        nSamples       % Number of samples in data
        nSeconds       % Number of seconds in data
        endtime        % Time of last sample in channel (NOT DRIFT CORRECTED)
        endtime_dc     % Time of last sample in channel (DRIFT CORRECTED)
    end
    properties (Access=private)
        PrivateSourcefile=ts_file();    % Private so that can set other properties
                                        % when setting sourcefile.  Dependent
                                        % property "sourcefile" is its public face
    end

CREATOR METHOD

obj=ts_channel(data,samprate)
obj=ts_channel(data,ts_file)
    methods
        function obj=ts_channel(data,arg2)
            % Create a ts_channel object
            %
            % Usage:
            %   obj=ts_channel(data,samprate) % Sets only data and sampling rate
            %       -or-
            %   obj=ts_channel(data,ts_file)  % Copies over all the parameters in the ts_file
            %                                 % ts_file.readchans must be an integer,  which
            %                                 % tells ts_channel which metadata to use
            %
            % Example
            %      tchan = ts_channel(sin([1:200]/20),10);
            ts_init ;	% Sets up global variables
            if nargin==0; return; end   % Handle zero-argument case
            if nargin~=2, error('ts_channel() needs two arguments, the data and the sampling rate'); end
            if isnumeric(arg2) && length(arg2)==1,
                obj.samprate=arg2;
            elseif isa(arg2,'ts_file');
                obj.sourcefile=arg2;  % set.sourcefile also sets the related obj properties
            else
                error('samprate must be a single numerical value or ts_file object');
            end
            if isnumeric(data)
                [M,N]=size(data);
                if M==1, obj.data=data;
                elseif N==1, obj.data=data';
                else error('data is %dx%d, must be a row or column vector',M,N);
                end
            else
                error('data is not numeric: type %s',class(data));
            end
        end
ans = 

       samprate: 
      starttime: NULL ('ts_datetime' object)
       nSamples: 0
           info: 'ts_metadata' object
                    station: '      ' 
                    network: '      ' 
                information: '      ' 
                   timezero:      NULL
                   tzoffset:        0 
                 clockdrift:        0 
                        lon:        0 
                        lat:        0 
                       elev:        0 
                     wdepth:        0 
                  component: '      ' 
               response (PxZ):  0x0     

     sourcefile: 'ts_file' object

METHODS THAT MODIFY THE DATA

filter(obj,ftype,npoles,lfreq,hfreq) % Butterworth filter
resample(obj,newsamprate)            % Resample data
trueAmplitude(obj,minfreq)           % Correct for instrument response
demean(obj)                          % Remove mean of data
detrend(obj,N)                       % Remove mean plus linear (or higher order) trend
        obj = filter(obj,ftype,npoles,lfreq,hfreq); % DEFINED IN FILE
        [obj,leftover] = resample(obj,newsamprate); % DEFINED IN FILE
        function obj = trueAmplitude(obj,minfreq)
            %Convert data to true amplitude
            %
            % Uses the object's stored instrument response
            % Usage:
            %	obj = trueAmplitude(obj);
            %	 -or-
            %	obj = trueAmplitude(obj,minfreq);
            %
            %	Sets new instrument response to unity
            %	Convolves the data with the inverse fft of the filter parameters.
            %	If minfreq isn't specified, the ifft is calculated for frequencies
            %   down to 1/2000 of the sampling rate.

            % 2/2006 Wayne Crawford

            if nargout ==0, error('Must have an output argument'); end
            ts_verbose(1,'Converting data to True Amplitude.');
            sr = obj.samprate;
            datalen = obj.nSamples;
            if exist('minfreq','var'); numFreqs=ceil((sr/2)/minfreq);
            else numFreqs=1000; end

            freqs = (1:numFreqs)*(sr/(2*numFreqs));
            oldResp = obj.info.response;
            d = obj.data;

            % Calculate the convolution operator to apply to the channel
            filtresp = oldResp.response(freqs).^(-1);
            filtresp2side = [0 filtresp fliplr(conj(filtresp(1:end-1)))];
            ifilter = (fftshift(ifft(filtresp2side)));	% Time-domain equivalent of the transfer function
            filtlen = length(ifilter);
            offset = filtlen/2;
            iadds = offset+1:datalen+offset;

            % Remove mean and trend from original data. Reduces low freq noise
            x=1:length(d); p=polyfit(x,d,1); d = d-polyval(p,x);

            truAmpl = conv(d,ifilter);		% Result is length(drive0)+length(ifilter)-1
            truAmpl = truAmpl(iadds)-mean(truAmpl(iadds));	% Extract only the original times

            obj.data = real(truAmpl);
            units = oldResp.inUnits;
            obj.info.response=ts_response('PZ',units,units,1,[],[]);

            ts_verbose(1,'Done');
        end %trueAmplitude
        function obj = demean(obj)
            % Remove mean of data
            obj.data=obj.data-mean(obj.data);
        end
        function obj = detrend(obj,N)
            % Remove mean and linear (or higher order) trend
            % Usage: obj=detrend(obj)   % removes mean and linear trend
            %        obj=detrend(obj,N) % removes trend up to order N
            %        obj=detrend(obj,1) % same as detrend(obj)
            if ~exist('N','var'); N=1; end
            x=1:length(obj.data);
            P=polyfit(x,obj.data,N);
            obj.data=obj.data-polyval(P,x);
        end

ACTION METHODS

plot(obj,...)                     % Plot data
calcNumSamps(obj,plotLength)      % Calculate the samples in a given number of seconds, hours, etc
calcSpectrogram(obj,windowlen,stepsize,...) % Make a spectrogram
cut(obj,start,len)                % Return a section of the data
cat(obj,objB)                     % Concatenate two data channels
        function [numSamps,suffix] = calcNumSamps(obj,plotLength)
            %Calculate the samples in a given number of seconds, hours, etc
            %
            % Usage: numSamps = calcnSamples(tdata,plotLength)
            % Input:
            %	tdata	  ts_data object
            %	plotLength character array with a number terminated by 's' (seconds),
            %	          'm' (minutes), 'h' (hours) or 'd' (days). If there is no
            %	          termination, the value is assumed to be in seconds
            %				If a scalar, assumed to be in seconds
            %
            %	Doesn't check to make sure that we don't go beyond the end of data
            %
            %	returns -1 if the input is screwy
            %
            % Output:
            %   the number of samples
            suffix = '';
            sampfreq = obj.samprate;
            if isa(plotLength,'char'),
                [A] = sscanf(plotLength,'%f %c');
                if length(A) > 1,
                    switch char(A(2)),
                        case 'm', numSamps = round(A(1) * sampfreq*60); suffix = 'm';
                        case 'h', numSamps = round(A(1) * sampfreq*3600); suffix = 'h';
                        case 'd', numSamps = round(A(1) * sampfreq*86400); suffix = 'd';
                        case 'p', numSamps = round(A(1)); suffix = 'p';
                        case 's', numSamps = round(A(1) * sampfreq);	 suffix = 's';
                        otherwise, ts_warning ('bad plot Length'); numSamps = -1;
                    end
                else numSamps = round(A(1))* sampfreq;
                end
            else numSamps = plotLength* sampfreq;
            end
            ts_debug(1,'ts_data/calcNumSamps: %s = %d samps',plotLength,numSamps);
        end
        [axhand,linehand] = plot(obj,varargin); % DEFINED IN FILE
        [sgrams,info]=calcSpectrogram(tdata,windowlen,stepsize,varargin); % MAKE A SPECTROGRAM
        obj=cut(obj,start,len); % DEFINED IN FILE
       function obj=cat(obj,obj2)
            % Concatenate two data channels
            %   fails if the offset between the data is more than the data length
            %   in the first objet
            if obj.samprate ~= obj2.samprate,
                error('Cannot concatenate objects with different sampling rates');
            end
            if obj2.starttime==ts_datetime(),    % time not set
                obj2.starttime=calcTime(obj,1,obj.nSamples+1);
            end
            if obj2.starttime <= obj.starttime
                error('The second data starts BEFORE the first!')
            end
            timea = calcTime(obj,obj.nSamples);
            timeb = calcTime(obj2,1);
            timediff = timeb-timea;
            dgap = timediff*obj.samprate - 1;
            if dgap > 0,
                ts_verbose(0,'%g sample excessive gap between datasets',dgap);
            elseif dgap<0,
                ts_verbose(0,'%g sample overlap between ',-dgap);
            end
            dgap = round(dgap);
            if dgap> obj.nSamples, error('Gap between data is more than nSamples'); end

            % Concatenate data
            totsamps = obj.nSamples+obj2.nSamples+dgap;
            newdata = zeros(1,totsamps);	% Predefine vector, already has zeros for gap (if any)
            istart2=totsamps-obj2.nSamples+1;
            % Put first data at beginning
            newdata(1:obj.nSamples) = obj.data;
            % Put second data at appropriate position
            newdata(istart2 : totsamps) = obj2.data;

            obj.data=newdata;
        end

INFORMATION METHODS

isempty(obj)                      % Returns 1 if data vector is empty
overLapSeconds(obj,objB)          % Returns overlap between two channels, in seconds
overLapSamples(obj,objB)          % Gives overlap between two channels, in samples
isConsecutive(obj,objB,tolerance) % Returns 1 if two data channels are consecutive
isNotConsecutive(obj,objB)        % Returns 1 if two data channels are NOT consecutive
cat(obj,objB)                     % Concatenate two data channels
        function out=isempty(obj)
            % Returns 1 if the data array is empty
            out=isempty(obj.data);
        end
        function out=overLapSeconds(objA,objB)
            % returns offset in seconds between perfect data overlap
            % objA is assumed to precede objB
            %   Positive values mean the start of objB is later than expected
            %   Negative values mean the start of objB is earlier than expected
            out=objB.starttime-objA.endtime_dc-objA.sampInterval;
        end
        function out=overLapSamples(objA,objB)
            % returns offset in samples from perfect data overlap
            % objA is assumed to precede objB
            %   Positive values mean the start of objB is later than expected
            %   Negative values mean the start of objB is earlier than expected
            % NaN means objA and objB have different sample rates, or
            % 	objB starts before objA
            if objA.samprate~=objB.samprate,
                ts_warning('the two channels have different sample rates');
                out=NaN;
                return;
            end
            if objB.startime < objA.starttime
                ts_warning('the second channel starts before the first');
                out=NaN;
                return;
            end
            out=overLapSeconds(objA,objB)*objA.samprate;
        end

        function out=isConsecutive(obj,obj2,tolerance)
            % Returns 1 if data have same samplerate and obj2 starts 1
            % sample after obj1
            % tolerance: tolerance for gap ~=1 (in samples) [0.5]

            if ~exist('tolerance','var'); tolerance = 0.5; end
            if tolerance < 0, error('tolerance cannot be less than 0'); end
            if tolerance > 100, ts_warning('you chose a tolerance > 100 samples!'); end
            % Returns 1 if data are consecutive in time, 0 if they are not or we can't tell (e.g., if starttime wasn't set)
            if obj.starttime==ts_datetime() || obj2.starttime==ts_datetime()    % Not consecutive if time is unset
                out=0; return;
            end
            if obj.samprate~=obj2.samprate, out=0; return; end  % Not consecutive if diff samp rates
            timediff = calcTime(obj2,1)-calcTime(obj,obj.nSamples);
            dgap = timediff*obj.samprate;    % # of data samples off in gap
            if dgap~=1,
                ts_verbose(2,'%g-sample gap between data',dgap);
            end
            out = abs(dgap-1)<tolerance;
        end

        function out=isNotConsecutive(obj,obj2)
            % Returns 1 if data are not consecutive in time, 0 if they are or we can't tell
            if obj.starttime==ts_datetime() || obj2.starttime==ts_datetime()
                out=0; return;
            end
            out=~isConsecutive(obj,obj2);
        end

REPLACEMENT METHODS

end(obj,k)                        % Implements 'end' in .() notation
char(obj,padding);                % Returns character string describing the object
disp(obj)                         % Called when you type the object without a semicolon
        function iend = end(obj,k)
            %Implement use of 'end' in .() notation
            if k == 1,
                iend = obj.nChannels;
            else
                iend = obj.nSamples;
            end
        end
        function otxt=char(obj,padding)
            % Make character string describing the object
            % padding pads the left side, to give an offset
            if ~exist('padding','var'), padding=0; end
            pformat=sprintf('%%%ds',padding);
            pad=sprintf(pformat,'');

            otxt=      sprintf('%s%15s: %g\n',                 pad,'samprate',obj.samprate);
            otxt=[otxt sprintf('%s%15s: %s (''%s'' object)\n',pad,'starttime',char(obj.starttime),class(obj.starttime))];
            otxt=[otxt sprintf('%s%15s: %d\n',           pad,'nSamples',obj.nSamples)];
            otxt=[otxt sprintf('%s%15s: ''%s'' object\n',pad,'info',class(obj.info))];
            otxt=[otxt sprintf('%s\n',                 char(obj.info,15+padding))];
            otxt=[otxt sprintf('%s%15s: ''%s'' object\n',pad,'sourcefile',class(obj.sourcefile))];
            %otxt=[otxt sprintf('%s',                   '',char(obj.sourcefile,15+padding))];
        end
        function disp(obj)
            % Displays object properties when object is typed without semicolon
            disp(char(obj));
        end

        % GET METHODS (FOR DEPENDENT VARIABLES)
        function obj = set.samprate(obj,val)
            if length(val)>1, error('samprate must be single numeric value'); end
            obj.samprate=val;
        end
        function obj = set.starttime(obj,val)
            val=ts_datetime(val);   % Converts strings, verifies the rest
            obj.starttime=val;
        end
        function obj = set.data(obj,val)
            if ~isnumeric(val), error('value must be numeric'); end
            [M,N]=size(val);
            if M~=1,
                if N==1, val=val';  % transpose if it was a column matrix
                else error('value must be a row or column vector');
                end
            end
            obj.data=val;
        end
        function obj = set.info(obj,value)
            if ~isa(value,'ts_metadata')
                error('value is not a ts_metadata object')
            end
            obj.info=value;
        end
        function sfile = get.sourcefile(obj)
            sfile=obj.PrivateSourcefile;
        end
        function obj = set.sourcefile(obj,finfo)
            % Reads ts_file  into obj.sourcefile and directly copies samprate and metadata
            % MATLAB doesn't like me coping directly into another property
            if ~isa(finfo,'ts_file'); error('value must be a ts_file object'), end
            i=finfo.readchans;
            if length(i)~=1,
                error('when setting a ts_channel''s sourcefile parameter, the passed ts_file readchans should be a scalar');
            end
            obj.PrivateSourcefile=finfo;
            %finfo   % DEBUG

            if ~isempty(finfo.samprate); obj.samprate=finfo.samprate; end
            if ~isempty(finfo.channelinfo(i)); obj.info=finfo.channelinfo(i); end
        end
        function obj = set.subplotinfo(obj,val)
            if ~isa(val,'ts_subplotparm'), error('value must be a ts_subplotparm object'); end
            obj.subplotinfo=val;
        end

        % GET METHODS (FOR DEPENDENT VARIABLES)
        function out =  get.sampinterval(obj)
            out=1/obj.samprate;
        end
        function ns =  get.nSamples(obj)
            ns=length(obj.data);
        end
        function ns =  get.nSeconds(obj)
            ns=obj.nSamples/obj.samprate;
        end
        function out = get.endtime(obj) % Time of last sample (not drift corrected)
            out=obj.starttime+(obj.nSamples-1)/obj.samprate;
        end
        function out = get.endtime_dc(obj) % Drift corrected time of last sample
            out=obj.GPSTime(obj.endtime);
        end
    end % Public methods

STATIC METHODS

Called using ts_channel.methodname

calcOffset(timeOffset,driftRate,samprate)
    methods (Static=true) % Does not require an object argument. Call using ts_channel.methodname
        function addrOffset= calcOffset(timeOffset,driftRate,samprate)
            % Return the address offset equalling a given time offset
            %  , for the specified drift rate.
            %
            % Usage: caddrOffset= calcAddrOffset(timeOffset,driftRate,samprate)
            %
            % Inputs:
            %	timeOffset: time (in seconds from current position)
            %	driftRate is the instrument drift rate:   (instElapsed-GPSElapsed)/GPSElapsed
            %   sampRate is the instrument sampling rate (samples/second)
            % Output:
            %	addrOffset: the address offset from the current posision

            instTimeOffset = ts_datatime.driftCorrecSecs(timeOffset,driftRate,'backwards'); % Convert to UNCORRECTED seconds
            addrOffset = round(instTimeOffset*samprate);
        end
    end % static methods

PRIVATE METHODS

Accessed by class members only. Use only if you are changing a TiSKit class

obj=cut_addr(obj,istart,iend)
currtime = calcTime(obj,address)
[address,delay] = calcAddress(obj,time)
addrList = makeAddressList(obj,addressString,addBuffer,option)
    methods (Access=private) % Access by class methods only
        obj=cut_addr(obj,istart,iend)  % DEFINED IN FILE
        function currtime = calcTime(obj,address)
            % Calculate the time at a given address on a given channel
            %
            % currtime = calcTime(obj,address);
            % where address the address of the current sample if we assume that address 1
            %   equals the starttime
            %
            % Assumes that obj.starttime corresponds to the CORRECT (drift
            % corrected) time at address 1

            % NEW WAY USES TS_FILE (WCC 02/2011)
            insttime=obj.starttime+(address-1)/obj.samprate;
            temp=obj.info;
            temp.timezero=obj.starttime;
            currtime=temp.GPSTime(insttime);

        end
        function [address,delay] = calcAddress(obj,time)
            % Calculate the closest address to the given time
            %
            % Usage: [address, delay] = calcAddress(obj,time);
            % Input:
            %   time    ts_datetime object giving the time to look for (will drift correct)
            %			-or-
            %			scalar giving seconds from start of data (won't drift correct)
            %
            % Output:
            %   address closest to the requested time (first address is number 1)
            %	delay is the # of seconds the first sample is AFTER the requested time
            %
            %  Note: the address may be outside of the current file/record!
            %
            % Corrects for the clock drift (SHOULD I REALLY DO THIS HERE?)

            ts_debug(1,'Starting...');
            if isa(time,'ts_datetime')
                temp=obj.info;
                temp.timezero=obj.starttime;
                insttime=temp.instrumentTime(time);
                toffset=insttime-obj.starttime;
            else
                toffset=time;
            end

            srate = obj.samprate;
            faddoffset=toffset*srate;
            address = 1 + floor(faddoffset+.5);
            delay = (address - faddoffset - 1)/srate;
            ts_debug(1,'Ending, address=%d, delay=%g...',address,delay);
        end
        function addrList = makeAddressList(obj,addressString,addBuffer,option)
            %MAKEADDRESSLIST GRANDFATHERED Convert a Matlab address string to numbers (should be private?)
            %
            % Usage: addrList = makeaddresslist(obj,addressString,addBuffer)
            %
            % Takes a string of addresses in Matlab format and convert it to numerical values.
            % Optionally, adds a buffer to the back end of the addresses (for causal filtering)
            % This is the only smart option, because you could just use the numerical
            % array otherwise (except that the string list is shorter and takes into
            % account the data length)
            %
            % For example, if addressString = '1:1000 10000:40000 60000:end',
            % length(obj)=100000 and addBuffer=0, then addrList will be
            %	[1:1000 10000:40000 60000:100000];
            % for the same inputs, but addBuffer=1000, then addrList will be
            %	[1:2000 9000:41000 59000:100000];
            %
            %	option[0]:	0, output results as a numerical array
            %				1, output results as an array of cell arrays containing
            %					{startaddress, length}

            if ~exist('addBuffer','var'), addBuffer=0; end
            if ~exist('option','var'), option=0; end

            addressString = strrep(addressString,'end',num2str(obj.nSamples));
            if addBuffer > 0, addressString = ['1 ' badaddrs]; end	% Will fill out to 1:addBuffer below
            iElem=0;
            if ~isempty(addressString)
                if addBuffer > 0 || option==1,
                    A=textscan(addressString,'%s');
                    addressString = [];
                    for i=1:length(A{1}),
                        T=A{1}{i};
                        if ~isempty(T),
                            k = strfind(T,':');
                            if isempty(k),		% If just a single value
                                minaddr = str2double(T);
                                maxaddr = str2double(T)+addBuffer;
                                if option==1,
                                    addrList{iElem} = {minaddr,maxaddr-minaddr+1};
                                end
                            else
                                B = sscanf(T,'%d:%d');
                                minaddr = B(1);
                                maxaddr = B(2)+addBuffer;
                            end
                        end
                        if option==1,
                            iElem=iElem+1;
                            addrList{iElem} = {minaddr,maxaddr-minaddr+1};
                        else
                            addressString = [addressString sprintf(' %d:%d',minaddr,maxaddr)];
                        end
                    end
                end
            else
                if option==1,
                    addrList = {};
                end
            end
            if option==0,
                eval(['addrList = [' addressString '];']);
            end
        end
    end % Private methods
end % classdef