% simple_irreg_demo
%
% most simple demo of spectral analysis for an irregularly sampled problem
%
% Author Piet M.T. Broersen, October 2008
%
% Two different solutions can be given to estimate the spectral density,
% depending on the highest frequency that is desired in the spectrum.
%
% -- The first is for low frequencies, lower than f0/6 Hz,
% where f0 is the mean data rate, the average irregular frequency.
% It uses equidistant Nearest Neighbor interpolation
% to obtain an equidistant signal that is analysed by ARMAsel
% which is available at http://www.mathworks.com/matlabcentral/
% The highest more or less reliable frequency is f0/6 Hz, but spectral
% details can be seen with considerable bias until f0/2 Hz.
%
% -- The second method is for frequencies higher than f0/2 Hz.
% ARMAsel-irreg uses multi-shift slotted Nearest Neighbor resampling MSSNNR
% that gives a missing data signal and requires a non-linear solution.
% Computation times are longer.
%
% The three statements :
% ***********************************************
% * [asel bsel]=ARMAsel_mis(ti,xi,ARmax,Tr,w); *
% * psdsel=arma2psd(asel,bsel,Lspec,Tr); *
% * corsel=arma2cor(asel,bsel,Lcor); *
% ***********************************************
% compute the selected estimated spectrum and autocorrelation function.
%
% ti : observation times [s]
% xi : observation amplitudes
% ARmax : highest 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
disp('Number of irregular 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('Mean data rate ( = 1/T0 )')
f0=1/T0
disp(' ')
Lspec=1000 % number of frequencies for spectral computation
Lcor=50 % number of lags for autocorrelation computation
% NN resampling with ARMAsel until frequency f0/6
% The estimated spectrum may have a bias due to aliasing of
% about a factor 2 at the highest frequency f0/6
disp(' ')
disp('*****')
disp('NN interpolation with frequency f0/3, for reliable but aliased spectrum')
disp('The estimated spectrum may be influenced by aliasing')
disp('Strong aliasing may cause a white noise spectrum')
disp(' ')
disp('The resampled signal has no missing values and can be analyzed')
disp('with the armasel program for contiguous equidistant data')
[xnn Tr]=nnresample(ti,xi,3);
[ann bnn]=armasel(xnn')
[psdnn fasnn]=arma2psd(ann,bnn,Lspec,Tr);
% NN resampling with ARMAsel until frequency f0/2
% The estimated spectrum will have a stronger resampling bias
% but spectral peaks may be visible in the range f0/6 - f0/2
% The aliasing bias is reduced
disp(' ')
disp('*****')
disp('NN interpolation with frequency f0, for spectrum with resampling bias')
disp('Strong resampling bias ia caused by multiple substitution of the same')
disp('irregular observation at different resampling instants.')
[xnn Tr]=nnresample(ti,xi,1);
[ann1 bnn1]=armasel(xnn')
[psdnn1 fasnn1]=arma2psd(ann1,bnn1,Lspec,Tr);
figure
semilogy(fasnn,psdnn,fasnn1,psdnn1)
title(['PSD of ',int2str(N),' irregular benchmark data, with NN resampling and armasel'])
xlabel('\rightarrow frequency [Hz]')
ylabel('\rightarrow Logarithm of power spectral density')
legend('strong alias, until \itf\rm_0/6','resampling bias, until \itf\rm_0/2',3)
axis tight
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(' ')
% ARMAsel-irreg can be used for higher frequencies
% This is a missing data non-linear estimation program
% A simple demo with low order estimates is given
% The highest frequency is equal to the mean data rate f0
% by using T0/2 as the resampling distance
%
disp('*****')
disp('MSSNNR resampling with ARMAsel-irreg and missing data')
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,Lspec,Tr);
cor_ir=arma2cor(air,bir,Lcor);
figure
loglog(f_ir,h_ir)
legend('\it{T_r = T}\rm_{0}/2 s, \itw\rm = 1/2 * \it{T_r}\rm s',2)
title(['PSD of ',int2str(N),' irregular benchmark data with ARMAsel-irreg.'])
xlabel('\rightarrow frequency [Hz]')
ylabel('\rightarrow Logarithm of power spectral density')
axis tight
figure
plot([0:50]*Tr,cor_ir)
title(['Autocorrelation function of ',int2str(N),' irregular benchmark data with ARMAsel-irreg.'])
xlabel('\rightarrow time lag [s]')
axis tight
