% irreg_demo
%
% most simple demo of an irregularly sampled problem
%
% Author Piet M.T. Broersen, September 2008
%
% The three statements :
%           *********************************************** 
%           *  [asel bsel]=ARMAsel_mis(ti,xi,ARmax,Tr,w); *
%           *  psdsel=arma2psd(asel,bsel);                *
%           *  corsel=arma2cor(asel,bsel,50);             *
%           *********************************************** 
%  compute the selected estimated spectrum and autocorrelation function.  
%
%  ti : observation times [s]
%  xi : observation amplitudes
%  ARmax : higherst candidate AR order
%  Tr : equidistant resampling distance [s]
%  w ; slot width as fraction of Tr [s]

clear all, close all, clc, echo off
   
load irreg_data
% row vectors of times ti and amplitudes xi of irregularly sampled observations

N=length(ti)

figure
plot(ti,xi,'p',ti,xi,':r')
title([int2str(N),' irregular benchmark  data.'])
xlabel('\rightarrow time axis [s]')
axis tight

disp('Mean time between observations')
T0=(ti(N) - ti(1))/(N-1)

disp(' ')
disp('*****')
disp('Some warning messages from the OPTIM toolbox cannot be suppressed easily')
disp('It is not necessary to provide gradient information')
disp('Warnings generated by the MATLAB OPTIM routine fminunc can be ignored')
disp('*****')
disp(' ')

% A simple demo with low order estimates.
% The highest frequency is equal to the mean data rate f0
% by using T0/2 as the resampling distance
%

ARmax=3
Tr=T0/2    % half of mean time between observations
w=1/2      % * Tr : slot width, as fraction of Tr

[air bir] = ARMAsel_irreg(ti,xi,ARmax,Tr,w)

[h_ir f_ir] = arma2psd(air,bir,1000,Tr);
cor_ir=arma2cor(air,bir,50);

figure
loglog(f_ir,h_ir)
legend('\it{Tr = T}\rm_{0}/2 s,  \itw\rm = 1/2 * \it{Tr}\rm s',2)
title(['PSD of ',int2str(N),' irregular benchmark data.'])
xlabel('\rightarrow frequency [Hz]')
ylabel('\rightarrow Logarithm of power spectral density')
axis tight

figure
plot(0:50,cor_ir)
title(['Autocorrelation function of ',int2str(N),' irregular benchmark data.'])
xlabel('\rightarrow time lag [s]')
axis tight