% Training session
%
% 1) Topic 1 - The basics of LTPDA
% 2) Topic 2 - Pre-processing of data
% 3) Topic 3 - Spectral Analysis
% 4) Topic 4 - Transfer function models and digital filtering
% 5) Topic 5 - Model fitting
%
% HISTORY:
%
% $Id: TrainingSessionAll.m,v 1.10 2010/02/24 17:55:39 ingo Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                      Topic 1 - The basics of LTPDA                      %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Clear all variables and figures
mc()
%%
% Reading the interferometer data
pl_file = plist('FILENAME', 'ifo_temp_example/ifo_training.dat', ...
                'TYPE',     'tsdata', ...
                'COLUMNS', [1 2],     ...
                'XUNITS', 's',        ...
                'YUNITS', '',         ...
                'ROBUST', 'no',       ...
                'DESCRIPTION', 'Interferometer data');
ifo = ao(pl_file);
ifo.setName(); % Set the object name to the variable name (here: 'ifo')

% Calibrating the interferometer data
lambda   = 1064e-9;
pl_scale = plist('factor', lambda/(2*pi), 'yunits', 'm');
ifo.scale(pl_scale);
ifo.setName(); % Set the object name to the variable name (here: 'ifo')

% Plot ifo
ifo.iplot(plist('XUNITS', 'h'));

% Save ifo to 'ifo_temp_example/ifo_disp.xml'
ifo.save('ifo_temp_example/ifo_disp.xml');


% Reading the interferometer data
pl_fileT = plist('FILENAME', 'ifo_temp_example/temp_training.dat', ...
                 'TYPE',     'tsdata', ...
                 'COLUMNS', [1 2],     ...
                 'XUNITS', 's',        ...
                 'YUNITS', 'degC',     ...
                 'ROBUST', 'no',       ...
                 'DESCRIPTION', 'Temperature data');
temp = ao(pl_fileT);

% Add offset
temp.offset(plist('offset', 273.15));
temp.setYunits('K');
temp.setName(); % Set the object name to the variable name (here: 'temp')

% Plot Tcel
temp.iplot(plist('XUNITS', 'h'));

% Save Tcel
temp.save(plist('filename', 'ifo_temp_example/temp_kelvin.xml'));


%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                     Topic 2 - Pre-processing of data                    %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Clear all variables and figures
mc()

% Load data from topic 1
ifo  = ao('ifo_temp_example/ifo_disp.xml');
temp = ao('ifo_temp_example/temp_kelvin.xml');

% plot the data
pl_plot1 = plist('arrangement', 'subplots');
iplot(ifo, temp, pl_plot1)

pl_plot2 = plist('ARRANGEMENT', 'subplots', ...
                 'LINESTYLES',  {'none','none'}, ...
                 'MARKERS',     {'+','+'}, ...
                 'XRANGES',     {'all', [200 210]}, ...
                 'YRANGES',     {[2e-7 3e-7], [200 350]});

iplot(ifo, temp, pl_plot2)

% The temperature data is unevenly sampled.
dt = diff(temp.x);
min(dt)
max(dt)

% Run 'data fixer' method ao/consolidate
[temp_fixed ifo_fixed] = consolidate(temp, ifo, plist('fs',1));

% Plot fixed data
iplot(ifo_fixed, temp_fixed, pl_plot1);
iplot(ifo_fixed, temp_fixed, pl_plot2);

% Save fixed data
save(temp_fixed,'ifo_temp_example/temp_fixed.xml');
save(ifo_fixed,'ifo_temp_example/ifo_fixed.xml');

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                       Topic 3 - Spectral Analysis                       %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Clear all variables and figures
mc()

% Get the consolidated data
% Using the xml format
T_filename = 'ifo_temp_example/temp_fixed.xml';
x_filename = 'ifo_temp_example/ifo_fixed.xml';

pl_load_T = plist('filename', T_filename);
pl_load_x = plist('filename', x_filename);

% Build the data aos
T = ao(pl_load_T);
x = ao(pl_load_x);

% PSD
x_psd = lpsd(x);
x_psd.setName('Interferometer');

T_psd = lpsd(T);
T_psd.setName('Temperature');

% Plot estimated PSD
pl_plot = plist('Arrangement', 'subplots', 'LineStyles', {'-','-'},'Linecolors', {'b', 'r'});
iplot(sqrt(x_psd), sqrt(T_psd), pl_plot);

% Skip some IFO glitch from the consolidation
pl_split = plist('split_type', 'interval', ...
                 'start_time', x.t0 + 40800, ...
                 'end_time', x.t0 + 193500);
x_red = split(x, pl_split);
T_red = split(T, pl_split);

% PSD
x_red_psd = lpsd(x_red);
x_red_psd.setName('Interferometer');

T_red_psd = lpsd(T_red);
T_red_psd.setName('Temperature');

% Plot estimated PSD
pl_plot = plist('Arrangement', 'stacked', 'LineStyles', {'-','-'},'Linecolors', {'b', 'r'});
iplot(sqrt(x_psd), sqrt(x_red_psd), pl_plot);
iplot(sqrt(T_psd), sqrt(T_red_psd), pl_plot);

% CPSD estimate
CTx = lcpsd(T_red, x_red);
CxT = lcpsd(x_red, T_red);

% Plot estimated CPSD
iplot(CTx);
iplot(CxT);

% Coherence estimate
coh = lcohere(T_red, x_red);

% Plot estimated cross-coherence
iplot(coh, plist('YScales', 'lin'))

% transfer function estimate
tf = ltfe(T_red, x_red);

% Plot estimated TF
iplot(tf);

