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