function gen_spectrogram(stn,chan,comp,net,tt1,tt2,ff1,ff2,caxis1,caxis2,color_flag)
% function gen_spectrogram(stn,chan,comp,net,tt1,tt2,ff1,ff2,caxis1,caxis2,color_flag)
% plot the raw, bandpass-filtered, and spectrogram of the seismogram
% this version could generate a ASCII table of the spectrogram
% that will be used in producing a GMT grd file
% db is from http://en.wikipedia.org/wiki/Decibel#Amplitude.2C_voltage_and_current
% applied a 0.5 Hz filter before computing spectrogram.
% see more at Peng et al. (SRL, 2012).
%  Peng, Z., L. T. Long, and P. Zhao* (2011), The relevance of high-frequency analysis artifacts to remote triggering, Seismol. Res. Lett., 82(5), 856-662, doi: 10.1785/gssrl.82.5.656.
% http://geophysics.eas.gatech.edu/people/zpeng/zpeng_paper/Peng_etal_SRL_2011.pdf
% need: eqfiltfilt.m at http://geophysics.eas.gatech.edu/people/zpeng/Software/matlab/eqfiltfilt.m
% this program can be downloaded at
% http://geophysics.eas.gatech.edu/people/zpeng/Software/matlab/gen_spectrogram.m
% update by zpeng on 08/07/2013

if nargin < 11, color_flag = 1;  end
if nargin < 10, caxis2 = -10;  end
if nargin < 9, caxis1 = -60; end
if nargin < 8, ff2 = 30;  end
if nargin < 7, ff1 = 0; end
if nargin < 6, tt2 = 2000;  end
if nargin < 5, tt1 = 0; end
if nargin < 4, net = 'BK'; end
if nargin < 3, comp = 'T'; end
if nargin < 2, chan = 'HH'; end
if nargin < 1, stn = 'PKD'; end

% add a new flag to apply a 0.5 Hz filtering to the data first
% before computing spectrogram
hp_flag = 1;

%
% do the test
ifsave = 0; % save the spectrogram output or not, default is not
event = 'Denali_20021103';
event1 = 'Denali\_20021103';
prefix = './';
postfix = 'SAC';

sacfile = sprintf('%s%s.%s.%s%s.%s',prefix,net,stn,chan,comp,postfix);
[time,data,SAChdr] = fget_sac(sacfile);
data = data - mean(data);
if (hp_flag), 
ch1 = 0.5;
ch2 = -inf;
datahp=eqfiltfilt(data,ch1,ch2,SAChdr.times.delta,4);
timehp = time;
end;

% compute the bp data (2-16 Hz band-pass filter)
ch1 = 2; ch2 = 16;
databp = eqfiltfilt(data,ch1,ch2,SAChdr.times.delta,4);
timebp = time;
% or replace it with the envelope function
% sacfilebp = sprintf('%s%s.%s.%s%s.%s.env',prefix,net,stn,chan,comp,postfix);
% [timebp,databp,SAChdrbp] = fget_sac(sacfilebp);


if (color_flag),
out_fig_all = sprintf('%s_%s_%s%s_spec_color.eps',event,stn,chan,comp);
else,
out_fig_all = sprintf('%s_%s_%s%s_spec_bw.eps',event,stn,chan,comp);
end;

% [timebp,databp,SAChdrbp] = fget_sac(sacfilebp);

if (tt1<SAChdr.times.b), tt1 = SAChdr.times.b; end
if (tt2>SAChdr.times.e), tt2 = SAChdr.times.e; end

sample_rate = round(1/SAChdr.times.delta);
if (ff2>sample_rate/2), ff2 = sample_rate/2; end

nfft = 256;
window = nfft;
noverlap = nfft*0.75;

if (hp_flag),
% [time1,data1,SAChdr1] = fget_sac(sacfilehp1);
[B,F,T] = specgram(datahp,nfft,sample_rate,window,noverlap);
disp('enter here');
else,
[B,F,T] = specgram(data,nfft,sample_rate,window,noverlap);
end;

% [B,F,T] = spectrogram(data,256,[],256);
Bmax = max(max(B));

% set(gcf,'PaperOrientation','landscape');
figure(1);
% set(gcf,'PaperOrientation','portrait','FontName','Times');
set(gcf,'PaperOrientation','portrait');
newpp = [0.3 1.0 7.5 9.5];
set(gcf,'paperPosition',newpp);

xx = 0.79;
yy = 0.88;
xx1 = 0.12;
yy1 = 1.10;
xx2 = 0.01;
yy2 = 0.92;

