% irreg_demo
%
% demo of an irregularly sampled problem
%
% Author Piet M.T. Broersen, September 2008
%
% The three statements :
%           ****************************************************** 
%           *  [asel bsel sellog]=ARMAsel_mis(ti,xi,ARmax,Tr,w); *
%           *  psdsel=arma2psd(asel,bsel);                       *
%           *  corsel=arma2cor(asel,bsel,50);                    *
%           ****************************************************** 
%  compute the selected estimated spectrum and autocorrelation function.
%  Sellog gives additional information and plots spectra of all estimated
%  AR models
%
%  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(' ')

ARmax=3
w=1/2            %*Tres   slot width, fraction of Tres
Tr=1/2000

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

w=1/8
[air2 bir2 sellogir2] = ARMAsel_irreg(ti,xi,ARmax,Tr,w)

w=1/2            %*Tres   slot width, fraction of Tres
Tr2=1/4000

[air3 bir3 sellogir3] = ARMAsel_irreg(ti,xi,ARmax,Tr2,w)

[h_ir f_ir] = arma2psd(air,bir,1000,Tr);
h_ir2 = arma2psd(air2,bir2,1000,Tr);
[h_ir3 f_ir3] = arma2psd(sellogir3.ar_sel_corrected,1,1000,Tr2);

disp(' ')
disp('  rc AR_sel for w = 1/2 and for w = 1/8, Tr = 0.0005 s')
[dum rc_sel]=ar2arset(sellogir.ar);
[dum rc_sel2]=ar2arset(sellogir2.ar);
[dum rc_sel3]=ar2arset(sellogir3.ar_sel_corrected);
disp(rc_sel)
disp(rc_sel2)
disp('  rc AR_sel-corrected for w = 1/2, Tr = 0.00025 s')
disp(rc_sel3)

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