Mercurial > hg > ltpda
view m-toolbox/test/LTPDA_training/TrainingSessionAll.m @ 44:409a22968d5e default
Add unit tests
author | Daniele Nicolodi <nicolodi@science.unitn.it> |
---|---|
date | Tue, 06 Dec 2011 18:42:11 +0100 |
parents | f0afece42f48 |
children |
line wrap: on
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)