% simple matlab script to add 100 earthquakes together into a possible
% tremor signals
% homework 1, problem III

% read the template data
[t,data,hdr] = fget_sac('YJ.BV.EHN.SAC');

% normalize the data
data=data./max(data);
nn=length(t);

t_shift = 0.5;

N = 200; % 200 s data
n = round(N./hdr.times.delta)+1;
% generate a new time 
t0=0;
hd1 = newSacHeader(n,hdr.times.delta,t0);
newdata = zeros(n,1);
NN = round((N-t(end))./t_shift)+1;
for ii =1:NN,
	tnew=t0+(ii-1).*t_shift;
	nshift=round(tnew./hdr.times.delta);
%	disp([ii nshift tnew]);
	newdata(nshift+1:nshift+nn)=newdata(nshift+1:nshift+nn)+data;
%	disp([ii nshift tnew]);
end;
newsac=['tremor-',num2str(t_shift),'.sac'];
wtSac(newsac,hd1,newdata);

% convert it into wavefile

sac2wav(newsac,20,1,0,200);
time=[0:n-1].*hdr.times.delta;
% generate the spectrogram plot
hp_flag = 1;
if (hp_flag),
% sacfilehp1 = sprintf('%s%s.%s.%s%s.%s.hp1',prefix,net,stn,chan,comp,postfix);
ch1 = 0.5;
ch2 = -inf;
datahp=eqfiltfilt(newdata,ch1,ch2,hdr.times.delta,4);
timehp = time;
end;

nfft = 256;
window = nfft;
noverlap = nfft*0.75;
sample_rate = round(1/hdr.times.delta);

[B,F,T] = specgram(datahp,nfft,sample_rate,window,noverlap);
Bmax = max(max(B));

figure(1)
subplot 211;
plot(time,newdata);
ylabel('Amplitude','FontSize',16);
set(gca,'Position',[0.10 0.65 0.80 0.2],'FontName','Times')
title(['Stack time interval ',num2str(t_shift),' s'],'FontSize',18);
subplot 212;
imagesc(timehp,F,log10(abs(B./Bmax)).*(10));
colormap(jet);
caxis([-40 0])
axis xy,
hc = colorbar('Location','EastOutside');
ylabel('Frequency (Hz)','FontSize',16);
xlabel('Time (s)','FontSize',16);
set(gca,'Position',[0.10 0.10 0.80 0.50],'FontName','Times')

eval(['print -depsc tremorspec-',num2str(t_shift),'.eps']);
