% Test read and plot WILBER seismograms  Z component only, test spectrum calculation
%  SAGE 2000, L. Braile (modified 12/02/03)
% Mag 6.5 event off coast of central Chile
% Origin Time = 1998 sept 3 (day 246) 17:37:59.0
% -29.2934, -71.5471
% Station NNA (Nana, Peru) -11.9875, -76.8422
% Dist = 17.93 deg (1993 km), Az = 343 deg
% Start time = 17:40:54.0 (00:2:55 after event)
% dt = 0.05s, Nyquist = 10 HZ, npts = 10520
% Z component
load NNAz1.txt -ascii
x=NNAz1;
mx=mean(x);
x=x-mx;  % remove mean from data
% plot x (seismogram) with x axis number of samples
plot(x)
title('Original Seismogram')
xlabel('Number of Samples')
ylabel('Relative Amplitude')
pause
%
% resample; note this process risks aliasing effects
%
s = x(1202:2:9392);  %  s is now 4096 points long, has a sample interval
%  of 0.10 s and a Nyquist frequency of 5 Hz.
n = max(size(s));
N = 8192;
dt = 0.1;  % resampled data sample interval
fs = 1/dt; % sample rate
S=fft(s,N);  % set fft length to be larger than length of s 
fnyq=1/(2*dt);  % Nyquist frequency
freq = (0:N/2)/(N/2)*fnyq; %  set up frequency axis values
Sxx = S.*conj(S)/N;  % Calculate Power Spectral Density (Periodogram method,
% no smoothing or windowing)
%  denominator is scale factor
loglog(freq,Sxx(1:((N/2)+1)),'r-','linewidth',2)  % Plot spectrum from f = 0 to f = Nyquist
axis([0.01 10 10^0 10^12]);
hold on
% extract signals corresponding to noise, P arrivals and Rayleigh
% waves and compute and plot spectra
%
s1 = x(2:2:1400); % noise
n1 = max(size(s1));
S=fft(s1,N);  % set fft length to be larger than length of x 
fnyq=1/(2*dt);  % Nyquist frequency
freq = (0:N/2)/(N/2)*fnyq; %  set up frequency axis values
Sxx = S.*conj(S)/N;  % Calculate Power Spectral Density (Periodogram method,
% no smoothing or windowing)
%  denominator is scale factor
loglog(freq,Sxx(1:((N/2)+1)),'b-')  % Plot spectrum from f = 0 to f = Nyquist
%
s2 = x(1400:2:4000); % P arrivals
n2 = max(size(s2));
S=fft(s2,N);  % set fft length to be larger than length of x 
fnyq=1/(2*dt);  % Nyquist frequency
freq = (0:N/2)/(N/2)*fnyq; %  set up frequency axis values
Sxx = S.*conj(S)/N;  % Calculate Power Spectral Density (Periodogram method,
% no smoothing or windowing)
%  denominator is scale factor
loglog(freq,Sxx(1:((N/2)+1)),'g--')  % Plot spectrum from f = 0 to f = Nyquist
%
s3 = x(6200:2:10520); % Rayleigh wave
n3 = max(size(s3));
S=fft(s3,N);  % set fft length to be larger than length of x 
fnyq=1/(2*dt);  % Nyquist frequency
freq = (0:N/2)/(N/2)*fnyq; %  set up frequency axis values
Sxx = S.*conj(S)/N;  % Calculate Power Spectral Density (Periodogram method,
% no smoothing or windowing)
%  denominator is scale factor
loglog(freq,Sxx(1:((N/2)+1)),'k:')  % Plot spectrum from f = 0 to f = Nyquist
title('Power Spectral Density: Solid red = Seismogram; Solid blue = noise; Dashed = P Arrivals; Dotted = Rayleigh')
xlabel('Frequency (Hz)')
ylabel('Relative Amplitude')
hold off
%
pause
% Plot the time series
subplot(2,2,1), plot(s)
title('Resampled and Truncated Seismogram')
xlabel('Number of Samples')
ylabel('Relative Amplitude')
subplot(2,2,2), plot(s1)
title('Noise Sample')
xlabel('Number of Samples')
ylabel('Relative Amplitude')
subplot(2,2,3), plot(s2)
title('P Arrivals')
xlabel('Number of Samples')
ylabel('Relative Amplitude')
subplot(2,2,4), plot(s3)
title('Rayleigh Wave')
xlabel('Number of Samples')
ylabel('Relative Amplitude')