% 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 filetn = 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 correctlytnxx1 = tn.psd(plist('Nfft',2000));% Removing first binstnxx1r = split(tnxx1,plist('split_type', 'frequencies', 'frequencies', [2e-3 5e-1]));iplot(tnxx1r)% finding proper weights for noisy datastnxx = smoother(tnxx1r);wgh = 1./(abs(stnxx.data.y));%% 2.2) Fittingplfit1 = 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 noisea1 = 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 differenceac1 = noisegen1D(a1, plng);%% 4) Checking results% Calculating psd of generated dataacxx1 = ac1.psd(plist('Nfft',2000));% Comparing with starting dataiplot(tnxx1,acxx1)