% Noise projection in frequency domain
proj = T_red_psd.*(abs(tf)).^2;
proj.simplifyYunits;
proj.setName('temp. contrib. projection')

% Plotting the noise projection in frequency domain
iplot(x_red_psd, proj);

% Save the PSD data
% Plists for the xml format
pl_save_x_PSD     = plist('filename', 'ifo_temp_example/ifo_psd.xml');
pl_save_T_PSD     = plist('filename', 'ifo_temp_example/T_psd.xml');
pl_save_xT_CPSD   = plist('filename', 'ifo_temp_example/ifo_T_cpsd.xml');
pl_save_xT_cohere = plist('filename', 'ifo_temp_example/ifo_T_cohere.xml');
pl_save_xT_TFE    = plist('filename', 'ifo_temp_example/T_ifo_tf.xml');

x_red_psd.save(pl_save_x_PSD);
T_red_psd.save(pl_save_T_PSD);
CxT.save(pl_save_xT_CPSD);
coh.save(pl_save_xT_cohere);
tf.save(pl_save_xT_TFE);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%         Topic 4 - Transfer function models and digital filtering        %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mc()

% Temperature noise PZMODEL
TMP = pzmodel(10,1e-5,[]);
TMP.setName();
TMP.setOunits('K');

% Interferometer noise PZMODEL
IFO = pzmodel(1e-3, {0.4}, []);
IFO.setName();
IFO.setOunits('rad');

% Temperature to interferometer coupling PZMODEL
K2RAD = pzmodel(1e-1, {5e-4}, []);
K2RAD.setName();
K2RAD.setOunits('rad');
K2RAD.setIunits('K');

% Plot the response
pl = plist('f1',1e-5,'f2',0.01);
resp(K2RAD*TMP,IFO,pl);

% Discretize the three transfer (TMP,IFO,K2RAD) with the MIIR constructor
pl_miir = plist('fs', 1);
TMPd   = miir(TMP, pl_miir);
IFOd   = miir(IFO, pl_miir);
K2RADd = miir(K2RAD, pl_miir);

% Generate white noise with the AO constructor
pl_ao = plist('tsfcn', 'randn(size(t))', ...
              'fs',    1, ...
              'nsecs', 250000);
WN1 = ao(pl_ao);
WN2 = ao(pl_ao);

% Filter white noise WN1 with the TMP filter 
T = filter(WN1,TMPd);

% Filter white noise WN2 with the IFO filter
T2 = filter(WN2, IFOd);

% Filter white noise WN2 with the TMP and the K2RAD filter
T3 = filter(WN1, TMPd, K2RADd, plist('bank','serial'));

% Add Noise
IFO = T2 + T3;

% Split data stream
pl_split = plist('times', [1e5 2e5]);
IFO = IFO.split(pl_split);
IFO.setName('Interferometer');
T   = T.split(pl_split);
T.setName('Temperature');

% Plot
pl_plot1 = plist('ARRANGEMENT', 'subplots');
IFO.iplot(pl_plot1, T);

% Compute power spectral estimates for the temperature and interferometric data 
pl_lpsd = plist('order', 1, 'scale', 'ASD');
lpsd_T   = lpsd(T, pl_lpsd);
lpsd_ifo = lpsd(IFO, pl_lpsd);
iplot(lpsd_T, lpsd_ifo, plist('Arrangement', 'subplots'))

% Compute transfer function estimate for the temperature and interferometric data
tfe_T = ltfe(T, IFO);
tfe_T.setName('Transfer function');

pl = plist('f1',1e-5,'f2',1);
iplot(tfe_T, resp(K2RAD,pl));

% Compute projection
Projection = abs(tfe_T).*lpsd_T;
Projection.simplifyYunits;
Projection.setName;
% Plot against interferometer noise
iplot(lpsd_ifo,Projection);

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                         Topic 5 - Model fitting                         %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Clear all variables and figures
mc()

% Load test data from topic 1
ifo = ao(plist('filename', 'ifo_temp_example/ifo_fixed.xml'));
ifo.setName;
T = ao(plist('filename', 'ifo_temp_example/temp_fixed.xml'));
T.setName;

% Split out the good part of the data
pl_split = plist('split_type', 'interval', ...
                 'start_time', ifo.t0 + 40800, ...
                 'end_time',   ifo.t0 + 193500);

ifo_red = split(ifo, pl_split);
T_red   = split(T, pl_split);

% Plot
iplot(ifo_red,T_red,plist('arrangement', 'subplots'))

% Load transfer function from topic 3
tf = ao('ifo_temp_example/T_ifo_tf.xml');

% split the transfer function to extract only meaningful data
tfsp = split(tf,plist('frequencies', [2e-5 1e-3]));
iplot(tf,tfsp)

% force zDomainFit to fit a stable model
plfit = plist('FS',1,               ...
              'AutoSearch','off',   ...
              'StartPolesOpt','clin',...
              'maxiter',20,         ...
              'minorder',3,         ...
              'maxorder',3,         ...
              'weightparam','abs',  ...
              'Plot','on',          ...
              'ForceStability','on',...
              'CheckProgress','off');

fobj = zDomainFit(tfsp,plfit);

fobj.filters.setIunits('K');
fobj.filters.setOunits('m');

% Detrend after the filtering
ifoT = filter(T_red,fobj,plist('bank','parallel'));
ifoT.detrend(plist('order',0));
ifoT.simplifyYunits;
ifoT.setName;

% Substract temperature contribution from measured interferometer data
ifonT = ifo_red - ifoT;
ifonT.setName;

% Plot data
iplot(ifo_red,ifoT,ifonT)

% LPSD
ifoxx = ifo_red.lpsd;
ifonTxx = ifonT.lpsd;
iplot(ifoxx,ifonTxx)