System identification in z-domain


System identification in Z-domain is performed with the function ao/zDomainFit. It is based on a modified version of the vector fitting algorithm that was changed to fit in the Z-domain. Details of 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 generated by a pzmodel
  3. Extract the transfer function from data
  4. Fit the transfer function with ao/zDomainFit
  5. Check results

Let's start by generating 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 with 'yunits' of metres ('m').

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

Run the commands...

  
    pzm = pzmodel(1, [0.005 2], [0.05 4]);

    filt   = miir(pzm, plist('fs', 1));
    filt.setIunits('m');
    filt.setOunits('V');

    % Filter the data
    ac = filter(a,filt);
    ac.simplifyYunits;

We can calcualte PSD of the data in order to check the difference between the coloured noise and the white noise.

    axx  = lpsd(a);
    acxx = lpsd(ac);
    iplot(axx,acxx)

You should see something similar to this picture.

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

    tf = ltfe(a,ac);
    tf = tf.index(1,2);
    tfsp = split(tf,plist('frequencies', [5e-4 5e-1]));
    
    iplot(tf,tfsp)

The plot look something like the following:

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 the fit accuracy and to exit the fitting loop.
We can start checking "Residual Spectral Flatness" and "Root Mean Squared Error". It is worth remembering that:

Key Description

ResLogDiff

This makes zDomainFit check that the normalized difference between the data and the fit residuals (in log scale) is larger than the value specified. If d is the input value for the parameter, then the algorithm checks that log10(data) - log10(res) > d.

RMSE

If r is the parameter value, then the algorithm checks that the step-by-step Root Mean Squared Error variation is less than 10^-r.

Now let's run the fit: we set 'ResLogDiff' to a small number (0.5) because with noisy data we do not want/expect a perfect match.

    % Set up the parameters
    plfit = 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',0.5,...      % Residuals log difference (must not assign this parameter for the moment)
      'ResFlat',[],...          % Residuals spectral flatness
      'RMSE',5,...              % Root Mean Squared Error Variation
      'Plot','on',...           % set the plot on or off
      'ForceStability','on',... % force to output a stable poles model
      'CheckProgress','off');   % display fitting progress on the command window
    
    % Do the fit
    fobj = zDomainFit(tfsp,plfit);
    
    % Set the input and output units for fitted model
    fobj.setIunits('m');
    fobj.setOunits('V');

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

We can now check the result of our fitting procedures. We will calculate the responses of the fitted models (filters) and compare them with the starting IIR filter, then we will plot the percentage error on the filters magnitudes.

Note that the result of the fitting procedure is a parallel bank of 3 IIR filters.

    % set plist for filter response
    plrsp = plist('bank','parallel','f1',1e-5,'f2',0.5,'nf',100,'scale','log');    
    
    % compute the response of the original noise-shape filter
    rfilt = resp(filt,plrsp);
    rfilt.setName;
    
    % compute the response of our fitted filter bank
    rfobj = resp(fobj,plrsp);
    rfobj.setName;
    
    % compare the responses    
    iplot(rfilt,rfobj)
    
    % and the percentage error on the magnitude        
    pdiff = 100.*abs((rfobj-rfilt)./rfilt);
    pdiff.simplifyYunits;
    iplot(pdiff,plist('YRanges',[1e-2 100]))

The first plot shows the response of the original filter and the fitted filter bank, whereas the second plot reports the percentage difference between fitted model and target filter magnitude.
As can be seen, the difference between filters magnitude is at most 10%.




©LTP Team