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