function [Px, Prob] = lomb(t, x, f)
% [Px Prob] = lomb(t, x, freq)
%
% Calculates the Lomb-Scargle normalized periodogram values
% "Px" as a function of the supplied vector of frequencies
% "freq" for input vectors "t" (time) and "x" (observations).
% Also returns the probability "Prob" that the null hypothesis
% is valid (same length as Px and freq). Time stamps, t and
% amplitudes "x" must be the same length.
%
% See Scargle J.D.:"Studies in astronomical time series analysis. II.
% Statistical aspects of spectral analysis of unevenly spaced data,"
% Astrophysical Journal, vol 263, pp. 835-853, 1982. ... and
% Lomb N.R: "Least-squares frequency analysis of unequally spaced data",
% Astrophysical and Spcae Science, vol 39, pp. 447-462, 1976.
% This file is an adaptation of the mysterious lomb.m which was emailed
% to me some time ago by a colleague who obtained it from a forgotten source.
% I claim no responsibility for its accuracy, although it seems to
% correspond with lomb.c from NRinC (but not fasper.c)
% Any information you have on this file, please email me:
% gari AT physionet DOT org
if nargin < 2
error('must have an amplitude for each time stamp')
end
% If no frequency vector is supplied, invent a default up to the
% highest frequency available (> Average Nyquist)
if nargin < 3
maxfreq = 1/min(diff(t));
f = [1/512:1/512:maxfreq];
end
% check length of inputs
if length(t) ~= length(x);
error('t and x not same length');
end;
% subtract mean, compute variance, initialize Px
z = x - mean(x);
var = std(x);
N = length(f);
Px = zeros(size(f));
% compute power by looping over all frequencies
for i=1:length(f)
w=2*pi*f(i);
if w > 0
twt = 2*w*t;
tau = atan2(sum(sin(twt)),sum(cos(twt)))/2/w;
wtmt = w*(t - tau);
Px(i) = (sum(z.*cos(wtmt)).^2)/sum(cos(wtmt).^2) + ...
(sum(z.*sin(wtmt)).^2)/sum(sin(wtmt).^2);
else
Px(i) = (sum(z.*t).^2)/sum(t.^2);
end
end
% normalize by variance and compute probabilities at each frequency
Prob = zeros(1,length(Px));
for i=1:length(Px)
if var~=0 % check for divide by zero
Px(i)=Px(i)/2/var.^2;
Prob(i) = 1-(1-exp(-Px(i)))^N;
else
Px(i)=inf;
Prob(i)=1;
end;
if Prob(i) < .001 % allow for possible roundoff error
Prob(i) = N*exp(-Px(i));
end
end;