Geoscience Reference
In-Depth Information
h e data contain periodicities of 100, 40 and 20 kyrs, as well as additive
Gaussian noise, and are unevenly spaced about the time axis. We dei ne two
new vectors
t
and
x
that contain the original time vector and the synthetic
oxygen-isotope data sampled at times
t
.
clear
series3 = load('series3.txt');
t = series3(:,1);
x = series3(:,2);
We then generate a frequency axis
f
. Since the Lomb-Scargle method is not
able to deal with the frequency of zero (i.e., with an ini nite period) we start
at a frequency value that is equivalent to the spacing of the frequency vector.
h evariable
ofac
is the oversampling parameter that inl uences the resolution
of the frequency axis about the
N(frequencies)=N(datapoints)
case. We also
need the highest frequency
fhi
that can be analyzed by the Lomb-Scargle
algorithm: the Nyquist frequency
fnyq
that would be obtained if the
N
data
points were evenly spaced over the same time interval is commonly used for
fhi
. h e following code uses the input parameter
hifac
, which is dei ned by
Press et al. (1992) as
hifac=fhi/fnyq
.
int = mean(diff(t));
ofac = 4; hifac = 1;
f = ((2*int)^(-1))/(length(x)*ofac): ...
((2*int)^(-1))/(length(x)*ofac): ...
hifac*(2*int)^(-1);
where
int
is the mean sampling interval. We normalize the data by subtracting
the mean.
x = x - mean(x);
We can now compute the normalized Lomb-Scargle periodogram
px
as a
function of the angular frequency
wrun
using the translation of Scargle's
FORTRAN code into MATLAB code.
for k = 1:length(f)
wrun = 2*pi*f(k);
px(k) = 1/(2*var(x)) * ...
((sum(x.*cos(wrun*t - ...
atan2(sum(sin(2*wrun*t)),sum(cos(2*wrun*t)))/2))).^2) ...
/(sum((cos(wrun*t - ...
atan2(sum(sin(2*wrun*t)),sum(cos(2*wrun*t)))/2)).^2)) + ...
((sum(x.*sin(wrun*t - ...
atan2(sum(sin(2*wrun*t)),sum(cos(2*wrun*t)))/2))).^2) ...
/(sum((sin(wrun*t - ...
atan2(sum(sin(2*wrun*t)),sum(cos(2*wrun*t)))/2)).^2));
end
