--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/LTPDA_training/TrainingSessionAll.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,361 @@
+% 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)
+
+
+
+
+
+
+
+
+
+
+
+