Generation of noise with given psd


Generation of model noise is performed with the function ao/noisegen1D. Details on the agorithm can be found in noisegen1D help page.

Generation of noise with given psd - Exercise

During this exercise we will:

  1. Load fsdata object from file
  2. Fit psd of test data with zDomainFit
  3. Genarate noise from fitted model
  4. Compare results

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

We can now start the noise generation process. The first step is to generate a white time series AO

  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




©LTP Team