xx3 = tt1-(tt2-tt1)*0.05;
max_data = max(data(time>=tt1 & time<=tt2));
subplot 311;
%plot(timeT,dataT-max_data,'k');
%hold on;
plot(time,data,'k');
%plot(timeZ,dataZ+max_data,'k');
xlim([tt1 tt2]);
h = ylabel('Amplitude (counts)','FontSize',17);
set(h,'Position',[xx3   max_data*(-1.0)    1.0001]);
ylim([-1.1 1.1].*max_data);
% title(['Event ',event1,', Station ',stn],'FontSize',16);
text(xx1,yy1,['Event ',event1,', Station ',net,'.',stn],'FontSize',18,'units','normalized');
% text(xx1,yy1,'2002 Mw7.8 Denali earthquake triggered tremors at Parkfield, CA','FontSize',17,'units','normalized');
%text(xx,0.75,[stn,' ',chan,compZ],'FontSize',18,'units','normalized');
text(xx,0.65,[stn,' ',chan,comp],'FontSize',18,'units','normalized');
%text(xx,0.35,[stn,' ',chan,compT],'FontSize',18,'units','normalized');
% text(xx,yy-0.15,'raw data','FontSize',16,'units','normalized');
text(xx2,yy2,'(a)','FontSize',18,'units','normalized');
% text(0.15,0.52,'\downarrow','FontSize',12,'units','normalized');
% text(0.08,0.60,'\it{P} \rm{wave}','FontSize',18,'units','normalized');
% text(0.33,0.42,'\it{S} \rm{wave}','FontSize',18,'units','normalized');
% text(0.40,0.37,'\downarrow','FontSize',12,'units','normalized');
% text(0.32,0.07,'Love wave','FontSize',18,'units','normalized');
% text(0.45,0.95,'Surface wave','FontSize',18,'units','normalized');
set(gca,'Position',[0.1 0.65 0.80 0.30],'XTickLabel',[],'FontName','Times');
% text(0.56,0.10,'Rayleigh','FontSize',16,'units','normalized');

subplot 312;
plot(timebp,databp,'k');
max_databp = max(databp(timebp>=tt1 & timebp<=tt2));
min_databp = min(databp(timebp>=tt1 & timebp<=tt2));
ylim([min_databp 1.05*max_databp]);
xlim([tt1 tt2]);
% ylabel('Amplitude (cm/s)','FontSize',16);
text(xx,yy,[stn,' ',chan,comp],'FontSize',18,'units','normalized');
% text(xx,yy-0.15,'hp 5 Hz','FontSize',18,'units','normalized');
text(xx,yy-0.15,'bp 2-16 Hz','FontSize',18,'units','normalized');
text(xx2,yy2,'(b)','FontSize',18,'units','normalized');
% text(0.05,0.25,'\it{P} \rm{wave}','FontSize',16,'units','normalized');
% text(0.30,0.10,'Triggered activity','FontSize',18,'units','normalized');
set(gca,'Position',[0.1 0.50 0.80 0.15],'XTickLabel',[],'FontName','Times')

subplot 313;
%specgram(data);
if (hp_flag),
imagesc(timehp,F,log10(abs(B./Bmax)).*(10));
else,
imagesc(time,F,log10(abs(B./Bmax)).*(10));
end;
if (color_flag), 
colormap(jet);
else,
colormap(gray);
cm=colormap(gray);
cm=cm(end:-1:1,:);
colormap(cm)
end;
% caxis([-60 -10])
caxis([caxis1 caxis2])
axis xy,
hc = colorbar('Location','EastOutside');
% set(hc,'YTick',[-30 -20 -10])

% h = colorbar('Location','East');

% hold on;
% plot(time,data./max_data*(ff2-ff1)/2+10,'r','LineWidth',1);
% h = liney(10);
% set(h,'Color',[1 1 1]);
% hold off;

xlabel('Time (s)','FontSize',17);
h = ylabel('Frequency (Hz)','FontSize',17);
yy3 = (ff1+ff2)/2;
set(h,'Position',[xx3 yy3 1.00005]);

xx = 0.77;
xlim([tt1 tt2]);
ylim([ff1 ff2]);
text(xx2,yy2+0.02,'(c)','FontSize',18,'units','normalized','BackgroundColor',[1 1 1]);
text(xx,yy2+0.02,[stn,' ',chan,comp,'     '],'FontSize',18,'units','normalized','BackgroundColor',[1 1 1]);
text(xx,yy2-0.08,'spectrogram','FontSize',18,'units','normalized','BackgroundColor',[1 1 1]);
set(gca,'Position',[0.10 0.20 0.80 0.30],'FontName','Times')
% text(0.03,0.93,'2002 Mw7.8 Denali EQ triggered tremors at Parkfield, CA','FontSize',16,'units','normalized','BackgroundColor',[1 1 1]);
%set(h,'Position',[0.88 0.1 0.08 0.5])

% text(0,-0.2,'Figure 2','FontSize',18,'units','normalized');

if (color_flag),
figure_id = 'Figure 1';
text(0,-0.3,figure_id,'FontSize',20,'units','normalized','BackgroundColor',[1 1 1]);
else,
figure_id = 'Figure 1';
text(0,-0.3,figure_id,'FontSize',20,'units','normalized','BackgroundColor',[1 1 1]);
end;

eval(['print -depsc2 ',out_fig_all]);

if ifsave,
out_file = sprintf('%s_%s_spectrogram.dat',event,stn);
fid = fopen(out_file,'w');
[N,M] = size(B);
for ii = 1:M,
   for jj = 1:N,
	fprintf(fid,'%9.4f %9.4f %9.4f\n',T(ii),F(jj),log10(abs(B(jj,ii)./Bmax).*10));
   end
end
fclose(fid);
end
