| LTPDA Toolbox | contents | ![]() |
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.
During this exercise we will:
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:
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:
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:
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.
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.
During this exercise we will:
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
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.
|
Topic 5 - Model fitting | Generation of noise with given psd | ![]() |
©LTP Team