% matlab code originally written by Xiaofeng Meng
% modified by zpeng@gatech.edu

sacfile = 'BK.PKD.HHT.SAC.vel';
[t,data,SAChdr] = fget_sac(sacfile);

%sampling rate
Fs = 80; 

%remove the mean
M = mean(data);
data_rm = data-M;

%cut data with boxcar
w1 = (t<300 | t>=1800).*0+(t>=300 & t<1800).*1;

%apply cosine taper
w2 = (t<300 | t>=1800).*0+...
     (t>=300 & t<375).*(0.5+0.5*cos(pi*(t+675)/75))+...
     (t>=375 & t<1725).*1+...
     (t>=1725 & t<1800).*(0.5+0.5*cos(pi*(t-675)/75));
data_taper = data_rm.*w2;

%extend time series to the nearest power of 2
NFFT = 2^nextpow2(1500*Fs);                                                

%cut the data
data1 = data_taper(300*Fs:(300*Fs+NFFT)-1);

%fast fourier transform
L=NFFT;
y1 = fft(data1,L)/NFFT;
f = Fs/2*linspace(0,1,L/2+1);

figure
%amplitude spectrum
subplot 211
plot(f,(abs(y1(1:L/2+1))),'k');
set(gca,'xscale','log');
ylim([-0.001 0.015]);
title('Amplitude Spectrum of y(t)');
ylabel('log10(|Y(f)|)');
set(gca,'yscale','log');
%phase spectrum
subplot 212
plot(f,atan(imag(y1(1:L/2+1))./real(y1(1:L/2+1))),'k');
title('Phase Spectrum of y(t)');
xlabel('Frequency (Hz)');
ylabel('Phase');

%compare original data and inverse fft data
%inverse fourier transform
figure
xx = 0.69;
yy = 0.88;
y2=ifft(y1,L);
subplot 211
plot(t(300*Fs:(300*Fs+NFFT)-1),y2,'k');
text(xx,yy-0.1,'inverse data','FontSize',18,'units','normalized');
subplot 212
plot(t(300*Fs:(300*Fs+NFFT)-1),data1(:),'k');
xlabel('Time (s)');
text(xx,yy-0.1,'original data','FontSize',18,'units','normalized');

%band pass filter of 2-16 Hz
f1 = Fs/2*linspace(0,1,L);
bp=(f1<2 | f1>=16).*0+...
   (f1>=2 & f1<2.7).*(0.5+0.5*cos(pi*(f1+6.3)/0.7))+...
   (f1>=2.7 & f1<15.3).*1+...
   (f1>=15.3 & f1<16).*(0.5+0.5*cos(pi*(f1-6.3)/0.7));

%low pass filter of 1/30 Hz
lp=(f1>=0.0333).*0+...
   (f1>=0 & f1<0.00333/2).*(0.5+0.5*cos(pi*(2*f1+0.03)/0.00333))+...
   (f1>=0.00333/2 & f1<(0.0333-0.00333/2)).*1+...
   (f1>=(0.0333-0.00333/2) & f1<0.0333).*(0.5+0.5*cos(pi*(2*f1-0.03)/0.00333));

%apply band pass filter
f_bp=y1.*bp';
y_bp=ifft(f_bp);

%apply low pass filter
f_lp=y1.*lp';
y_lp=ifft(f_lp);

figure
subplot 311
plot(t(300*Fs:(300*Fs+NFFT)-1),y_bp,'k');
text(xx,yy-0.15,'bp 2-16 Hz','FontSize',18,'units','normalized');
subplot 312
plot(t(300*Fs:(300*Fs+NFFT)-1),y_lp,'k');
text(xx,yy-0.1,'lp 30s','FontSize',18,'units','normalized');
subplot 313
plot(t(300*Fs:(300*Fs+NFFT)-1),data1(:),'k');
text(xx,yy-0.1,'original data','FontSize',18,'units','normalized');

%plot the spectrum
figure
noverlap=128;
nfft=256;
L1=nfft;
f_spec = Fs/2*linspace(0,1,L1)';
y_spec=ones(L1,1280);
[B,F,T]=specgram(data,nfft,Fs,nfft,noverlap);
Bmax = max(max(B));

subplot 311
plot(t(300*Fs:(300*Fs+NFFT)-1),y_bp,'k');
text(xx,yy-0.15,'bp 2-16 Hz','FontSize',18,'units','normalized');

subplot 312
imagesc(T,F,log10(abs(B./Bmax))*10);
xlim([200 2000]);
axis xy,
ylim([0 20]);
caxis([-120 0]);
ylabel('Frequency (Hz)');

subplot 313
for i=1:1:1279
    y_spec(:,i) = fft(data(i*128:i*128+nfft-1),L1);
end
t_spec=0*1.6:1.6:1278*1.6;
ymax=max(max(y_spec));
imagesc(t_spec,f_spec,log10(abs(y_spec./ymax))*10);
ylim([0 20]);
xlim([200 2000]);
caxis([-80 0])
axis xy,
xlabel('Time (s)');
ylabel('Frequency (Hz)');
