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