filename='2005.295.00.36.25.0000.YZ.BOOE..BHZ.D.SAC';
% new 24 data set
[sac_trace_structure] = msac_read(filename);
timevec = sac_trace_structure.time; 
% have determined that is 23.393 hours
trace = sac_trace_structure.trace;  % trace
dt = timevec(2)-timevec(1);
% sampling is 50Hz again
Fs=20000;
fac = max(trace)-min(trace)+1;  % 
%wavwrite(trace/fac,Fs,'new.wav');  % make a sound file

day = sac_trace_structure.nzjday;   % timestamp of beginning
hour = sac_trace_structure.nzhour;
min = sac_trace_structure.nzmin;
sec = sac_trace_structure.nzsec;
msec = sac_trace_structure.nzmsec;
jsec= dhmstos(day,hour,min,sec,msec); % 
%jsec is seconds from start of year
fs = 1.0/dt; % sampling rate Hz

trace = trace-median(trace); % remove DC offset
[sb,sa]=butter(4,1/(200.0*fs),'high'); %high pass at 200s
trace2 = filter(sb,sa,trace); % cut low,  filter out DC
[tb,ta]=butter(4,1.0/(5*fs),'low');  % low pass for power
sy=filter(tb,ta,abs(trace2));  % rectify and smooth
% sy looks real good.  identify all events that go above 0.1e5
% and don't let another event happen for 1000 samples?
% powerdb = 10*log10(sy);  % power as a function of time
%
% 

cuts=[0.25 2 5 10 18 25]; % filterbank cutoff frequencies
framehz=0.25;  % frequency in Hz we are taking frames in filterbank
[e,time]=auditoryfbank(trace2,fs,framehz,cuts);  % filterbanks



%[day,hour,min,sec,msec] = stodhms(ysec)
ns = length(time);  %  number of frame samples]]
% it is straightforward to only plot color/color for points 
% above an index ii = e(:,k) > median(e(:,k)) + 2*std(e(:,k));
% then plot(color14(ii),color23(ii),'d')