added mean and peak frequency option to spectrograms

geoscience-community-codes · Jun 24, 2016 · 8582f5b · 8582f5b
1 parent cbe641b
commit 8582f5b
Show file tree

Hide file tree

Showing 7 changed files with 138 additions and 73 deletions.
diff --git a/applications/+iceweb/spectrogram_iceweb.m b/applications/+iceweb/spectrogram_iceweb.m
@@ -63,7 +63,7 @@
 
 nfft = round(get(s,'nfft'));
 overlap = floor(get(s, 'over'));
-dBlims = get(s, 'dBlims'); 
+dBlims = get(s, 'dBlims')
 fmax = get(s, 'freqmax');
 
 % Set appropriate date ticks
@@ -81,7 +81,9 @@
         [S,F,T] = spectrogram(data, nfft, nfft/2, nfft, fsamp);
 
         Y = 20*log10(abs(S)+eps);
-        index = find(F <= fmax); 
+        fmax
+        max(F)
+        index = find(F <= fmax)
         if F(1)==0,
             F(1)=0.001;
         end
@@ -91,7 +93,24 @@
         T = wt.start + T/86400;
         F = F(1:max(index));
         Y = Y(1:max(index),:);
-        imagesc(T,F,Y,dBlims); 
+        S = S(1:max(index),:);
+        if isempty(dBlims)
+            imagesc(T,F,abs(S));
+            % mean frequency
+            numerator = abs(S)' * F;
+            denominator = sum(abs(S),1);
+            meanF = numerator./denominator';
+            hold on; plot(T,meanF,'k','LineWidth',3);
+            % peak frequency
+            [maxvalue,maxindex] = max(abs(S));
+            size(maxindex)
+            fmax = F(maxindex);
+            size(fmax)
+            hold on; plot(T,fmax,'r','LineWidth',3);
+
+        else
+            imagesc(T,F,Y,dBlims); 
+        end
         axis xy;
         colormap(mycolormap);
 

diff --git a/core/@spectralobject/spectralobject.m b/core/@spectralobject/spectralobject.m
@@ -54,6 +54,8 @@
 
         if ~isempty(dBlims)
             s = set(s,'dblims',dBlims);
+        else
+            s.dBlims = [];
         end
 
         if exist('scaling','var')

diff --git a/core/@waveform/get.m b/core/@waveform/get.m
@@ -225,7 +225,8 @@
 %%%%%%%%%%%%%%%%
 function val = grabEndTime(w)
    if isempty(w)
-      val = [];
+      %val = [];
+      val = w.start;
       return
    end
    dlens = get(w,'data_length');

diff --git a/core/@waveform/pad.m b/core/@waveform/pad.m
@@ -10,14 +10,16 @@
 %       w2 = pad(w, min(snum), max(enum), 0)
 
 % Glenn Thompson USF
-
     if ~exist('padvalue','var')
         padvalue=0; % could also be nanmean?
     end
 
 	for waveform_num=1:numel(w)
 
         thisw = w(waveform_num);
+        if isempty(thisw)
+            continue;
+        end
 		y = get(thisw,'data');
 		fs = get(thisw,'freq');
 		[wsnum wenum] = gettimerange(thisw);

diff --git a/core/@waveform/private/load_antelope.m b/core/@waveform/private/load_antelope.m
@@ -185,8 +185,10 @@
       end
 
       % Glenn 2016/05/12 Get the trace objects
-      [wcell, database, fdb] = get_traces(startDatenums, endDatenums, expr, database, combineWaves);
-      w=wcell{:};
+      [w, database, fdb] = get_traces(startDatenums, endDatenums, expr, database, combineWaves);
+      if iscell(w)
+          w=w{:};
+      end
 
       % Glenn 2016/05/12 Close db if open
       if fdb.database ~= -102   % fdb ~= dbinvalid
@@ -277,11 +279,11 @@
 end
 
 %%
-function [w, rawDb, filteredDb] =  get_traces(startDatenums, endDatenums, expr, database, combineWaves)
+function [w, rawDb, filteredDb] =  get_traces(startDatenum, endDatenum, expr, database, combineWaves)
    % GET_TRACES gets interesting data from database, and returns the tracebuf object
-   % [tr(1:end)] =  get_antelope_trace(startDatenums, endDatenums, criteriaList, database)
-   %    STARTDATE is a matlab datenum
-   %    ENDDATE is also in matlab datenum
+   % [tr(1:end)] =  get_antelope_trace(startDatenum, endDatenum, criteriaList, database)
+   %    STARTDATENUM is a matlab datenum
+   %    ENDDATENUM is also in matlab datenum
    %    EXPR is the expression list formed from the station & channel
    %    combinations
    %    DATABASE can either be a database name, or an open antelope database
