+
−% Training session Topic 5 exercise 03
+
−% 
+
−% Generation of noise with given psd
+
−% 
+
−% 1) Load fsdata object from file
+
−% 2) Fit psd of test data with zDomainFit
+
−% 3) Genarate noise from fitted model
+
−% 4) Fit psd of test data with curvefit
+
−% 5) Genarate noise from fitted model
+
−% 6) Compare results
+
−% 
+
−% L FERRAIOLI 22-02-09
+
−%
+
−% $Id: TrainigSession_T5_Ex03.m,v 1.3 2009/02/25 18:18:45 luigi Exp $
+
−%
+
−%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
−
+
−%% 1) load fsdata
+
−
+
−% load test noise AO from file
+
−tn = ao(plist('filename', 'topic5\T5_Ex03_TestNoise.xml'));
+
−tn.setName;
+
−
+
−%% 2) Fitting psd of test data with zDomainFit
+
−
+
−%% 2.1) We try to identify a proper model for the psd of loaded data
+
−
+
−% making psd ot test data - we need some average otherwise the fitting
+
−% algorithm is not able to run correctly
+
−tnxx1 = tn.psd(plist('Nfft',2000));
+
−
+
−% Removing first bins
+
−tnxx1r = split(tnxx1,plist('split_type', 'frequencies', 'frequencies', [2e-3 5e-1]));
+
−iplot(tnxx1r)
+
−
+
−% finding proper weights for noisy data
+
−stnxx = smoother(tnxx1r);
+
−wgh = 1./(abs(stnxx.data.y));
+
−
+
−%% 2.2) Fitting
+
−
+
−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',10,... % minimum model order
+
−  'maxorder',45,... % maximum model order
+
−  'weights',wgh,... % assign externally calculated weights
+
−  'ResLogDiff',[],... % Residuals log difference (no need to assign this parameter for the moment)
+
−  'ResFlat',0.77,... % Rsiduals spectral flatness
+
−  'RMSE',5,... % Root Mean Squared Error Variation
+
−  'Plot','on',... % set the plot on or off
+
−  'ForceStability','off',... % we do not need a stable filter
+
−  'CheckProgress','off'); % display fitting progress on the command window
+
−
+
−% Do the fit
+
−[param1,fmod1] = zDomainFit(tnxx1r,plfit1);
+
−
+
−%% 3) Generating noise from model psd
+
−
+
−% noisegen1D is a noise coloring tool. It accept as input a white noise
+
−% time series and output a colored noise time series with one sided psd
+
−% matching the model specifiend into the parameters
+
−
+
−% Generate white noise
+
−a1 = ao(plist('tsfcn', 'randn(size(t))', 'fs', 1, 'nsecs', 10000,'yunits','m'));
+
−
+
−plng = plist(...
+
−    'model', abs(fmod1), ... % 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
+
−
+
−ac1 = noisegen1D(a1, plng);
+
−
+
−%% 4) Checking results
+
−
+
−% Calculating psd of generated data
+
−acxx1 = ac1.psd(plist('Nfft',2000));
+
−
+
−% Comparing with starting data
+
−iplot(tnxx1,acxx1)
+
−