--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_ao_zDomainFit_3.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,261 @@
+% Test script for zDomainFit on noisy data
+% 
+% L. Ferraioli 02-12-08
+% 
+% $Id: test_ao_zDomainFit_3.m,v 1.3 2009/04/20 14:31:14 luigi Exp $
+% 
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Smooth psd model
+% 
+% We define a smooth psd model and generate colored noise with such a psd
+% shape. Then we fit sqrt(psd) of noisy data following two different
+% procedure to decide when the fit is satisfactory and we can exit from the
+% fitting routine.
+% The first procedure is based on the check of residuals spectral flatness.
+% The second procedure is based on the check of residuals log distance to
+% fitted data.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Useful constants
+
+fs = 10; % sampling frequency Hz
+Nsecs = 1e5; % Number of seconds in the data stream
+Nfft = 1e4; % Number of samples in the fft
+pls  = plist('Nfft', Nfft,'Order',0); % plist for spectra
+
+%% Building a frequency response AO
+
+% Create a frequency-series AO
+% pl_data = plist('fsfcn', '0.01./(0.01+f)', 'f1', 1e-6, 'f2', 5, 'nf', 1000);
+pl_data = plist('fsfcn', '(1e-3./(f).^2 + 1e3./(0.001+f) + 1e5.*f.^2).*1e-10', 'f1', 1e-6, 'f2', 5, 'nf', 300);
+mod = ao(pl_data);
+mod.setName('psd model');
+
+iplot(mod)
+
+%% Building white noise
+
+a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
+a.setName('white noise')
+a.setYunits('m')
+
+%% Calling the noise generator
+
+pl = plist(...
+    'model', mod, ...
+    'MaxIter', 50, ...
+    'PoleType', 2, ...
+    'MinOrder', 15, ...
+    'MaxOrder', 25, ...
+    'Weights', 2, ...
+    'Plot', false,...
+    'Disp', false,...
+    'RMSEVar', 8,...
+    'FitTolerance', 2);
+
+ac = noisegen1D(a, pl);
+
+%% Spetrum and Transfer function
+
+acxx = ac.psd(pls); % spectrum of data
+etf = tfe(a,ac,pls); % calculating transfer functions from data
+
+iplot(acxx)
+iplot(etf(1,2))
+
+%% Fitting Sqrt Spectrum - Residuals spectral flatness checking
+
+% Fitting parameter list
+pl_fit1 = plist('FS',fs,...
+  'AutoSearch','on',...
+  'StartPoles',[],...
+  'StartPolesOpt','c2',...
+  'maxiter',60,...
+  'minorder',19,...
+  'maxorder',25,...
+  'weights',[],...
+  'weightparam','abs',...
+  'ResLogDiff',[],...
+  'ResFlat',0.7,...
+  'RMSE',5,...
+  'Plot','off',...
+  'ForceStability','off',...
+  'CheckProgress','on');
+
+% Do fit
+[filt2,bsp1,bspres1,bsprmse1] = zDomainFit(sqrt(acxx./2), pl_fit1); % division by 2 is necessary because acxx is the one-sided psd
+bsp1.setName('Fit1');
+
+%% Checing spectral flatness
+% checking the spectral flatness of the normalized residuals in order to be
+% sure that the procedure works correctly
+kk = abs(bspres1./bsp1);
+rf = utils.math.spflat(kk.data.y);
+iplot(abs(bspres1./bsp1))
+
+%% Fitting Spectrum - Residuals log difference checking
+
+% Fitting parameter list
+pl_fit2 = plist('FS',fs,...
+  'AutoSearch','on',...
+  'StartPoles',[],...
+  'StartPolesOpt','c2',...
+  'maxiter',60,...
+  'minorder',17,...
+  'maxorder',20,...
+  'weights',[],...
+  'weightparam','abs',...
+  'ResLogDiff',0.5,...
+  'ResFlat',0.7,...
+  'RMSE',5,...
+  'Plot','off',...
+  'ForceStability','off',...
+  'CheckProgress','on');
+
+% Do fit
+[filt2,bsp2,bspres2,bsprmse2] = zDomainFit(sqrt(acxx./2), pl_fit2); % division by 2 is necessary because acxx is the one-sided psd
+bsp2.setName('Fit2');
+
+%% Model calculated at the same frequencies of the fit AOs
+
+% This step is useful in order to compare the results of the two fitting
+% procedures
+f = acxx.data.x;
+pl_datar = plist('fsfcn', '(1e-3./(f).^2 + 1e3./(0.001+f) + 1e5.*f.^2).*1e-10', 'f', f);
+modr = ao(pl_datar);
+modr.setName('psd model');
+
+%% Checking the two fit results
+% plot the normalized difference between fit results and model data in
+% order to check the fit accuracy
+plpl = plist('Legends',{'Fit1','Fit2'});
+iplot(abs(bsp1-sqrt(modr))./sqrt(modr),abs(bsp2-sqrt(modr))./sqrt(modr),plpl)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Peaked psd model
+% 
+% We define a psd model with some resonance peaks and generate colored
+% noise with such a psd shape. Then we fit sqrt(psd) of noisy data
+% following two different procedure to decide when the fit is satisfactory
+% and we can exit from the fitting routine.
+% The first procedure is based on the check of residuals spectral flatness.
+% The second procedure is based on the check of residuals log distance to
+% fitted data.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% Building a frequency response AO
+
+% Create a psd model with two peak resonance at 5e-2 and 7e-2 Hz
+func = '(1e-3./(f).^2 + 1e3./(0.001+f) + 1e5.*f.^2).*1e-10 + 1e-6./((f./7e-2).^2-1).^2 + 1e-5./((f./5e-2).^2-1).^2';
+pl_data = plist('fsfcn', func, 'f1', 1e-6, 'f2', 5, 'nf', 300);
+mod = ao(pl_data);
+mod.setName('psd model');
+
+iplot(mod)
+
+%% Building white noise
+
+a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
+a.setName('white noise')
+
+%% Calling the noise generator
+
+pl = plist(...
+    'model', mod, ... % Multiplication by fs needed to preserve energy
+    'MaxIter', 70, ...
+    'PoleType', 2, ...
+    'MinOrder', 30, ...
+    'MaxOrder', 45, ...
+    'Weights', 2, ...
+    'Plot', false,...
+    'Disp', false,...
+    'RMSEVar', 8,...
+    'FitTolerance', 2);
+
+ac = noisegen1D(a, pl);
+
+%% Spetrum and Transfer function
+
+acxx = ac.psd(pls); % spectrum of data
+
+iplot(acxx)
+
+%% Fitting Sqrt Spectrum - Residuals spectral flatness checking
+
+% Fitting parameter list
+pl_fit1 = plist('FS',fs,...
+  'AutoSearch','on',...
+  'StartPoles',[],...
+  'StartPolesOpt','c2',...
+  'maxiter',60,...
+  'minorder',15,...
+  'maxorder',25,...
+  'weights',[],...
+  'weightparam','abs',...
+  'ResLogDiff',[],...
+  'ResFlat',0.7,...
+  'RMSE',6,...
+  'Plot','off',...
+  'ForceStability','off',...
+  'CheckProgress','on');
+
+% Do fit
+[filt1,bsp1,bspres1,bsprmse1] = zDomainFit(sqrt(acxx./2), pl_fit1); % division by 2 is necessary because acxx is the one-sided psd
+bsp1.setName('Fit1');
+
+%%
+
+iplot(sqrt(acxx./2),bsp1)
+
+%% Checing spectral flatness
+% checking the spectral flatness of the normalized residuals in order to be
+% sure that the procedure works correctly
+kk = abs(bspres1./bsp1);
+rf = utils.math.spflat(kk.data.y);
+iplot(abs(bspres1./bsp1))
+
+%% Fitting Spectrum - Residuals log difference checking
+
+% Fitting parameter list
+pl_fit2 = plist('FS',fs,...
+  'AutoSearch','on',...
+  'StartPoles',[],...
+  'StartPolesOpt','c2',...
+  'maxiter',60,...
+  'minorder',15,...
+  'maxorder',25,...
+  'weights',[],...
+  'weightparam','abs',...
+  'ResLogDiff',0.5,...
+  'ResFlat',0.7,...
+  'RMSE',6,...
+  'Plot','off',...
+  'ForceStability','off',...
+  'CheckProgress','on');
+
+% Do fit
+[filt2,bsp2,bspres2,bsprmse2] = zDomainFit(sqrt(acxx./2), pl_fit2); % division by 2 is necessary because acxx is the one-sided psd
+bsp2.setName('Fit2');
+
+%%
+
+iplot(sqrt(acxx./2),bsp2)
+
+%% Model calculated at the same frequencies of the fit AOs
+
+% This step is useful in order to compare the results of the two fitting
+% procedures
+f = acxx.data.x;
+func = '(1e-3./(f).^2 + 1e3./(0.001+f) + 1e5.*f.^2).*1e-10 + 1e-6./((f./7e-2).^2-1).^2 + 1e-5./((f./5e-2).^2-1).^2';
+pl_datar = plist('fsfcn', func, 'f', f);
+modr = ao(pl_datar);
+modr.setName('psd model');
+
+%% Checking the two fit results
+% plot the normalized difference between fit results and model data in
+% order to check the fit accuracy
+plpl = plist('Legends',{'Fit1','Fit2'});
+iplot(abs(bsp1-sqrt(modr))./sqrt(modr),abs(bsp2-sqrt(modr))./sqrt(modr),plpl)
+