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