| LTPDA Toolbox | contents | ![]() |
Generation of model noise is performed with the function ao/noisegen1D. Details on the agorithm can be found in noisegen1D help page.
During this exercise we will:
Let's open a new editor window and load test data. Run...
tn = ao(plist('filename', 'topic5\T5_Ex03_TestNoise.xml'));
tn.setName;
This command will load an analysis object containing a test time series
10000 s long sampled at 1 Hz.
The command setName set the name of the AO to tn.
Now do the psd of our data. we need some average otherwise the fitting algorithm is not able to run correctly.
tnxx = tn.psd(plist('Nfft',2000));
In order to extract a reliable model from psd data we need to discharge the first bins.
tnxxr = split(tnxx,plist('split_type', 'frequencies', 'frequencies', [2e-3 5e-1]));
iplot(tnxx,tnxxr)
Resul should look like
Now it is the moment to fit our psd to extract a smooth model to pass to the noise generator. First of all we could define a set of proper weights for our fit process. We sooth our psd data and then define the weights as 1 over the absolute value of smoothed psd. This should help fit function to do a good job with noisy data.
stnxx = smoother(tnxxr); wgh = 1./(abs(stnxx.data.y));
Now let us run an automatic search for the proper model and pass the set of externally defined weights. First output of zDomainFit is a miir filter model where the second output is the model response.
plfit = plist('FS',1,...
'AutoSearch','on',...
'StartPolesOpt','c1',...
'maxiter',50,...
'minorder',10,...
'maxorder',45,...
'weights',wgh,... % assign externally calculated weights
'ResLogDiff',[],...
'ResFlat',0.77,...
'RMSE',5,...
'Plot','on',...
'ForceStability','off',...
'CheckProgress','off');
% Do the fit
[param,fmod] = zDomainFit(tnxxr,plfit);
Fit result should look like
a = ao(plist('tsfcn', 'randn(size(t))', 'fs', 1, 'nsecs', 10000,'yunits','m'));
Then start noise coloring process calling noisegen1D
plng = plist(...
'model', abs(fmod), ... % model for colored noise psd
'MaxIter', 50, ... % maximum number of fit iteration per model order
'PoleType', 2, ... % generates complex poles distributed in the unitary circle
'MinOrder', 20, ... % minimum model order
'MaxOrder', 50, ... % maximum model order
'Weights', 2, ... % weight with 1/abs(model)
'Plot', false,... % on to show the plot
'Disp', false,... % on to display fit progress on the command window
'RMSEVar', 7,... % Root Mean Squared Error Variation
'FitTolerance', 2); % Residuals log difference
ac = noisegen1D(a, plng);
Lut us check the result. Calculate psd of generated noise and compare it with psd of starting data.
acxx = ac.psd(plist('Nfft',2000));
iplot(tnxx,acxx)
As can be seen, result is in satisfactory agreement with original data
|
System identification in z-domain | Fitting time series with polynimials | ![]() |
©LTP Team