%    Tests the Lomb Periodogram algorithm and plots results
%
%	first, create a random set of times from 0 to 1 seconds with
%		the rand function, and sort in ascending order
%
t=rand(100,1); t=sort(t);
%			choose two frequencies to build data
f1=50; f2=130;
%			then create data, with some noise
y=sin(2*pi*f1*t)+2*sin(2*pi*f2*t)+0.5*randn(size(t));

%	next, plot the data

subplot(221); plot(t,y,'r*');		% plot data
set(gca,'LineWidth',2','FontSize',16);
xlabel('Time (sec)'); ylabel('Observation'); title('Input Data')

%	now plot the time interval distributions

subplot(222); hist(log10(diff(t)));
xlabel('Log10(Delta-T) (s)'); title('Intervals');
%
%		Now do compute the Lomb normalized periodogram
%
%	first create a vector of frequency bins, 1 per Hz
f=[0:200];
[Pn Prob]=lomb(t,y,f);

%	and plot the periodogram
%
subplot(223)
h=plot(f,Pn,'r'); set(h,'LineWidth',2);
xlabel('Frequency (Hz)'); ylabel('Normalized PSD'); 
title('Lomb Periodogram'); set(gca,'LineWidth',2,'FontSize',16); 
%
%	and plot the probabilities
%
subplot(224)
h=plot(f,Prob,'r'); set(h,'LineWidth',2);
xlabel('Frequency (Hz)'); ylabel('Probability');
set(gca,'LineWidth',2,'FontSize',16); 
%
% 	then sort for smallest values first 
%		(ie. those points least likely to be random)
%
[p,ind]=sort(Prob);
most=ind(1:5);				% just the first 5 values
ANS=[f(most)' Pn(most)' Prob(most)']