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.
During this exercise we will:
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:
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. |
% 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%.