--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/LTPDA_training/topic5/TrainigSession_T5_Ex03.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,89 @@
+% 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)
+