+
−% 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)
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−
+
−