diff --git a/classes/@SeismicTrace/legend.m b/classes/@SeismicTrace/legend.m index 0d06064..93301a6 100644 --- a/classes/@SeismicTrace/legend.m +++ b/classes/@SeismicTrace/legend.m @@ -1,7 +1,7 @@ function varargout = legend(T, varargin) - %legend Create a legend for a waveform graph. + %legend Create a legend for the graph of a trace. % legend(traces) attempts to automatically create a legend based upon - % unique values within the waveforms. in order, the legend will + % unique values within the traces. in order, the legend will % preferentially use station, channel, start time. % % legend(traces, field1, [field2, [..., fieldn]]) will create a legend, @@ -12,8 +12,8 @@ % location, etc.) % % For additional control, use matlab's legend function by passing it - % cells & strings instead of a waveform. - % (hint:useful functions include waveform/get, strcat, sprintf, num2str) + % cells & strings instead of a trace. + % (hint:useful functions include strcat, sprintf, num2str) % % See also plot, legend diff --git a/classes/@SeismicTrace/plot.m b/classes/@SeismicTrace/plot.m index d8ad963..ccc425d 100644 --- a/classes/@SeismicTrace/plot.m +++ b/classes/@SeismicTrace/plot.m @@ -142,15 +142,12 @@ p.KeepUnmatched = true; p.CaseSensitive = false; p.StructExpand = false; - if ismethod(p,'addParameter') + if verLessThan('matlab','8.2'); %r2013b + addParameter = @addParamValue; + end addParameter(p,'autoscale', false); addParameter(p,'xunit', 's'); addParameter(p,'fontsize', 10); - else % older usage: pre r2013b - addParamValue(p,'autoscale', false); %#ok<*NVREPL> - addParamValue(p,'xunit', 's'); - addParamValue(p,'fontsize', 10); - end p.parse(argList{:}); end @@ -268,4 +265,4 @@ newParams = horzcat(fieldList,v); end end -end \ No newline at end of file +end diff --git a/classes/@TraceFilter/TraceFilter.m b/classes/@TraceFilter/TraceFilter.m index 71f7ddf..5a0157f 100644 --- a/classes/@TraceFilter/TraceFilter.m +++ b/classes/@TraceFilter/TraceFilter.m @@ -1,17 +1,28 @@ classdef TraceFilter %TraceFilter Simple butterworth filtering for Traces + % + % TraceFilter Properties: + % type - describes the passband type of filter. ('B', 'L', or 'H') + % cutoff - upper and lower bounds for the filter + % poles - number of poles for the filter + % + % TraceFilter Methods: + % filtfilt - Zero-phase forward and reverse butterworth filtering. + % % replaces filterobject + % + % see also: butter, filtfilt properties type = 'B'; % could be 'B'andpass, 'H'ighpass, or 'L'lowpass cutoff = [0.8 5]; % filter cutoffs [lower bound, upper bound] - poles = 2; - + poles = 2; % number of poles end methods function f = TraceFilter(filtType, cutoff_, poles_) - %TRACEFILTER constructor for a filter object + %TRACEFILTER construct a TraceFilter + % switch nargin case 1 if isa(filtType,'filterobject') @@ -82,6 +93,10 @@ function assertCutoffCountMatchesType(f) % See also FILTFILT, BUTTER % Check, if w is a number array, just quick & dirty it... + if isa(trace,'waveform') + disp('converting waveform(s) to SeismicTrace(s)'); + trace = SeismicTrace(trace); + end if isnumeric(trace) WN = f.cutoff / NYQ; [b, a] = getButter(f,WN); @@ -89,7 +104,7 @@ function assertCutoffCountMatchesType(f) elseif isa(trace,'SeismicTrace') assert(numel(f) == 1, 'can only filter using a single TraceFilter at a time'); - assertCutoffCountMatchesType(); + assertCutoffCountMatchesType(f); HASGAP = arrayfun(@(x) any(isnan(x.data)), trace); if any(HASGAP(:)) warning('Filterobject:filtfilt:hasNaN',... @@ -161,4 +176,7 @@ function disp(f) end end end + methods(Static) + cookbook() %run the filter cookbook + end end diff --git a/classes/@TraceFilter/cookbook.m b/classes/@TraceFilter/cookbook.m new file mode 100644 index 0000000..624c234 --- /dev/null +++ b/classes/@TraceFilter/cookbook.m @@ -0,0 +1,8 @@ +function cookbook() + %cookbook Cookbook for TraceFilter + % Demonstrates the usage of TraceFilter + + disp('TraceFilter Cookbook') + disp(' - not created, yet -'); +end + diff --git a/classes/@TraceSpectra/TraceSpectra.m b/classes/@TraceSpectra/TraceSpectra.m index e9b120d..05289d1 100644 --- a/classes/@TraceSpectra/TraceSpectra.m +++ b/classes/@TraceSpectra/TraceSpectra.m @@ -1,16 +1,41 @@ classdef TraceSpectra - %TraceSpectra makes examining Trace spectra easy - % replaces spectralobject + %TraceSpectra Simple spectral analysis for SeismicTrace(s) + % + % most uses will require the Signal Processing Toolbox + % + % TraceSpectra properties: + % nfft - length of the fft + % over - amount of overlap (number of samples) + % freqmax - maximum frequency for display + % dBlims - decibel limits for display [lower upper] + % scaling - scaling factor (one of 's','m','h','d','date','doy') + % + % TraceSpectra methods: + % fft - Discrete Fourier transform. + % ifft - Inverse discrete Fourier transform. + % pwelch - Power Spectral Density estimate via Welch's method. (best to instead use psd, below) + % + % Newer (more powerful) representations: + % psd - Power Spectral Density, using newer matlab techniques + % spectrogram - spectrogram Spectrogram using a Short-Time Fourier Transform (STFT). + % + % display methods: + % specgram - plot the spectra + % specgram2 - plot the spectra and trace together + % + % s.specgram(trace) + % + % See also SeismicTrace, fft, ifft, pwelch, specgram properties nfft = 1024; % length of the fft - over = 1024 * 0.8; %a mount of overlap + over = 1024 * 0.8; % amount of overlap (number of samples) freqmax = 8; % maximum frequency for display - dBlims = [50 100]; % decibel limits + dBlims = [50 100]; % decibel limits for display [lower upper] (affects clim) scaling = 's'; % scaling factor (one of 's','m','h','d','date','doy') end properties(Dependent, Hidden) - map + map % colormap end properties(Hidden) SPECTRAL_MAP = TraceSpectra.loadmap; @@ -67,325 +92,6 @@ s.SPECTRAL_MAP = alternateMap; end %% plotting and spectra-taking functions - function h = specgram(s, ws, varargin) - %specgram - plots spectrogram of waveforms - % h = TraceSpectra.specgram(waveforms) Generates a spectrogram from the - % waveform(s), overwriting the current figure. The return value is a - % handle to the spectrogram, and is optional. The spectrograms will be - % created in the same shape as the passed waveforms. ie, if W is a 2x3 - % matrix of waveforms, then specgram(TraceSpectra,W) will generate a 2x3 - % plot of spectra. - % - % Many additional behaviors can be modified through the passing of - % additional parameters, as listed further below. These parameters are - % always passed in pairs, as such: - % - % TraceSpectra.specgram(waveforms,'PARAM1',VALUE1,...,'PARAMn',VALUEn) - % Any number of these parameters may be passed to specgram. - % - % specgram(..., 'axis', AXIS_HANDLE) - % Specify the axis AXIS_HANDLE within which the spectrogram will be - % generated. The boundary of the axis becomes the boundary for the - % entire spectra plot. For a matrix of waveforms, this area is - % subdivided into NxM subplots, where N and M are the size of the - % waveform matrix. - % - % specgram(..., 'xunit', XUNIT) - % Spedifies the x-unit scale to be used with the spectrogram. The - % default unit is 'seconds'. - % valid xunits: - % 'seconds','minutes','hours','days','doy' (day of year),and 'date' - % - % specgram(..., 'colormap', ALTERNATEMAP) - % Instead of using the default colormap, any colormap may be used. An - % alternate way of setting the global map is by using the SETMAP - % function. ALTERNATEMAP will either be a name (eg. grayscale) or an - % Nx3 numeric. Type HELP GRAPH3D to see additional useful colormaps. - % - % - % specgram(..., 'colorbar', COLORBAR_OPTION) - % Generates a spectrogram from the waveform and uses a specific map - % valid COLORBAR_OPTION values: 'horiz' (default),'vert','none', - % 'HORIZ' places a single colorbar below all plots - % 'VERT' places a single colorbar to the right of all plots - % 'NONE' supresses the colorbar placement - % - % specgram(..., 'yscale', YSCALE) - % Choosing 'log' Allows the y-axis to be generated on a log-frequency - % scale, with uneven vertical cell spacing. The default value is - % 'normal', and provides the standard spectrogram view. - % valid yscales: 'normal', 'log' (see NOTE below) - % - % NOTE: In order to use the log scale, UIMAGESC needs to be available on - % the matlab path. This routine was created by Frederic Moisy, and may - % be downloaded from the maltabcentral fileexchange (File ID: 11368). - % If this routine is not found,then the original spectrogram will be - % created. - % - % h = specgram(..., 'fontsize', FONTSIZE) - % Specify the font size for a spectrogram. The default font size is 8. - % - % specgram2(..., 'innerLabels', SHOWINNERLABELS) - % Suppress the labling of the inside graphs by setting SHOWINNERLABELS - % to false. If this is false, then the frequency label only shows on - % the leftmost spectrograms, and the X-unit label only shows on the - % bottommost spectrograms. - % - % The following plots a waveform using an alternate mapping, an xunit of - % 'hours', and with the y-axis plotted using a log scale. - % ex. specgram(TraceSpectra, waveform,... - % 'colormap', alternateMap,'xunit','h','yscale','log') - % - % - % Example: - % % create an arbitrary subplot, and then plot multiple spectra - % a = subplot(3,2,1); - % specgram(TraceSpectra,waves,'axis',a); % waves is an NxM waveform - % - % - % See also SPECTRALOBJECT/SPECGRAM2, WAVEFORM/PLOT, WAVEFORM/FILLGAPS - - - % Thanks to Jason Amundson for providing the way to do log scales - - currFontSize = 8; - %enforce input arguments. - if ~isa(ws,'TraceData') - error('Second input argument should be a TraceData, not a %s',class(ws)); - end - - hasExtraArg = mod(numel(varargin),2); - if hasExtraArg - error('Spectralobject:specgram:DepricatedArgument',... - ['%s/n%s/n%s','spectralobject/specgram now requires parameter pairs.',... - 'See help spectralobject/specgram for new usage instructions',... - 'most likely you tried to call specgram with an alternate ',... - 'color map without specifying the property ''colormap''']); - else - proplist= TraceSpectra.parseargs(varargin); - end - - %% search for relevent property pairs passed as parameters - - % AXIS: a handle to the axis, defining the area to be used - [~,myaxis,proplist] = TraceSpectra.getproperty('axis',proplist,0); - - % POSITION: 1x4 vector, specifying area in which to plot - % [left, bottom, width, height] - %[~,~,proplist] = TraceSpectra.getproperty('position',proplist,[]); - - % XUNIT: Specify the time units for the plot. eg, hours, minutes, doy, etc. - [~, xChoice, proplist] =... - TraceSpectra.getproperty('xunit',proplist,s.scaling); %'s' is default value for xChoice - - % FONTSIZE: specify the font size to be used for all labels within the plot - [~, currFontSize, proplist] =... - TraceSpectra.getproperty('fontsize',proplist,currFontSize);%default font size or override? - - %%take care of additional parsing that affects all waveforms - - [~, alternateMap, proplist] = ... - TraceSpectra.getproperty('colormap',proplist, s.SPECTRAL_MAP); %default map - - % YSCALE: either 'normal', or 'log' - [~, yscale, proplist] = ... - TraceSpectra.getproperty('yscale',proplist,'normal'); %default yscale to 'normal' - - logscale = strcmpi(yscale,'log'); - - % COLORBAR: Dictate the position of the colorbart relative to the plot - [~,colorbarpref,proplist] = TraceSpectra.getproperty('colorbar',proplist,'horiz'); - - % SUPRESSINNERLABELS: true, false - [~, suppressLabels, proplist] = ... - TraceSpectra.getproperty('innerlabels',proplist,false); %only show outside - - % SUPRESSXLABELS: true, false - [~, useXlabel, proplist] = ... - TraceSpectra.getproperty('useXlabel',proplist,true); %show no x, y labels - % SUPRESSXLABELS: true, false - [~, useYlabel, proplist] = ... - TraceSpectra.getproperty('useYlabel',proplist,true); %show no x, y labels - - - %% figure out exactly WHERE to plot the spectrogram(s) - %find out area(axis) in which the spectrograms will be plotted - clabel= 'Relative Amplitude (dB)'; - - if myaxis == 0, - clf; - myaxis = gca; - get(gca,'position'); - else - get(myaxis,'position'); - end - - %% If there are multiple waveforms... - % subdivide the axis and loop through specgram2 with individual waveforms. - - if numel(ws) > 1 - %create the colorbar if desired - TraceSpectra.createcolorbar(s, colorbarpref, clabel, currFontSize) - h = TraceSpectra.subdivide_axes(myaxis,size(ws)); - remainingproperties = TraceSpectra.property2varargin(proplist); - for n=1:numel(h) - keepYlabel = ~suppressLabels || (n <= size(h,1)); - keepXlabel = ~suppressLabels || (mod(n,size(h,2))==0); - specgram(s,ws(n),... - 'xunit',xChoice,... - 'axis',h(n),... - 'fontsize',currFontSize,... - 'useXlabel',keepXlabel,... - 'useYlabel',keepYlabel,... - 'colorbar','none',... - remainingproperties{:}); - end - return - end - - axes(myaxis); - - if any(isnan(double(ws))) - warning('Spectralobject:specgram:nanValue',... - ['This waveform has at least one NaN value, which will blank',... - 'the related spectrogram segment. ',... - 'Remove NaN values by either splitting up the ',... - 'waveform into non-NaN sections or by using waveform/fillgaps',... - ' to replace the NaN values.']); - end - - [xunit, xfactor] = TraceSpectra.parse_xunit(xChoice); - - switch lower(xunit) - case 'date' - % we need the actual times... - Xvalues = get(ws,'timevector'); - - case 'day of year' - startvec = datevec(get(ws,'start')); - Xvalues = get(ws,'timevector'); - dec31 = datenum(startvec(1)-1,12,31); % 12/31/xxxx of previous year in Matlab format - Xvalues = Xvalues - dec31; - xunit = [xunit, ' (', datestr(startvec,'yyyy'),')']; - - otherwise, - dl= 1:numel(ws.data); %dl : DataLength - Xvalues = dl .* ws.period ./ xfactor; - end - - - %% once was function specgram(d, NYQ, nfft, noverlap, freqmax, dBlims) - - nx = length(ws.data); - window = hanning(s.nfft); - nwind = length(window); - if nx < nwind % zero-pad x if it has length less than the window length - ws.data(nwind) = 0; - nx = nwind; - end - - ncol = fix( (nx - s.over) / (nwind - s.over) ); - - %added "floor" below - colindex = 1 + floor(0:(ncol-1))*(nwind- s.over); - rowindex = (1:nwind)'; - if length(ws.data)<(nwind+colindex(ncol)-1) - ws.data(nwind+colindex(ncol)-1) = 0; % zero-pad x - end - - y = zeros(nwind,ncol); - - % put x into columns of y with the proper offset - % should be able to do this with fancy indexing! - A_ = colindex( ones(nwind, 1) ,: ) ; - B_ = rowindex(:, ones(1, ncol) ) ; - y(:) = ws.data(fix(A_ + B_ -1)); - clear A_ B_ - - for k = 1:ncol; % remove the mean from each column of y - y(:,k) = y(:,k)-mean(y(:,k)); - end - - % Apply the window to the array of offset signal segments. - y = window(:,ones(1,ncol)).*y; - - % USE FFT - % now fft y which does the columns - y = fft(y,s.nfft); - if ~any(any(imag(ws.data))) % x purely real - if rem(s.nfft,2), % nfft odd - select = 1:(s.nfft+1)/2; - else - select = 1:s.nfft/2+1; - end - y = y(select,:); - else - select = 1:s.nfft; - end - f = (select - 1)' * ws.samplerate / s.nfft; - - NF = s.nfft/2 + 1; - nf1=round(f(1) / ws.nyquist * NF); %frequency window - if nf1==0, nf1=1; end - nf2=NF; - - y = 20*log10(abs(y(nf1:nf2,:))); - - F = f(f <= s.freqmax); - - - if F(1)==0, - F(1)=0.001; - end - - if logscale - t = linspace(Xvalues(1), Xvalues(end),ncol); % Replaced by TCB - 01/18/2012 - try - h = uimagesc(t, log10(F), y(nf1:length(F),:), s.dBlims); - catch exception - if strcmp(exception.identifier, 'MATLAB:UndefinedFunction') - warning('Spectralobject:specgram:uimageNotInstalled',... - ['Cannot plot with log spacing because uimage, uimagesc ',... - ' not installed or not visible in matlab path.']); - h = imagesc(Xvalues,F,y(nf1:length(F),:),s.dBlims); - logscale = false; - else - rethrow(exception) - end - end - else - h = imagesc(Xvalues, F, y(nf1:length(F),:), s.dBlims); - end - set(gca, 'fontsize', currFontSize); - if strcmpi(xunit, 'date') - datetick('x', 'keepticks'); - end - titlename = [ws.station '-' ws.channel ' from:' ws.start]; - title (titlename); - axis xy; - colormap(alternateMap); - shading flat - axis tight; - if useYlabel - if ~logscale - ylabel ('Frequency (Hz)') - else - ylabel ('Log Frequency (log Hz)') - end - - end - if useXlabel - xlabel(['Time - ',xunit]); - end - - %create the colorbar if desired - TraceSpectra.createcolorbar(s,colorbarpref, clabel, currFontSize); - - %% added a series of functions that help with argument parsing. - % These were ported from my waveform/plot function. - - - end %specgram %%located in another file: h = specgram(s, ws, varargin) % spectral plot @@ -454,7 +160,7 @@ end end function [varargout] = pwelch(s,w, varargin) - %PWELCH overloaded pwelch for spectralobjects + %PWELCH overloaded Power Spectral Density estimate via Welch's method. % pwelch(TraceSpectra, waveform) - plots the spectral density % Pxx = pwelch(TraceSpectra, waveform) - returns the Power Spectral % Density (PSD) estimate, Pxx, of a discrete-time signal @@ -469,7 +175,7 @@ % NOTE: voltage offsets may cause a large spike for lowest Pxx value. % NOTE: NaN values will result in blank % - % See also pwelch, waveform/fillgaps + % See also pwelch, TraceSpectra.fillgaps if nargin < 2 error('Spectralobject:pwelch:insufficientArguments',... @@ -477,17 +183,17 @@ end if ~isscalar(w) - error('Spectralobject:pwelch:nonScalarWaveform',... + error('TraceSpectra:pwelch:nonScalarWaveform',... 'waveform must be scalar (1x1)'); end - if ~isa(w,'waveform') - error('Spectralobject:pwelch:invalidArgument',... - 'second argument expected to be WAVEFORM, but was [%s]', class(w)); + if ~isa(w,'TraceData') + error('TraceSpectra:pwelch:invalidArgument',... + 'second argument expected to be trace, but was [%s]', class(w)); end - if any(isnan(double(w))) - warning('Spectralobject:pwelch:nanValue',... + if any([w.hasnan]) + warning('TraceSpectra:pwelch:nanValue',... ['This waveform has at least one NaN value. ',... 'Remove NaN values by either splitting up the',... ' waveform into non-NaN sections or by using ',... @@ -496,12 +202,92 @@ if nargin == 3 if strcmpi(varargin{1}, 'DEFAULT') %disp('defaulting') - window = []; - over = []; + s.window = []; + s.over = []; end end [varargout{1:nargout}] = pwelch(double(w),numel(w.data),s.over,s.nfft,w.samplerate); end + function [varargout] = spectrogram(s, w, varargin) + %spectrogram spectrogram Spectrogram using a Short-Time Fourier Transform (STFT). + % run(s, trace) will attempt to run any matlab spectral function + % of the form fn(X, window, noverlap, nfft, fs) with a default + % window of length nfft + % + % example: + % s = TraceSpectra; + % T = SeismicTrace.synthetic() + % spectrogram(s, T) + % + % See also: spectrogram, periodogram + assert(numel(w) == 1, 'This function requires a single SeismicTrace'); + switch nargout + case 0 + spectrogram(w.data, [], s.nfft, round(s.over), w.samplerate); + case 1 + a = spectrogram(w.data, [], s.nfft, round(s.over), w.samplerate); + varargout = {a}; + case 2 + [a, b] = spectrogram(w.data, [], s.nfft, round(s.over), w.samplerate); + varargout = {a,b}; + case 3 + [a, b, c] = spectrogram(w.data, [], s.nfft, round(s.over), w.samplerate); + varargout = {a,b,c}; + otherwise + error('TraceSpectra:run:unexpectedNumberOfOutputs',... + 'Unanticipated number of outputs. expected between 0 and 3'); + end + end + function [varargout] = psd(s, w, estimator, varargin) + %psd Power Spectral Density, using newer matlab techniques + % psd(s, trace, estimator) will attempt to run the Power Spectral Density function + % of the form fn(X, window, noverlap, nfft, fs) with a default + % window of length nfft + % + % only works on a single trace at a time + % + % a list of estimators can be found by looking at the help for spectrum.psd + % + % Valid psd estimators (as of 2012b:): + % periodogram mcov + % welch mtm + % burg yulear + % cov + % + % example: + % s = TraceSpectra; + % T = SeismicTrace.synthetic() + % estimator = spectrum.pwelch; %choose a valid estimator + % s.psd(T, estimator); % plot + % Hpsd = s.psd(T, estimator); % return values in a psd object + % + % See also: spectrum.psd, spectrum.welch, spectrum.burg, + % spectrum.cov, spectrum.mcov, spectrum.mtm, spectrum.yulear, + % spectrum.periodogram + + assert(numel(w) == 1, 'This function requires a single SeismicTrace'); + + % suppressed the warning about psd being oudated (FDEPR) because I am + % using the new one. Because estimator is a variable, mlint can't + % tell the difference. + + switch nargout + case 0 + psd(estimator, w.data, 'fs', w.samplerate,'nfft',s.nfft); %#ok<*FDEPR> + case 1 + a = psd(estimator, w.data, 'fs', w.samplerate,'nfft',s.nfft); + varargout = {a}; + case 2 + [a, b] = psd(estimator, w.data, 'fs', w.samplerate,'nfft',s.nfft); + varargout = {a,b}; + case 3 + [a, b, c] = psd(estimator, w.data, 'fs', w.samplerate,'nfft',s.nfft); + varargout = {a,b,c}; + otherwise + error('TraceSpectra:run:unexpectedNumberOfOutputs',... + 'Unanticipated number of outputs. expected between 0 and 3'); + end + end function handle = colorbar_axis(s,loc,clabel,rlab1,rlab2, fontsize) % COLORBAR - Display color bar (color scale). @@ -549,8 +335,6 @@ end if nargin>1, - lower=s.dBlims(1); - upper=s.dBlims(2); paxis = s.dBlims; end @@ -592,8 +376,8 @@ end else t = paxis; - cmin=t(1); - cmax=t(2); + % cmin=t(1); %apparently not used + % cmax=t(2); %apparently not used end h = gca; @@ -644,7 +428,7 @@ % Create color stripe n = size(colormap,1); - image([0 1],t,[1:n]'); set(ax,'Ydir','normal') + image([0 1],t,(1:n)'); set(ax,'Ydir','normal') if nargin>2, set(ax,'ylabel',text(0,0,clabel, 'FontSize', fontsize)); @@ -664,7 +448,8 @@ end end %d=date; - yshift=.12*(ypos(2)-ypos(1)); + + % yshift=.12*(ypos(2)-ypos(1)); % text(xpos(1)+2.0,ypos(1)-yshift,0.,d) % Create color axis @@ -679,19 +464,29 @@ labels = []; width = []; set (gca, 'FontSize', fontsize) %This sets the font size - for i=1:length(yticks), - labels = [labels;text(1+0*xspace,yticks(i),deblank(ylabels(i,:)), ... + for i=length(yticks):-1:1 + % unknown length for each, so cannot preallocat + labels(i) = text(1+0*xspace,yticks(i),deblank(ylabels(i,:)), ... + 'HorizontalAlignment','right', ... + 'VerticalAlignment','middle', ... + 'FontName',get(ax,'FontName'), ... + 'FontSize',get(ax,'FontSize'), ... + 'FontAngle',get(ax,'FontAngle'), ... + 'FontWeight',get(ax,'FontWeight')); + %{ + labels(i) = [text(1+0*xspace,yticks(i),deblank(ylabels(i,:)), ... 'HorizontalAlignment','right', ... 'VerticalAlignment','middle', ... 'FontName',get(ax,'FontName'), ... 'FontSize',get(ax,'FontSize'), ... 'FontAngle',get(ax,'FontAngle'), ... 'FontWeight',get(ax,'FontWeight'))]; - width = [width;get(labels(i),'Extent')]; + %} + width((i*4 - 3):i*4) = get(labels(i),'Extent'); end % Shift labels over so that they line up - [dum,k] = max(width(:,3)); width = width(k,3); + [~,k] = max(width(:,3)); width = width(k,3); for i=1:length(labels), pos = get(labels(i),'Position'); set(labels(i),'Position',[pos(1)+width pos(2:3)]) @@ -699,7 +494,7 @@ % If we need an exponent then draw one [ymax,k] = max(abs(yticks)); - if abs(abs(str2num(ylabels(k,:)))-ymax)>sqrt(eps), + if abs(abs(str2num(ylabels(k,:)))-ymax)>sqrt(eps), %#ok ex = log10(max(abs(yticks))); ex = sign(ex)*ceil(abs(ex)); l = text(0,ylim(2)+2*yspace,'x 10', ... @@ -726,12 +521,13 @@ if isempty(ax), pos = get(h,'Position'); %[left,bottom,width,height] stripe = 0.05; space = 0.1; %stripe = 0.075 - ori=get(gcf,'paperorientation'); if fsize<=16, sfact=1; else sfact=2; end + + % ori=get(gcf,'paperorientation'); % if ori=='landscape', stripe = 0.05; end % if fsize<=26, % set(h,'Position',... @@ -770,7 +566,7 @@ end % image(t,[0 1],[1:n]); set(ax,'Ydir','normal') - image(t,[0 1],[n1:n2]); set(ax,'Ydir','normal') + image(t,[0 1],n1:n2); set(ax,'Ydir','normal') set(ax,'yticklabelmode','manual') set(ax,'yticklabel','') @@ -791,8 +587,9 @@ text(xpos(2)-xshift,ypos(1)-1.0,0.,rlab2); end end - d=date; - xshift=.15*(xpos(2)-xpos(1)); + + % d=date; + % xshift=.15*(xpos(2)-xpos(1)); % text(xpos(1)-xshift,ypos(1)-3.5,0.,d) %set(gca,'xtick',[]) @@ -806,7 +603,11 @@ if nargout>0, handle = ax; end end + end %methods + methods(Access=protected) + p = parseSpecgramInputs(me, cellOfArgs) + end methods(Static, Access=protected) function [unitName, secondMultiplier] = parse_xunit(unitName) % PARSE_XUNIT returns a labelname and a multiplier for an incoming xunit @@ -1010,6 +811,20 @@ set(hbar,'fontsize',currFontSize) end end + + function newParams = buildParameterList(paramStruct) + %newParams creates argument list from a struct + fieldList = fieldnames(paramStruct); + if isempty(fieldList) + newParams = {}; + else + v = struct2cell(paramStruct); + newParams = horzcat(fieldList,v); + end + end + end + methods(Static) + cookbook(); % cookbook for TraceSpectra end end diff --git a/classes/@TraceSpectra/cookbook.m b/classes/@TraceSpectra/cookbook.m new file mode 100644 index 0000000..4f51bd3 --- /dev/null +++ b/classes/@TraceSpectra/cookbook.m @@ -0,0 +1,8 @@ +function cookbook() + %cookbook Cookbook for TraceSpectra + % Demonstrates the usage of TraceSpectra + + disp('TraceSpectra Cookbook') + disp(' - not created, yet -'); +end + diff --git a/classes/@TraceSpectra/parseSpecgramInputs.m b/classes/@TraceSpectra/parseSpecgramInputs.m new file mode 100644 index 0000000..0e8678f --- /dev/null +++ b/classes/@TraceSpectra/parseSpecgramInputs.m @@ -0,0 +1,38 @@ +function p = parseSpecgramInputs(me, cellOfArgs) + %parseSpecgramInputs handles inputs for specgram functions + % p = parseSpecgramInputs(obj, cell_of_arguments) will parse the + % cell of arguments into an inputParser. specgram-specific + % arguments will be returned in p.Results, while any additional + % arguments will be returned in p.Unmatched. + % + % specifically looks for: + % 'axis' - a handle to the axis, defining the area to be used + % 'position' - 1x4 vector, specifying area in which to plot [left, bottom, width, height] + % 'colorbar' - Position of the colorbar relative to the plot + % 'colormap' - alternate colormap + % 'xunit' - Specify the time units for the plot. eg, hours, minutes, doy, etc. + % 'fontsize' - specify the font size to be used for all labels within the plot + % 'yscale' - either 'normal', or 'log' + % 'innerlabels' - true/false + % 'useXlabel' - true/false + % 'useYlabel' - true/false + p = inputParser; + p.StructExpand = true; + p.KeepUnmatched = true; + p.CaseSensitive = false; + if verLessThan('matlab','8.2') %2013b + addParameter = @addParamValue; % usage is the same as of 2015b + end + %NOTE: older version of matlab requires addParamValue instead of addParamter + addParameter(p,'axis', 0); %myaxis + addParameter(p,'position', []); + addParameter(p,'colorbar', 'horiz'); + addParameter(p,'colormap', me.SPECTRAL_MAP); + addParameter(p,'xunit', me.scaling); + addParameter(p,'fontsize', 10); + addParameter(p,'yscale', 'normal'); + addParameter(p,'innerlabels', false); + addParameter(p,'useXlabel', true); + addParameter(p,'useYlabel', true); + p.parse(cellOfArgs{:}); +end diff --git a/classes/@TraceSpectra/specgram.m b/classes/@TraceSpectra/specgram.m index fd1a49d..f04b7bb 100644 --- a/classes/@TraceSpectra/specgram.m +++ b/classes/@TraceSpectra/specgram.m @@ -1,25 +1,25 @@ - function h = specgram(s, ws, varargin) - %specgram - plots spectrogram of waveforms - % h = TraceSpectra.specgram(waveforms) Generates a spectrogram from the - % waveform(s), overwriting the current figure. The return value is a + function h = specgram(s, T, varargin) + %specgram - plots spectrogram of traces + % h = TraceSpectra.specgram(traces) Generates a spectrogram from the + % trace(s), overwriting the current figure. The return value is a % handle to the spectrogram, and is optional. The spectrograms will be - % created in the same shape as the passed waveforms. ie, if W is a 2x3 - % matrix of waveforms, then specgram(TraceSpectra,W) will generate a 2x3 + % created in the same shape as the passed traces. ie, if W is a 2x3 + % matrix of traces, then specgram(TraceSpectra,W) will generate a 2x3 % plot of spectra. % % Many additional behaviors can be modified through the passing of % additional parameters, as listed further below. These parameters are % always passed in pairs, as such: % - % TraceSpectra.specgram(waveforms,'PARAM1',VALUE1,...,'PARAMn',VALUEn) + % TraceSpectra.specgram(traces,'PARAM1',VALUE1,...,'PARAMn',VALUEn) % Any number of these parameters may be passed to specgram. % % specgram(..., 'axis', AXIS_HANDLE) % Specify the axis AXIS_HANDLE within which the spectrogram will be % generated. The boundary of the axis becomes the boundary for the - % entire spectra plot. For a matrix of waveforms, this area is + % entire spectra plot. For a matrix of traces, this area is % subdivided into NxM subplots, where N and M are the size of the - % waveform matrix. + % trace matrix. % % specgram(..., 'xunit', XUNIT) % Spedifies the x-unit scale to be used with the spectrogram. The @@ -35,7 +35,7 @@ % % % specgram(..., 'colorbar', COLORBAR_OPTION) - % Generates a spectrogram from the waveform and uses a specific map + % Generates a spectrogram from the trace and uses a specific map % valid COLORBAR_OPTION values: 'horiz' (default),'vert','none', % 'HORIZ' places a single colorbar below all plots % 'VERT' places a single colorbar to the right of all plots @@ -62,30 +62,34 @@ % the leftmost spectrograms, and the X-unit label only shows on the % bottommost spectrograms. % - % The following plots a waveform using an alternate mapping, an xunit of + % The following plots a trace using an alternate mapping, an xunit of % 'hours', and with the y-axis plotted using a log scale. - % ex. specgram(TraceSpectra, waveform,... + % ex. specgram(TraceSpectra, trace,... % 'colormap', alternateMap,'xunit','h','yscale','log') % % % Example: % % create an arbitrary subplot, and then plot multiple spectra % a = subplot(3,2,1); - % specgram(TraceSpectra,waves,'axis',a); % waves is an NxM waveform + % specgram(TraceSpectra,waves,'axis',a); % waves is an NxM trace % % - % See also SPECTRALOBJECT/SPECGRAM2, WAVEFORM/PLOT, WAVEFORM/FILLGAPS + % See also SPECTRALOBJECT/SPECGRAM2, trace/PLOT, trace/FILLGAPS % Thanks to Jason Amundson for providing the way to do log scales - currFontSize = 8; + % currFontSize = 8; % replaced in the parser %enforce input arguments. - if ~isa(ws,'TraceData') - error('Second input argument should be a TraceData, not a %s',class(ws)); + if ~isa(T,'TraceData') + try + T = SeismicTrace(T); + catch er + error('Second input argument should be a TraceData, not a %s',class(T)); + end end - p = parseInputs(s, varargin); + p = parseSpecgramInputs(s, varargin); %% search for relevent property pairs passed as parameters myaxis = p.Results.axis; % area to be used @@ -113,18 +117,18 @@ get(myaxis,'position'); end - %% If there are multiple waveforms... - % subdivide the axis and loop through specgram2 with individual waveforms. + %% If there are multiple traces... + % subdivide the axis and loop through specgram2 with individual traces. - if numel(ws) > 1 + if numel(T) > 1 %create the colorbar if desired TraceSpectra.createcolorbar(s, colorbarpref, clabel, currFontSize) - h = TraceSpectra.subdivide_axes(myaxis,size(ws)); - remainingproperties = buildParameterList(p.Unmatched); + h = TraceSpectra.subdivide_axes(myaxis,size(T)); + remainingproperties = TraceSpectra.buildParameterList(p.Unmatched); for n=1:numel(h) keepYlabel = ~suppressLabels || (n <= size(h,1)); keepXlabel = ~suppressLabels || (mod(n,size(h,2))==0); - specgram(s,ws(n),... + specgram(s,T(n),... 'xunit',xChoice,... 'axis',h(n),... 'fontsize',currFontSize,... @@ -138,12 +142,12 @@ axes(myaxis); - if ws.hasnan % only 1 trace guaranteed at this point... + if T.hasnan % only 1 trace guaranteed at this point... warning('TraceSpectra:specgram:nanValue',... ['This trace has at least one NaN value, which will blank',... 'the related spectrogram segment. ',... 'Remove NaN values by either splitting up the ',... - 'waveform into non-NaN sections or by using TraceData/fillgaps',... + 'trace into non-NaN sections or by using TraceData/fillgaps',... ' to replace the NaN values.']); end @@ -152,28 +156,28 @@ switch lower(xunit) case 'date' % we need the actual times... - Xvalues = get(ws,'timevector'); + Xvalues = get(T,'timevector'); case 'day of year' - startvec = datevec(get(ws,'start')); - Xvalues = get(ws,'timevector'); + startvec = datevec(get(T,'start')); + Xvalues = get(T,'timevector'); dec31 = datenum(startvec(1)-1,12,31); % 12/31/xxxx of previous year in Matlab format Xvalues = Xvalues - dec31; xunit = [xunit, ' (', datestr(startvec,'yyyy'),')']; otherwise, - dl= 1:ws:nsamples(); %dl : DataLength - Xvalues = dl .* ws.period ./ xfactor; + dl= 1:T.nsamples(); %dl : DataLength + Xvalues = dl .* T.period ./ xfactor; end %% once was function specgram(d, NYQ, nfft, noverlap, freqmax, dBlims) - nx = length(ws.data); + nx = length(T.data); window = hanning(s.nfft); nwind = length(window); if nx < nwind % zero-pad x if it has length less than the window length - ws.data(nwind) = 0; + T.data(nwind) = 0; nx = nwind; end @@ -182,8 +186,8 @@ %added "floor" below colindex = 1 + floor(0:(ncol-1))*(nwind- s.over); rowindex = (1:nwind)'; - if length(ws.data)<(nwind+colindex(ncol)-1) - ws.data(nwind+colindex(ncol)-1) = 0; % zero-pad x + if length(T.data)<(nwind+colindex(ncol)-1) + T.data(nwind+colindex(ncol)-1) = 0; % zero-pad x end y = zeros(nwind,ncol); @@ -192,7 +196,7 @@ % should be able to do this with fancy indexing! A_ = colindex( ones(nwind, 1) ,: ) ; B_ = rowindex(:, ones(1, ncol) ) ; - y(:) = ws.data(fix(A_ + B_ -1)); + y(:) = T.data(fix(A_ + B_ -1)); clear A_ B_ for k = 1:ncol; % remove the mean from each column of y @@ -205,7 +209,7 @@ % USE FFT % now fft y which does the columns y = fft(y,s.nfft); - if ~any(any(imag(ws.data))) % x purely real + if ~any(any(imag(T.data))) % x purely real if rem(s.nfft,2), % nfft odd select = 1:(s.nfft+1)/2; else @@ -215,10 +219,10 @@ else select = 1:s.nfft; end - f = (select - 1)' * ws.samplerate / s.nfft; + f = (select - 1)' * T.samplerate / s.nfft; NF = s.nfft/2 + 1; - nf1=round(f(1) / ws.nyquist * NF); %frequency window + nf1=round(f(1) / T.nyquist * NF); %frequency window if nf1==0, nf1=1; end nf2=NF; @@ -253,7 +257,7 @@ if strcmpi(xunit, 'date') datetick('x', 'keepticks'); end - titlename = [ws.station '-' ws.channel ' from:' ws.start]; + titlename = [T.station '-' T.channel ' from:' T.start]; title (titlename); axis xy; colormap(alternateMap); @@ -275,49 +279,8 @@ TraceSpectra.createcolorbar(s,colorbarpref, clabel, currFontSize); %% added a series of functions that help with argument parsing. - % These were ported from my waveform/plot function. + % These were ported from my trace/plot function. end %specgram - -function p = parseInputs(me, cellOfArgs) - %TODO: Unify(?) specgram2 and specgram parseInputs - p = inputParser; - p.StructExpand = true; - p.KeepUnmatched = true; - p.CaseSensitive = false; - %NOTE: older version of matlab requires addParamValue instead of addParamter - % AXIS: a handle to the axis, defining the area to be used - p.addParameter('axis', 0); %myaxis - % POSITION: 1x4 vector, specifying area in which to plot [left, bottom, width, height] - p.addParameter('position', []); - % alternateMap: COLORMAP: - p.addParameter('colormap', me.SPECTRAL_MAP); - % COLORBAR: Dictate the position of the colorbart relative to the plot - p.addParameter('colorbar', 'horiz'); - % XUNIT: Specify the time units for the plot. eg, hours, minutes, doy, etc. - p.addParameter('xunit', me.scaling); - % FONTSIZE: specify the font size to be used for all labels within the plot - p.addParameter('fontsize', 10); - % YSCALE: either 'normal', or 'log' - p.addParameter('yscale', 'normal'); - % SUPRESSINNERLABELS: true, false - p.addParameter('innerlabels', false); - % USEXLABELS: true, false - p.addParameter('useXlabel', true); - % USEYLABELS: true, false - p.addParameter('useYlabel', true); - p.parse(cellOfArgs{:}); -end -function newParams = buildParameterList(paramStruct) - %newParams creates argument list from a struct - %TODO: Unify(?) specgram2 and specgram buildParameterList - fieldList = fieldnames(paramStruct); - if isempty(fieldList) - newParams = {}; - else - v = struct2cell(paramStruct); - newParams = horzcat(fieldList,v); - end -end \ No newline at end of file diff --git a/classes/@TraceSpectra/specgram2.m b/classes/@TraceSpectra/specgram2.m index 1a484aa..a355dc9 100644 --- a/classes/@TraceSpectra/specgram2.m +++ b/classes/@TraceSpectra/specgram2.m @@ -79,11 +79,16 @@ % See also SPECTRALOBJECT/SPECGRAM if ~isa(T,'TraceData') - error('Spectralobject:specgram2:invalidArgument','Second input argument should be a trace (ex. TraceData, SeismicTrace), not a %s',class(T)); + try + T = SeismicTrace(T); + disp('successfully converted to a SeismicTrace'); + catch er + error('Spectralobject:specgram2:invalidArgument','Should work on a trace (ex. TraceData, SeismicTrace), not a %s',class(T)); + end end %% search for relevent property pairs passed as parameters - p = parseInputs(s, varargin); + p = parseSpecgramInputs(s, varargin); %% figure out exactly WHERE to plot the spectrogram(s) %find out area(axis) in which the spectrograms will be plotted @@ -108,7 +113,8 @@ %create the colorbar if desired TraceSpectra.createcolorbar(s,p.Results.colorbar, clabel, p.Results.fontsize); h = TraceSpectra.subdivide_axes(myaxis, size(T)); - remainingproperties = TraceSpectra.property2varargin(proplist); + % remainingproperties = TraceSpectra.property2varargin(proplist); + remainingproperties = TraceSpectra.buildParameterList( p.Unmatched); for n=1:numel(h) keepYlabel = ~p.Results.innerlabels || (n <= size(h,1)); keepXlabel = ~p.Results.innerlabels || (mod(n,size(h,2))==0); @@ -120,7 +126,6 @@ 'useYlabel',keepYlabel,... 'colorbar','none',... remainingproperties{:}); - end return end @@ -138,7 +143,7 @@ %plot the spectra ax_spectra = subplot('position',spectraPosition(pos)); - additionalParams = buildParameterList(p.Unmatched); + additionalParams = TraceSpectra.buildParameterList( p.Unmatched); specgram(s,T,... 'xunit',p.Results.xunit,... 'fontsize',p.Results.fontsize,... @@ -173,42 +178,3 @@ %spectraPosition spectra occupies the bottom 85% of the axis specPos = pos .* [1, 1, 1, 0.85 ]; end - - -function p = parseInputs(me, cellOfArgs) - p = inputParser; - p.StructExpand = true; - p.KeepUnmatched = true; - p.CaseSensitive = false; - %NOTE: older version of matlab requires addParamValue instead of addParamter - % AXIS: a handle to the axis, defining the area to be used - p.addParameter('axis', 0); %myaxis - % POSITION: 1x4 vector, specifying area in which to plot [left, bottom, width, height] - p.addParameter('position', []); - % COLORBAR: Dictate the position of the colorbart relative to the plot - p.addParameter('colorbar', 'horiz'); - % XUNIT: Specify the time units for the plot. eg, hours, minutes, doy, etc. - p.addParameter('xunit', me.scaling); - % FONTSIZE: specify the font size to be used for all labels within the plot - p.addParameter('fontsize', 10); - % YSCALE: either 'normal', or 'log' - p.addParameter('yscale', 'normal'); - % SUPRESSINNERLABELS: true, false - p.addParameter('innerlabels', false); - % USEXLABELS: true, false - p.addParameter('useXlabel', true); - % USEYLABELS: true, false - p.addParameter('useYlabel', true); - p.parse(cellOfArgs{:}); -end - -function newParams = buildParameterList(paramStruct) - %newParams creates argument list from a struct - fieldList = fieldnames(paramStruct); - if isempty(fieldList) - newParams = {}; - else - v = struct2cell(paramStruct); - newParams = horzcat(fieldList,v); - end -end \ No newline at end of file