@@ -335,8 +337,8 @@
    % do not close the database if a pointer is asked for in the return arguments
    onFinishCloseDB = nargout < 2;
 
-   start_epochs = mep2dep(startDatenums);
-   end_epochs = mep2dep(endDatenums);
+   start_epoch = mep2dep(startDatenum);
+   end_epoch = mep2dep(endDatenum);
 
    %if the database isn't already open, then open it for reading
    if ~useExistingDatabasePtr
@@ -367,46 +369,48 @@
       cleanUpFail('Waveform:load_antelope:dataNotFound', 'No records found for criteria [%s].', expr);
       return;
    end
-
+   debug.print_debug(1, sprintf('Found %d matching wfdisc records after sta/chan subset', safe_dbnrecs(mydb)));
+
    filteredDb = mydb;
 
    [wfdisctime, wfdiscendtime] = dbgetv(mydb,'time','endtime');
+
+
    %% Get the tracebuf object for this starttime, endtime
-   % Loop through all times.  Result is tr(1:numel(startDatenums) of all tracebuffers.
-   for mytimeIDX = 1:numel(start_epochs)
-      someDataExists = any(start_epochs(mytimeIDX)<= wfdiscendtime & end_epochs(mytimeIDX) >= wfdisctime);
-
-      if someDataExists
-
+    someDataExists = any(start_epoch<= wfdiscendtime & end_epoch >= wfdisctime);
+
+    if someDataExists
+
          % time subset
          nbefore = safe_dbnrecs(mydb);
-         expr_time = sprintf('time < %f  && endtime > %f ',end_epochs(mytimeIDX), start_epochs(mytimeIDX));
+         expr_time = sprintf('time < %f  && endtime > %f ',end_epoch, start_epoch);
          dbptr = dbsubset(mydb, expr_time);
+         debug.print_debug(1, sprintf('Found %d matching wfdisc records after time subset', safe_dbnrecs(mydb)));
+
 
          % Display matching records
          if debug.get_debug()>0
              nnow = safe_dbnrecs(dbptr);
-            fprintf('Found %d matching time records (had %d before time subset):',nnow,nbefore);
-            for c=1:safe_dbnrecs(dbptr)
-                dbptr.record = c-1;
-                [wfid, sta, chan, wftime, wfendtime]=dbgetv(dbptr , 'wfid','sta','chan','time', 'endtime');
-                fprintf('%12d %s %s %s %s\n',wfid, sta, chan, datestr(epoch2datenum(wftime)), datestr(epoch2datenum(wfendtime)));   
-            end
+             fprintf('Found %d matching time records (had %d before time subset):',nnow,nbefore);
+             for c=1:safe_dbnrecs(dbptr)                
+                    dbptr.record = c-1;                   
+                    [wfid, sta, chan, wftime, wfendtime]=dbgetv(dbptr , 'wfid','sta','chan','time', 'endtime');                   
+                    fprintf('%12d %s %s %s %s\n',wfid, sta, chan, datestr(epoch2datenum(wftime)), datestr(epoch2datenum(wfendtime)));             
+             end
          end
 
          %%% Glenn Thompson 2012/02/06: Occasionally the C program trload_css cannot even load the trace data.
          % This error needs to be handled. So adding a try..catch..end around the original instruction.
          try
-            debug.print_debug(1,sprintf('\nTRYING TO LOAD TRACES FROM WHOLE WFDISC SUBSET IN ONE GO\n'));
-            tr = trload_css(dbptr, start_epochs(mytimeIDX), end_epochs(mytimeIDX));
-            trsplice(tr,20);            
-            %st = tr2struct(tr)
-            w0 = trace2waveform(tr);  % Glenn's simple method
-            %w0 = traceToWaveform(tr);  % Celso's more sophisticated method
-            %w0 = cycleThroughTraces(tr, combineWaves); % This is even more
-            %sophisticated, and calls traceToWaveform and others, but
-            %causes seg faults.
-            trdestroy( tr );
+                debug.print_debug(1,sprintf('\nTRYING TO LOAD TRACES FROM WHOLE WFDISC SUBSET IN ONE GO\n'));
+                tr = trload_css(dbptr, start_epoch, end_epoch);
+                trsplice(tr,20);            
+                w0 = trace2waveform(tr)  % Glenn's simple method
+                %w0 = traceToWaveform(tr);  % Celso's more sophisticated method
+                %w0 = cycleThroughTraces(tr, combineWaves); % This is even more
+                %sophisticated, and calls traceToWaveform and others, but
+                %causes seg faults.
+                trdestroy( tr );
          catch
                 debug.print_debug(1,sprintf('\nBulk mode failed\nTRYING TO LOAD TRACES FROM ONE WFDISC ROW AT A TIME\n'));
                 w0 = [];
@@ -418,7 +422,7 @@
                     [wfid, sta, chan, st, et]=dbgetv(dbptr_one_record , 'wfid','sta','chan','time', 'endtime');
                     debug.print_debug(1,sprintf('%12d %s %s %f %f\n',wfid, sta, chan, st, et));
                     try
-                        tr = trload_css(dbptr_one_record , start_epochs(mytimeIDX), end_epochs(mytimeIDX));
+                        tr = trload_css(dbptr_one_record , start_epoch, end_epoch);
                         trsplice(tr,20);
                         w0 = [w0;trace2waveform(tr)];  % Glenn's simple method
                         %w0 = [w0;traceToWaveform(tr)];  % Celso's more sophisticated method
@@ -429,7 +433,7 @@
                     catch ME
                         if strcmp(ME.identifier, 'MATLAB:unassignedOutputs')
                             % no trace table returned by trload_css
-                            w0 = [w0; waveform(ChannelTag('',sta,'',chan), NaN, epoch2datenum(start_epochs(mytimeIDX)), [], '')];
+                            w0 = [w0; waveform(ChannelTag('',sta,'',chan), NaN, epoch2datenum(start_epoch), [], '')];
                         else
                             rethrow(ME)
                         end
@@ -442,14 +446,15 @@
                     end     
                 end             
          end
-      else
-         w0=waveform(); 
-      end
-      w{mytimeIDX} = combine(w0);
-      %w{mytimeIDX} = pad(w{mytimeIDX}, startDatenums(mytimeIDX), endDatenums(mytimeIDX), NaN);
-      %w{mytimeIDX} = w0;
-      clear w0
-   end %mytimeIDX
+    else
+        debug.print_debug(0, 'Time subset failed to find any matching records in wfdisc table')
+        w0=waveform(); 
+    end
+    w = combine(w0);
+    %w = pad(w, startDatenum, endDatenum, NaN);
+    %w = w0;
+    clear w0
+
    closeIfAppropriate(mydb, onFinishCloseDB);
 
    function cleanUpFail(varargin)
@@ -473,32 +478,64 @@ function cleanUpFail(varargin)
         fprintf('\n\nMatching trace objects are:\n');
     end
 
-    % create empty waveform variables
-    wt = repmat(waveform(),safe_dbnrecs(tr),1);
-
-    % load data and metadata from trace table into waveform objects
-    for cc=1:safe_dbnrecs(tr)
-        tr.record = cc-1;
-        [trnet, trsta, trchan, trtime, trendtime, trnsamp, trsamprate, trinstype, trcalib, trcalper, trresponse, trdatatype, trsegtype] = dbgetv(tr, 'net', 'sta', 'chan', 'time', 'endtime', 'nsamp', 'samprate','instype','calib','calper','response','datatype','segtype');
-        trtime = epoch2datenum(trtime);
-        trendtime = epoch2datenum(trendtime);
-        trdata=trextract_data( tr );
-        npts=length(trdata);
-        if debug.get_debug() > 0
-            fprintf('%s\t%s\t%s\t%s\t%d\t%f\t%d\n',trsta,trchan,datestr(trtime),datestr(trendtime),trnsamp,trsamprate,npts);       
+%     % create empty waveform variables
+%     wt = repmat(waveform(),safe_dbnrecs(tr),1);
+
+    % Load the trace table into a trace structure
+    st = tr2struct(tr);
+
+    for cc=1:numel(st)
+        y = st(cc).data;
+        calib = st(cc).calib;
+        if abs(calib) > 0
+            y = y * calib;
         end
-        if strcmp(trnet,'-')
-            trnet='';
-        end
-        [trunits, ~] = segtype2units(trsegtype);
-        wt(cc) = waveform(ChannelTag(trnet,trsta,'',trchan), trsamprate, trtime, trdata*trcalib, trunits);
-        wt(cc) = addfield(wt(cc),'calib', trcalib);
-        if trcalib~=0
+        [trunits, ~] = segtype2units(st(cc).segtype);
+        wt(cc) = waveform(ChannelTag(st(cc).net, st(cc).sta, '',st(cc).chan), st(cc).samprate, epoch2datenum(st(cc).time), y, trunits);
+
+        if calib~=0
+            wt(cc) = addfield(wt(cc),'calib', st(cc).calib);
             wt(cc) = addfield(wt(cc), 'calibration_applied', 'YES');
         else
             wt(cc) = addfield(wt(cc), 'calibration_applied', 'NO');
         end
     end
+
+%         [trunits, ~] = segtype2units(trsegtype);
+%         wt(cc) = waveform(ChannelTag(trnet,trsta,'',trchan), trsamprate, trtime, trdata*trcalib, trunits);
+%         wt(cc) = addfield(wt(cc),'calib', trcalib);
+%         if trcalib~=0
+%             wt(cc) = addfield(wt(cc), 'calibration_applied', 'YES');
+%         else
+%             wt(cc) = addfield(wt(cc), 'calibration_applied', 'NO');
+%         end    
+
+
+%     % load data and metadata from trace table into waveform objects
+%     for cc=1:safe_dbnrecs(tr)
+%         tr.record = cc-1;
+%         
+%         pause(5)
+%         [trnet, trsta, trchan, trtime, trendtime, trnsamp, trsamprate, trinstype, trcalib, trcalper, trresponse, trdatatype, trsegtype] = dbgetv(tr, 'net', 'sta', 'chan', 'time', 'endtime', 'nsamp', 'samprate','instype','calib','calper','response','datatype','segtype');
+%         trtime = epoch2datenum(trtime);
+%         trendtime = epoch2datenum(trendtime);
+%         trdata=trextract_data( tr );
+%         npts=length(trdata);
+%         if debug.get_debug() > 0
+%             fprintf('%s\t%s\t%s\t%s\t%d\t%f\t%d\n',trsta,trchan,datestr(trtime),datestr(trendtime),trnsamp,trsamprate,npts);       
+%         end
+%         if strcmp(trnet,'-')
+%             trnet='';
+%         end
+%         [trunits, ~] = segtype2units(trsegtype);
+%         wt(cc) = waveform(ChannelTag(trnet,trsta,'',trchan), trsamprate, trtime, trdata*trcalib, trunits);
+%         wt(cc) = addfield(wt(cc),'calib', trcalib);
+%         if trcalib~=0
+%             wt(cc) = addfield(wt(cc), 'calibration_applied', 'YES');
+%         else
+%             wt(cc) = addfield(wt(cc), 'calibration_applied', 'NO');
+%         end
+%     end
     w = combine(wt); % combine waveforms based on ChannelTag (I think)
     clear wt
 

diff --git a/core/@waveform/spectrogram.m b/core/@waveform/spectrogram.m
@@ -1,8 +1,11 @@
-function spectrogram( w )
+function spectrogram( w, s )
 %SPECTROGRAM Plot an IceWeb-style spectrogram
-%   spectrogram(w) Creates an IceWeb style spectrogram by wrapping the
-%   function iceweb.spectrogram_iceweb(). For greater control, call that
-%   function directly, or use spectralobject/specgram or
+%   spectrogram(w, s) Creates an IceWeb style spectrogram by wrapping the
+%   function iceweb.spectrogram_iceweb(). If s is omitted it defaults to:
+%       spectralobject(1024, 924, 10, [60 120]);
+%
+% For greater control, call that
+%   iceweb.spectrogram_iceweb() directly, or use spectralobject/specgram or
 %   spectralobject/specgram2 (not clear how these differ). Note that
 %   spectrogram_iceweb() is significantly faster.
 
@@ -11,7 +14,9 @@ function spectrogram( w )
 if numel(w)>1
     w = reshape(w, numel(w), 1);
 end
-s = spectralobject(1024, 924, 10, [60 120]);
+if ~exist('s','var')
+    s = spectralobject(1024, 924, 10, [60 120]);
+end
 iceweb.spectrogram_iceweb(s, w, 0.75);
 
 
diff --git a/libgismo/epoch2datenum.m b/libgismo/epoch2datenum.m
@@ -6,7 +6,6 @@
 % $Revision$
 
 function time = epoch2datenum(epoch)
-
     if admin.antelope_exists()
         time = datenum(epoch2str(epoch,'%m %d %Y %H %M %S.%s'),'mm dd yyyy HH MM SS.FFF');
     else % does not handle leap seconds