Mercurial > hg > ltpda

comparison m-toolbox/test/test_ao_zDomainFit_3.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % Test script for zDomainFit on noisy data
2 %
3 % L. Ferraioli 02-12-08
4 %
5 % $Id: test_ao_zDomainFit_3.m,v 1.3 2009/04/20 14:31:14 luigi Exp $
6 %
7 %
8 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
9 % Smooth psd model
10 %
11 % We define a smooth psd model and generate colored noise with such a psd
12 % shape. Then we fit sqrt(psd) of noisy data following two different
13 % procedure to decide when the fit is satisfactory and we can exit from the
14 % fitting routine.
15 % The first procedure is based on the check of residuals spectral flatness.
16 % The second procedure is based on the check of residuals log distance to
17 % fitted data.
18 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
19 %% Useful constants
20
21 fs = 10; % sampling frequency Hz
22 Nsecs = 1e5; % Number of seconds in the data stream
23 Nfft = 1e4; % Number of samples in the fft
24 pls = plist('Nfft', Nfft,'Order',0); % plist for spectra
25
26 %% Building a frequency response AO
27
28 % Create a frequency-series AO
29 % pl_data = plist('fsfcn', '0.01./(0.01+f)', 'f1', 1e-6, 'f2', 5, 'nf', 1000);
30 pl_data = plist('fsfcn', '(1e-3./(f).^2 + 1e3./(0.001+f) + 1e5.*f.^2).*1e-10', 'f1', 1e-6, 'f2', 5, 'nf', 300);
31 mod = ao(pl_data);
32 mod.setName('psd model');
33
34 iplot(mod)
35
36 %% Building white noise
37
38 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs));
39 a.setName('white noise')
40 a.setYunits('m')
41
42 %% Calling the noise generator
43
44 pl = plist(...
45 'model', mod, ...
46 'MaxIter', 50, ...
47 'PoleType', 2, ...
48 'MinOrder', 15, ...
49 'MaxOrder', 25, ...
50 'Weights', 2, ...
51 'Plot', false,...
52 'Disp', false,...
53 'RMSEVar', 8,...
54 'FitTolerance', 2);
55
56 ac = noisegen1D(a, pl);
57
58 %% Spetrum and Transfer function
59
60 acxx = ac.psd(pls); % spectrum of data
61 etf = tfe(a,ac,pls); % calculating transfer functions from data
62
63 iplot(acxx)
64 iplot(etf(1,2))
65
66 %% Fitting Sqrt Spectrum - Residuals spectral flatness checking
67
68 % Fitting parameter list
69 pl_fit1 = plist('FS',fs,...
70 'AutoSearch','on',...
71 'StartPoles',[],...
72 'StartPolesOpt','c2',...
73 'maxiter',60,...
74 'minorder',19,...
75 'maxorder',25,...
76 'weights',[],...
77 'weightparam','abs',...
78 'ResLogDiff',[],...
79 'ResFlat',0.7,...
80 'RMSE',5,...
81 'Plot','off',...
82 'ForceStability','off',...
83 'CheckProgress','on');
84
85 % Do fit
86 [filt2,bsp1,bspres1,bsprmse1] = zDomainFit(sqrt(acxx./2), pl_fit1); % division by 2 is necessary because acxx is the one-sided psd
87 bsp1.setName('Fit1');
88
89 %% Checing spectral flatness
90 % checking the spectral flatness of the normalized residuals in order to be
91 % sure that the procedure works correctly
92 kk = abs(bspres1./bsp1);
93 rf = utils.math.spflat(kk.data.y);
94 iplot(abs(bspres1./bsp1))
95
96 %% Fitting Spectrum - Residuals log difference checking
97
98 % Fitting parameter list
99 pl_fit2 = plist('FS',fs,...
100 'AutoSearch','on',...
101 'StartPoles',[],...
102 'StartPolesOpt','c2',...
103 'maxiter',60,...
104 'minorder',17,...
105 'maxorder',20,...
106 'weights',[],...
107 'weightparam','abs',...
108 'ResLogDiff',0.5,...
109 'ResFlat',0.7,...
110 'RMSE',5,...
111 'Plot','off',...
112 'ForceStability','off',...
113 'CheckProgress','on');
114
115 % Do fit
116 [filt2,bsp2,bspres2,bsprmse2] = zDomainFit(sqrt(acxx./2), pl_fit2); % division by 2 is necessary because acxx is the one-sided psd
117 bsp2.setName('Fit2');
118
119 %% Model calculated at the same frequencies of the fit AOs
120
121 % This step is useful in order to compare the results of the two fitting
122 % procedures
123 f = acxx.data.x;
124 pl_datar = plist('fsfcn', '(1e-3./(f).^2 + 1e3./(0.001+f) + 1e5.*f.^2).*1e-10', 'f', f);
125 modr = ao(pl_datar);
126 modr.setName('psd model');
127
128 %% Checking the two fit results
129 % plot the normalized difference between fit results and model data in
130 % order to check the fit accuracy
131 plpl = plist('Legends',{'Fit1','Fit2'});
132 iplot(abs(bsp1-sqrt(modr))./sqrt(modr),abs(bsp2-sqrt(modr))./sqrt(modr),plpl)
133
134 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
135
136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
137 % Peaked psd model
138 %
139 % We define a psd model with some resonance peaks and generate colored
140 % noise with such a psd shape. Then we fit sqrt(psd) of noisy data
141 % following two different procedure to decide when the fit is satisfactory
142 % and we can exit from the fitting routine.
143 % The first procedure is based on the check of residuals spectral flatness.
144 % The second procedure is based on the check of residuals log distance to
145 % fitted data.
146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
147
148 %% Building a frequency response AO
149
150 % Create a psd model with two peak resonance at 5e-2 and 7e-2 Hz
151 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';
152 pl_data = plist('fsfcn', func, 'f1', 1e-6, 'f2', 5, 'nf', 300);
153 mod = ao(pl_data);
154 mod.setName('psd model');
155
156 iplot(mod)
157
158 %% Building white noise
159
160 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', fs, 'nsecs', Nsecs, 'yunits', 'm'));
161 a.setName('white noise')
162
163 %% Calling the noise generator
164
165 pl = plist(...
166 'model', mod, ... % Multiplication by fs needed to preserve energy
167 'MaxIter', 70, ...
168 'PoleType', 2, ...
169 'MinOrder', 30, ...
170 'MaxOrder', 45, ...
171 'Weights', 2, ...
172 'Plot', false,...
173 'Disp', false,...
174 'RMSEVar', 8,...
175 'FitTolerance', 2);
176
177 ac = noisegen1D(a, pl);
178
179 %% Spetrum and Transfer function
180
181 acxx = ac.psd(pls); % spectrum of data
182
183 iplot(acxx)
184
185 %% Fitting Sqrt Spectrum - Residuals spectral flatness checking
186
187 % Fitting parameter list
188 pl_fit1 = plist('FS',fs,...
189 'AutoSearch','on',...
190 'StartPoles',[],...
191 'StartPolesOpt','c2',...
192 'maxiter',60,...
193 'minorder',15,...
194 'maxorder',25,...
195 'weights',[],...
196 'weightparam','abs',...
197 'ResLogDiff',[],...
198 'ResFlat',0.7,...
199 'RMSE',6,...
200 'Plot','off',...
201 'ForceStability','off',...
202 'CheckProgress','on');
203
204 % Do fit
205 [filt1,bsp1,bspres1,bsprmse1] = zDomainFit(sqrt(acxx./2), pl_fit1); % division by 2 is necessary because acxx is the one-sided psd
206 bsp1.setName('Fit1');
207
208 %%
209
210 iplot(sqrt(acxx./2),bsp1)
211
212 %% Checing spectral flatness
213 % checking the spectral flatness of the normalized residuals in order to be
214 % sure that the procedure works correctly
215 kk = abs(bspres1./bsp1);
216 rf = utils.math.spflat(kk.data.y);
217 iplot(abs(bspres1./bsp1))
218
219 %% Fitting Spectrum - Residuals log difference checking
220
221 % Fitting parameter list
222 pl_fit2 = plist('FS',fs,...
223 'AutoSearch','on',...
224 'StartPoles',[],...
225 'StartPolesOpt','c2',...
226 'maxiter',60,...
227 'minorder',15,...
228 'maxorder',25,...
229 'weights',[],...
230 'weightparam','abs',...
231 'ResLogDiff',0.5,...
232 'ResFlat',0.7,...
233 'RMSE',6,...
234 'Plot','off',...
235 'ForceStability','off',...
236 'CheckProgress','on');
237
238 % Do fit
239 [filt2,bsp2,bspres2,bsprmse2] = zDomainFit(sqrt(acxx./2), pl_fit2); % division by 2 is necessary because acxx is the one-sided psd
240 bsp2.setName('Fit2');
241
242 %%
243
244 iplot(sqrt(acxx./2),bsp2)
245
246 %% Model calculated at the same frequencies of the fit AOs
247
248 % This step is useful in order to compare the results of the two fitting
249 % procedures
250 f = acxx.data.x;
251 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';
252 pl_datar = plist('fsfcn', func, 'f', f);
253 modr = ao(pl_datar);
254 modr.setName('psd model');
255
256 %% Checking the two fit results
257 % plot the normalized difference between fit results and model data in
258 % order to check the fit accuracy
259 plpl = plist('Legends',{'Fit1','Fit2'});
260 iplot(abs(bsp1-sqrt(modr))./sqrt(modr),abs(bsp2-sqrt(modr))./sqrt(modr),plpl)
261