--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/LTPDA_training/topic5/TrainigSession_T5_Ex01.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,192 @@
+% Training session Topic 5 exercise 01
+% 
+% System identification in z-domain 1
+% 
+% 1) Generate white noise
+% 2) Filter white noise with a miir filter
+% 3) Extract transfer function from data
+% 4) Fit TF with zDomainFit
+% 5) Compare results
+% 6) Comments
+% 
+% L FERRAIOLI 22-02-09
+%
+% $Id: TrainigSession_T5_Ex01.m,v 1.2 2009/02/25 18:18:45 luigi Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% 1) Generate White Noise
+
+% Build a gaussian random noise data series
+a = ao(plist('tsfcn', 'randn(size(t))', 'fs', 1, 'nsecs', 10000,'yunits','m'));
+
+%% 2) Filtering White Noise
+
+% set a proper plist for a miir filter converting meters to volt
+plmiir = plist('a',[2 0.3 0.4],'b',[1 .05 .3 0.1],'fs',1,'iunits','m','ounits','V');
+
+% Building the filter
+filt = miir(plmiir);
+
+% filtering white noise
+ac = filter(a,filt);
+ac.simplifyYunits;
+
+%% Checking
+% set a plist for psd calculation
+plpsd = plist('Nfft',1000);
+
+% calculate psd of white data
+axx = psd(a,plpsd);
+
+% calculate psd of colored data
+acxx = psd(ac,plpsd);
+
+% plotting psd
+iplot(axx,acxx)
+
+%% 3) Extracting transfer function
+
+% settung plist for TF estimation
+pltf = plist('Nfft',1000);
+
+% TF estimation
+tf = tfe(a,ac,pltf);
+
+% Removing first bins form tf
+tf12 = split(tf(1,2),plist('split_type', 'frequencies', 'frequencies', [3e-3 5e-1]));
+
+% plotting
+iplot(tf(1,2),tf12)
+
+%% 4.1) Fitting TF - checking residuals log difference
+
+% Perform an automatic fit of TF data in order to identify a discrete
+% model. The parameters controlling fit accuracy are
+% 
+% 'ResLogDiff' check that the normalized difference between data and fit
+% residuals in log scale is larger than the value specified.
+% if d is the input value for the parameter, than the algorithm check that
+% log10(data) - log10(res) > d
+% In other words if d = 1 an accuracy of 10% is guaranteed
+% if d = an accuracy of 1% is guaranteed
+% Note1 - fit residuals are calculated as res = model - data
+% Note2 - In fitting noisy data, this option should be relaxed otherwise
+% the function try to fit accurately the random oscillations due to the
+% noise
+% 
+% 'RMSE' if r is the parameter value, then the algorithm check that the
+% step-by-step Root Mean Squared Error variation is less than 10^-r 
+
+plfit1 = plist('FS',1,... % Sampling frequency for the model filters
+  'AutoSearch','on',... % Automatically search for a good model
+  'StartPolesOpt','c1',... % Define the properties of the starting poles - complex distributed in the unitary circle
+  'maxiter',50,... % maximum number of iteration per model order
+  'minorder',2,... % minimum model order
+  'maxorder',9,... % maximum model order
+  'weightparam','abs',... % assign weights as 1./abs(data)
+  'ResLogDiff',2,... % Residuals log difference
+  'ResFlat',[],... % Rsiduals spectral flatness (no need to assign this parameter for the moment)
+  'RMSE',7,... % Root Mean Squared Error Variation
+  'Plot','on',... % set the plot on or off
+  'ForceStability','on',... % force to output a stable ploes model
+  'CheckProgress','off'); % display fitting progress on the command window
+
+% Do the fit
+fobj1 = zDomainFit(tf12,plfit1);
+% setting input and output units for fitted model
+fobj1(:).setIunits('m');
+fobj1(:).setOunits('V');
+
+%% 4.2) Fitting TF - checking residuals spectral flatness
+
+% Perform an automatic fit of TF data in order to identify a discrete
+% model. The parameters controlling fit accuracy are
+% 
+% 'ResFlat' Check that normalized fit residuals spectral flatness is larger
+% than the parameter value. Admitted values for the parameter are 0<r<1.
+% Suggested values are 0.5<r<0.7.
+% Note - normalized fit residuals are calculated as
+% res = (model - data)./data
+% Note2 - This option could be useful with noisy data, when ResLogDiff
+% option is not useful
+% 
+% 'RMSE' if r is the parameter value, then the algorithm check that the
+% step-by-step Root Mean Squared Error variation is less than 10^-r 
+
+plfit2 = plist('FS',1,... % Sampling frequency for the model filters
+  'AutoSearch','on',... % Automatically search for a good model
+  'StartPolesOpt','c1',... % Define the properties of the starting poles - complex distributed in the unitary circle
+  'maxiter',50,... % maximum number of iteration per model order
+  'minorder',2,... % minimum model order
+  'maxorder',9,... % maximum model order
+  'weightparam','abs',... % assign weights as 1./abs(data)
+  'ResLogDiff',[],... % Residuals log difference (no need to assign this parameter for the moment)
+  'ResFlat',0.65,... % Rsiduals spectral flatness
+  'RMSE',7,... % Root Mean Squared Error Variation
+  'Plot','on',... % set the plot on or off
+  'ForceStability','on',... % force to output a stable ploes model
+  'CheckProgress','off'); % display fitting progress on the command window
+
+fobj2 = zDomainFit(tf12,plfit2);
+% setting input and output units for fitted model
+fobj2(:).setIunits('m');
+fobj2(:).setOunits('V');
+
+%% 5) Evaluating accuracy
+
+% fobj1 and fobj2 are parallel bank of miir filters. Therofore we can
+% calculate filter response and compare with our starting model
+
+% set plist for filter response
+plrsp = plist('bank','parallel','f1',1e-5,'f2',0.5,'nf',100,'scale','log');
+
+% calculate filters responses
+rfilt = resp(filt,plrsp); % strating model response
+rfobj1 = resp(fobj1,plrsp); % result of first fit
+rfobj2 = resp(fobj2,plrsp); % result of second fit
+
+% plotting filters
+iplot(rfilt,rfobj1,rfobj2)
+
+% plotting difference
+iplot(abs(rfobj1-rfilt),abs(rfobj2-rfilt))
+
+%% 6) Comments
+
+% Decomposition in partial fractions of the starting filter provides
+% residues =
+% 
+%           0.531297296250259 -      0.28636710753776i
+%           0.531297296250259 +      0.28636710753776i
+%           0.937405407499481                         
+% 
+% 
+% poles =
+% 
+%           0.112984252455744 +     0.591265379634664i
+%           0.112984252455744 -     0.591265379634664i
+%          -0.275968504911487
+% 
+% If we look inside the fitted objects fobj1 and fobj2 the following
+% results can be extracted
+% rsidues =
+%           0.531391341497363 - i*0.286395031447695
+%           0.531391341497363 + i*0.286395031447695
+%           0.937201773620951
+% 
+% poles = 
+%           0.112971704620295 + i*0.591204610859687
+%           0.112971704620295 - i*0.591204610859687
+%          -0.276032231799774
+% 
+% Model order is correct and residues and poles values are accurate to the
+% third decimal digit.
+% 
+% With low order model the algorithm works fine. For higher order model the
+% algorithm is capable to provide a model within the prescribed accuracy
+% but the order could be not guaranteed. Lets see what happen with higher
+% order models opening the second exercise script:
+% TrainingSession_T5_Ex02.m
+
+