%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SEISMIC TRACE ANALYSIS
% Basic seismic data QC
% Use plotseismic.m to plot data, i.e. have plotseismic.m script in the same directory
% Display power spectra
% -------------------------------------------------------------------------
clear;
%------------------------------------------------------------------------
% Read in seismic line and plot in time
%------------------------------------------------------------------------
load data
dt = 0.004 ; % seismic data sample interval in ms
n_traces = size(data,2); % number of traces in data
ns = size(data,1); % number of samples in data
p = nextpow2(ns);
nt =2^p; % determine nt power or 2 >= to ns
fnyq = 1.0/(2.0*dt); % determine Nyquist frequency
df = 1.0/(nt*dt); % determine fundamental frequency
t = [(0:nt-1)*dt]; % Create time series
t=t';
%fft_trace = fft(data,nt); % compute the fft of matrix data
fft_trace = fft(data(:,2),nt); % compute the fft of trace #2
figure(1)
gain = 1.5;
plotseismic(data,gain,dt), xlabel('trace number'), ylabel('two-way time (s)')
grid
power_sp = (abs(fft_trace)).*(abs(fft_trace)); % compute power spectrum (frequency spectrum ^2)
max_power = (max(max(power_sp))); % maximum amplitude of power spectrum
figure(2) % Plot the normalized power spectrum
plot(0:df:2*fnyq-df , power_sp/max_power)
title 'Normalized Power Spectrum'
axis([0 fnyq 0 1])
grid