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 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

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);

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




©LTP Team