% 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