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 10000s 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('frequencies', [2e-3 5e-1]));
iplot(tnxx,tnxxr)
The result 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 smooth our PSD data and then define the weights as the inverse of the absolute value of smoothed PSD. This should help the fit function to do a good job with noisy data.
stnxx = smoother(tnxxr); iplot(tnxxr, stnxx) wgh = 1./abs(stnxx);
The result of the smoother method is shown in the plot below:
Now let us run an automatic search for the proper model and pass the set of externally defined weights. The first output of zDomainFit is a miir filter model; 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);
Let us check the result. Calculate the PSD of the generated noise and compare it with PSD of starting data.
acxx = ac.psd(plist('Nfft',2000));
iplot(tnxx,acxx)
As can be seen, the result is in satisfactory agreement with the original data
![]() |
System identification in z-domain | Fitting time series with polynimials | ![]() |
©LTP Team