System identification in z-domain


System identification in z-domain is performed with the function ao/zDomainFit. It is based on a modeified version of the vector fitting algorithm that was chenged to fit in z-domain. Details on the agorithm can be found in Z-domain fit help page.

System identification in z-domain - Exercise 1

During this exercise we will:

  1. Generate white noise
  2. Filter white noise with a miir filter
  3. Extract transfer function from data
  4. Fit TF with zDomainFit
  5. Check results

Let's open a new editor window and generate some white noise. Run...

  a = ao(plist('tsfcn', 'randn(size(t))', 'fs', 1, 'nsecs', 10000,'yunits','m'));

This command generates a time series of gaussian distributed random noise with a sampling frequency ('fs') of 1 Hz, 10000 seconds long ('nsecs') and assuming 'yunits' being meters ('m').

Now we are ready to move on the second step where we will:

Run the commands...

  plmiir = plist('a',[2 0.3 0.4],'b',[1 .05 .3 0.1],'fs',1,'iunits','m','ounits','V');

  filt = miir(plmiir);

  ac = filter(a,filt);
  ac.simplifyYunits;

We could calcualte data psd in order to check the difference with starting white noise. Run the commands...

  plpsd = plist('Nfft',1000);
  axx = psd(a,plpsd);

  acxx = psd(ac,plpsd);

  iplot(axx,acxx)

You should see something of similar to this picture.

Let's move to the third step. We will generate the transfer function from data and split them in order to remove the first 3 bins. The last operation is useful for the fitting process. Run the commands...

  pltf = plist('Nfft',1000);

  tf = tfe(a,ac,pltf);

  tf12 = split(tf(1,2),plist('split_type', 'frequencies', 'frequencies', [3e-3 5e-1]));

  iplot(tf(1,2),tf12)

The plot should report the as calculated and splitted transfer functions.

It is now the moment to start fitting with zDomainFit. As reported in the introduction above we can choose between two different criteria to check fit accuracy and to exit fitting loop.
We can start checking Log residuals difference and root mean squared error. It is worth to remember that:

Now let's run the commands...

  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');

When 'plot' parameter is set to on the function plot the fit progress.

It is now intersersting to run the fit and checking for Residuals spectral flatness and root mean squared error. It is worth to remember that:

Start a new fitting session...

  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');

It is now the moment to check the result of our fitting procedures. We will calculate the responses of the fitted models and compare them with the starting miir filter.
It is worth to note that fobj1 and fobj2 are parallel bank of miir filters.

  % 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))

The first plot show the respons of the three miir objects whereas the second plot reports the difference between fitted models and target filter. The result give an idea of the accuracy of the identification procedure.

The two fitting procedures came together to the same model.

In order to have more insight of the fit results, we look at the filters coefficients. The 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.
The algorithm works fine with low order models. For higher order model the algorithm is capable to provide a model within the prescribed accuracy but the order could be not guaranteed. Let's see what happen with higher order models creating a new editor window and running the second exercise.

System identification in z-domain - Exercise 2

During this exercise we will:

  1. Load fsdata object from file
  2. Fit loaded TF data with zDomainFit and fixed order
  3. Check results
Start loading AO...

  rfilt = ao(plist('filename', 'topic5\T5_Ex02_rfilt.xml'));
  iplot(rfilt)

The command loads a fsdata analysis object containing the frequency response of a 19 order iir filter expanded in partial fractions.
As we know the order of the model we want to fit we can set set 'Autosearch' to off status. In this case the function does not perform the accuracy test but simply run how far maximum number of iteration is reached. Model order is fixed by 'minorder' parameter.
Start fit with the commands...

  plfit1 = plist('FS',10,... % Sampling frequency for the model filters
    'AutoSearch','off',... % Automatically search for a good model
    'StartPolesOpt','c1',... % Define the properties of the starting poles - complex distributed in the unitary circle
    'maxiter',30,... % maximum number of iteration per model order
    'minorder',19,... % fixed model order
    'weightparam','abs',... % assign weights as 1./abs(data)
    '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
  fobj = zDomainFit(rfilt,plfit1);

Fit result should look like

Once the fit is done, we want to check results.
We need to extract residues and poles from the fitted model and compare them with the expected results. Extract residues and poles running...

  fRes = zeros(numel(fobj),1); % fit residue vector initialization
  fPoles = zeros(numel(fobj),1); % fit poles vector initialization

  % extracting data from fitted filters
  for ii = 1:numel(fobj)
    fRes(ii,1) = fobj(ii).a(1);
    fPoles(ii,1) = -1*fobj(ii).b(2);
  end
  [fRes,idx] = sort(fRes);
  fPoles = fPoles(idx);

Nominal residues and poles are well known

  mRes = [2.44554138162509e-011 - 1.79482547894083e-011i;
  2.44554138162509e-011 + 1.79482547894083e-011i;
  2.66402334803101e-009 +  1.1025122049153e-009i;
  2.66402334803101e-009 -  1.1025122049153e-009i;
  -7.3560293387644e-009;
  -1.82811618589835e-009 - 1.21803627800855e-009i;
  -1.82811618589835e-009 + 1.21803627800855e-009i;
  1.16258677367555e-009;
  1.65216557639319e-016;                         
  -1.78092396888606e-016;
  -2.80420398962379e-017;
  9.21305973049041e-013 - 8.24686706827269e-014i;
  9.21305973049041e-013 + 8.24686706827269e-014i;
  5.10730060739905e-010 - 3.76571756625722e-011i;
  5.10730060739905e-010 + 3.76571756625722e-011i;
  3.45893698149735e-009;
  3.98139182134446e-014 - 8.25503935419059e-014i;
  3.98139182134446e-014 + 8.25503935419059e-014i;
  -1.40595719147164e-011];
  [mRes,idx] = sort(mRes);

  mPoles = [0.843464045655194 - 0.0959986292915475i;
  0.843464045655194 + 0.0959986292915475i;
  0.953187595424927 - 0.0190043625473383i;
  0.953187595424927 + 0.0190043625473383i;
  0.967176277937188;
  0.995012027005247 - 0.00268322602801729i;
  0.995012027005247 + 0.00268322602801729i;
  0.996564761885673;
  0.999999366165445;
  0.999981722418555;
  0.999921882627659;
  0.999624431675213 - 0.000813407848742761i;
  0.999624431675213 + 0.000813407848742761i;
  0.997312006278751 - 0.00265611346834941i;
  0.997312006278751 + 0.00265611346834941i;
  0.990516544257531;
  0.477796923118318 - 0.311064085401834i;
  0.477796923118318 + 0.311064085401834i;
  0];
  mPoles = mPoles(idx);

Final comparison can be performed with a check on the relative difference between fitted and nominal residues and poles

  (mRes-fRes)./abs(mRes)
  (mPoles-fPoles)./abs(mPoles)

On the command window should not appear numbers larger than 1e-7. This indicating the tesults are accurate to the 7th decimal digit.




©LTP Team