Mercurial > hg > ltpda

comparison m-toolbox/test/LTPDA_training/topic5/TrainigSession_T5_Ex01.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 % Training session Topic 5 exercise 01
2 %
3 % System identification in z-domain 1
4 %
5 % 1) Generate white noise
6 % 2) Filter white noise with a miir filter
7 % 3) Extract transfer function from data
8 % 4) Fit TF with zDomainFit
9 % 5) Compare results
10 % 6) Comments
11 %
12 % L FERRAIOLI 22-02-09
13 %
14 % $Id: TrainigSession_T5_Ex01.m,v 1.2 2009/02/25 18:18:45 luigi Exp $
15 %
16 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17
18 %% 1) Generate White Noise
19
20 % Build a gaussian random noise data series
21 a = ao(plist('tsfcn', 'randn(size(t))', 'fs', 1, 'nsecs', 10000,'yunits','m'));
22
23 %% 2) Filtering White Noise
24
25 % set a proper plist for a miir filter converting meters to volt
26 plmiir = plist('a',[2 0.3 0.4],'b',[1 .05 .3 0.1],'fs',1,'iunits','m','ounits','V');
27
28 % Building the filter
29 filt = miir(plmiir);
30
31 % filtering white noise
32 ac = filter(a,filt);
33 ac.simplifyYunits;
34
35 %% Checking
36 % set a plist for psd calculation
37 plpsd = plist('Nfft',1000);
38
39 % calculate psd of white data
40 axx = psd(a,plpsd);
41
42 % calculate psd of colored data
43 acxx = psd(ac,plpsd);
44
45 % plotting psd
46 iplot(axx,acxx)
47
48 %% 3) Extracting transfer function
49
50 % settung plist for TF estimation
51 pltf = plist('Nfft',1000);
52
53 % TF estimation
54 tf = tfe(a,ac,pltf);
55
56 % Removing first bins form tf
57 tf12 = split(tf(1,2),plist('split_type', 'frequencies', 'frequencies', [3e-3 5e-1]));
58
59 % plotting
60 iplot(tf(1,2),tf12)
61
62 %% 4.1) Fitting TF - checking residuals log difference
63
64 % Perform an automatic fit of TF data in order to identify a discrete
65 % model. The parameters controlling fit accuracy are
66 %
67 % 'ResLogDiff' check that the normalized difference between data and fit
68 % residuals in log scale is larger than the value specified.
69 % if d is the input value for the parameter, than the algorithm check that
70 % log10(data) - log10(res) > d
71 % In other words if d = 1 an accuracy of 10% is guaranteed
72 % if d = an accuracy of 1% is guaranteed
73 % Note1 - fit residuals are calculated as res = model - data
74 % Note2 - In fitting noisy data, this option should be relaxed otherwise
75 % the function try to fit accurately the random oscillations due to the
76 % noise
77 %
78 % 'RMSE' if r is the parameter value, then the algorithm check that the
79 % step-by-step Root Mean Squared Error variation is less than 10^-r
80
81 plfit1 = plist('FS',1,... % Sampling frequency for the model filters
82 'AutoSearch','on',... % Automatically search for a good model
83 'StartPolesOpt','c1',... % Define the properties of the starting poles - complex distributed in the unitary circle
84 'maxiter',50,... % maximum number of iteration per model order
85 'minorder',2,... % minimum model order
86 'maxorder',9,... % maximum model order
87 'weightparam','abs',... % assign weights as 1./abs(data)
88 'ResLogDiff',2,... % Residuals log difference
89 'ResFlat',[],... % Rsiduals spectral flatness (no need to assign this parameter for the moment)
90 'RMSE',7,... % Root Mean Squared Error Variation
91 'Plot','on',... % set the plot on or off
92 'ForceStability','on',... % force to output a stable ploes model
93 'CheckProgress','off'); % display fitting progress on the command window
94
95 % Do the fit
96 fobj1 = zDomainFit(tf12,plfit1);
97 % setting input and output units for fitted model
98 fobj1(:).setIunits('m');
99 fobj1(:).setOunits('V');
100
101 %% 4.2) Fitting TF - checking residuals spectral flatness
102
103 % Perform an automatic fit of TF data in order to identify a discrete
104 % model. The parameters controlling fit accuracy are
105 %
106 % 'ResFlat' Check that normalized fit residuals spectral flatness is larger
107 % than the parameter value. Admitted values for the parameter are 0<r<1.
108 % Suggested values are 0.5<r<0.7.
109 % Note - normalized fit residuals are calculated as
110 % res = (model - data)./data
111 % Note2 - This option could be useful with noisy data, when ResLogDiff
112 % option is not useful
113 %
114 % 'RMSE' if r is the parameter value, then the algorithm check that the
115 % step-by-step Root Mean Squared Error variation is less than 10^-r
116
117 plfit2 = plist('FS',1,... % Sampling frequency for the model filters
118 'AutoSearch','on',... % Automatically search for a good model
119 'StartPolesOpt','c1',... % Define the properties of the starting poles - complex distributed in the unitary circle
120 'maxiter',50,... % maximum number of iteration per model order
121 'minorder',2,... % minimum model order
122 'maxorder',9,... % maximum model order
123 'weightparam','abs',... % assign weights as 1./abs(data)
124 'ResLogDiff',[],... % Residuals log difference (no need to assign this parameter for the moment)
125 'ResFlat',0.65,... % Rsiduals spectral flatness
126 'RMSE',7,... % Root Mean Squared Error Variation
127 'Plot','on',... % set the plot on or off
128 'ForceStability','on',... % force to output a stable ploes model
129 'CheckProgress','off'); % display fitting progress on the command window
130
131 fobj2 = zDomainFit(tf12,plfit2);
132 % setting input and output units for fitted model
133 fobj2(:).setIunits('m');
134 fobj2(:).setOunits('V');
135
136 %% 5) Evaluating accuracy
137
138 % fobj1 and fobj2 are parallel bank of miir filters. Therofore we can
139 % calculate filter response and compare with our starting model
140
141 % set plist for filter response
142 plrsp = plist('bank','parallel','f1',1e-5,'f2',0.5,'nf',100,'scale','log');
143
144 % calculate filters responses
145 rfilt = resp(filt,plrsp); % strating model response
146 rfobj1 = resp(fobj1,plrsp); % result of first fit
147 rfobj2 = resp(fobj2,plrsp); % result of second fit
148
149 % plotting filters
150 iplot(rfilt,rfobj1,rfobj2)
151
152 % plotting difference
153 iplot(abs(rfobj1-rfilt),abs(rfobj2-rfilt))
154
155 %% 6) Comments
156
157 % Decomposition in partial fractions of the starting filter provides
158 % residues =
159 %
160 % 0.531297296250259 - 0.28636710753776i
161 % 0.531297296250259 + 0.28636710753776i
162 % 0.937405407499481
163 %
164 %
165 % poles =
166 %
167 % 0.112984252455744 + 0.591265379634664i
168 % 0.112984252455744 - 0.591265379634664i
169 % -0.275968504911487
170 %
171 % If we look inside the fitted objects fobj1 and fobj2 the following
172 % results can be extracted
173 % rsidues =
174 % 0.531391341497363 - i*0.286395031447695
175 % 0.531391341497363 + i*0.286395031447695
176 % 0.937201773620951
177 %
178 % poles =
179 % 0.112971704620295 + i*0.591204610859687
180 % 0.112971704620295 - i*0.591204610859687
181 % -0.276032231799774
182 %
183 % Model order is correct and residues and poles values are accurate to the
184 % third decimal digit.
185 %
186 % With low order model the algorithm works fine. For higher order model the
187 % algorithm is capable to provide a model within the prescribed accuracy
188 % but the order could be not guaranteed. Lets see what happen with higher
189 % order models opening the second exercise script:
190 % TrainingSession_T5_Ex02.m
191
192