line source
% 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)