0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % BUILDWHITENER1D builds a whitening filter based on the input frequency-series.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % DESCRIPTION:BUILDWHITENER1D builds a whitening filter based on the input frequency-series.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % The filter is built by fitting to the model provided.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 % If no model is provided, a fit is made to a spectral-density estimate of the
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % input time-series (made using psd+bin_data or lpsd).
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % Note: The function assumes that the input model corresponds
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % to the one-sided psd of the data to be whitened.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 % ALGORITHM:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % 1) If no model provided, make psd+bin_data or lpsd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 % of time-series and take it as a model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 % for the data power spectral density
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 % 2) Fit a set of partial fraction z-domain filters using
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 % utils.math.psd2wf. The fit is automatically stopped when
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 % the accuracy tolerance is reached.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 % 3) Convert to array of MIIR filters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 % 4) Assemble into a parallel filterbank object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 % CALL: b = buildWhitener1D(a, pl)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 % [b1,b2,...,bn] = buildWhitener1D(a1,a2,...,an, pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 % INPUT:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 % - as is a time-series analysis object or a vector of
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 % analysis objects
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 % - pl is a plist with the input parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 % OUTPUT:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 % - b "whitening" filters, stored into a filterbank.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 % <a href="matlab:utils.helper.displayMethodInfo('ao', 'buildWhitener1D')">Parameters Description</a>
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 % VERSION: $Id: buildWhitener1D.m,v 1.11 2011/04/18 19:46:57 mauro Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 function varargout = buildWhitener1D(varargin)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 % Check if this is a call for parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 if utils.helper.isinfocall(varargin{:})
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 varargout{1} = getInfo(varargin{3});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44 return
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 import utils.const.*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 utils.helper.msg(msg.PROC3, 'running %s/%s', mfilename('class'), mfilename);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 % Collect input variable names
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 in_names = cell(size(varargin));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 for ii = 1:nargin,in_names{ii} = inputname(ii);end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 % Collect all AOs and plists
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55 [as, ao_invars] = utils.helper.collect_objects(varargin(:), 'ao', in_names);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 pl = utils.helper.collect_objects(varargin(:), 'plist', in_names);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 if nargout == 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 error('### buildWhitener cannot be used as a modifier. Please give an output variable.');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 % combine plists
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 if isempty(pl)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 model = 'psd';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 model = find(pl, 'model');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 if isempty(model)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68 model = 'psd';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 pl.pset('model', model);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 if ischar(model)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 pl = parse(pl, getDefaultPlist(model));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 pl = parse(pl, getDefaultPlist('Default'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 pl.getSetRandState();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80 % Collect input histories
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 inhists = [as.hist];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 % Initialize output objects
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 bs = filterbank.initObjectWithSize(1, numel(as));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 % Loop over input AOs
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 for jj = 1:numel(as)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 % 1) searching for input model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 switch class(as(jj).data)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 case 'tsdata'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 % Build the model based on input time-series
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 utils.helper.msg(msg.PROC1, 'user input tsdata object, estimating the model from it');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 model = estimateModel(as(jj), pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 case 'fsdata'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 % The input data are the model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 utils.helper.msg(msg.PROC1, 'user input fsdata object, taking it as the model');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 model = as(jj);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 warning('!!! %s expects ao/tsdata or ao/fsdata objects. Skipping AO %s', mfilename, ao_invars{jj});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 return;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 %-------------- Whiten this AO
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 % Extract necessary parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 % Tolerance for MSE Value
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 lrscond = find(pl, 'FITTOL');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 % give an error for strange values of lrscond
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 if lrscond < 0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 error('!!! Negative values for FITTOL are not allowed !!!')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 % handling data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 lrscond = -1 * log10(lrscond);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 % give a warning for strange values of lrscond
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 if lrscond<0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 warning('You are searching for a MSE lower than %s', num2str(10^(-1*lrscond)))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120 params.lrscond = lrscond;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122 % Tolerance for the MSE relative variation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123 msevar = find(pl, 'MSEVARTOL');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 % handling data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 msevar = -1 * log10(msevar);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126 % give a warning for strange values of msevar
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 if msevar<0
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128 warning('You are searching for MSE relative variation lower than %s', num2str(10^(-1*msevar)))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 params.msevar = msevar;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 if isempty(params.msevar)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133 params.ctp = 'chival';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135 params.ctp = 'chivar';
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 % Weights
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139 switch find(pl, 'Weights')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140 case {'equal', 'flat', 1}
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141 params.weightparam = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142 case {'1/abs', '1./abs', 2}
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143 params.weightparam = 2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144 case {'1/abs^2', '1/abs2', '1./abs^2', '1./abs2', 3}
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145 params.weightparam = 3;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147 warning('Unrecognized weights option %s', find(pl, 'Weights'))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150 % 2) Build filters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152 % Build input structure for psd2wf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153 params.idtp = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154 params.Nmaxiter = find(pl, 'MaxIter');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155 params.minorder = find(pl, 'MinOrder');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
156 params.maxorder = find(pl, 'MaxOrder');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
157 params.spolesopt = find(pl, 'PoleType');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
158 params.spy = find(pl, 'Disp');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
159
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
160
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
161 if (find(pl, 'plot'))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
162 params.plot = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
163 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
164 params.plot = 0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
165 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
166
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
167 fs = find(pl, 'fs');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
168 if isempty(fs) || fs <= 0 || ~isfinite(fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
169 if isempty(model.fs) || model.fs <= 0 || ~isfinite(model.fs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
170 error('### Invalid fs value %s. Please specify a meaningful fs, either via the model or in the plist', num2str(fs));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
171 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
172 fs = model.fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
173 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
174 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
175
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
176 params.fs = fs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
177 params.usesym = 0;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
178 params.dterm = 0; % it is better to fit without direct term
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
179 params.fullauto = 1;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
180
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
181 % call psd2wf
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
182 [res, poles, dterm, mresp, rdl] = ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
183 utils.math.psd2wf(model.y,[],[],[],model.x,params);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
184
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
185 % 3) Convert to MIIR filters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
186
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
187 % filtering with a stable model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
188 pfilts = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
189 for kk = 1:numel(res)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
190 ft = miir(res(kk), [ 1 -poles(kk)], fs);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
191 pfilts = [pfilts ft];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
192 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
193
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
194 % 4) Build the output filterbank object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
195 bs(jj) = filterbank(plist('filters', pfilts, 'type', 'parallel'));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
196 % set the input units to be the same as the model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
197 bs(jj).setIunits(sqrt(model.yunits * unit('Hz')));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
198 % set the output units to be empty
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
199 bs(jj).setOunits(unit());
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
200
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
201 % set the name for this object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
202 bs(jj).name = sprintf('buildWhitener1D(%s)', ao_invars{jj});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
203 % add history
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
204 bs(jj).addHistory(getInfo('None'), pl, ao_invars(jj), inhists(jj));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
205 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
206
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
207 % Set output
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
208 if nargout == numel(bs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
209 % List of outputs
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
210 for ii = 1:numel(bs)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
211 varargout{ii} = bs(ii);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
212 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
213 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
214 % Single output
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
215 varargout{1} = bs;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
216 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
217 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
218
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
219 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
220 % Get Info Object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
221 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
222 function ii = getInfo(varargin)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
223 if nargin == 1 && strcmpi(varargin{1}, 'None')
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
224 sets = {};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
225 pl = [];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
226 elseif nargin == 1 && ~isempty(varargin{1}) && ischar(varargin{1})
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
227 sets{1} = varargin{1};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
228 pl = getDefaultPlist(sets{1});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
229 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
230 sets = SETS();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
231 % get plists
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
232 pl(size(sets)) = plist;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
233 for kk = 1:numel(sets)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
234 pl(kk) = getDefaultPlist(sets{kk});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
235 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
236 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
237 % Build info object
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
238 ii = minfo(mfilename, 'ao', 'ltpda', utils.const.categories.sigproc, '$Id: buildWhitener1D.m,v 1.11 2011/04/18 19:46:57 mauro Exp $', sets, pl);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
239 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
240
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
241
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
242 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
243 % Defintion of Sets
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
244 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
245
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
246 function out = SETS()
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
247 out = {...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
248 'Default', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
249 'PSD', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
250 'LPSD' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
251 };
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
252 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
253
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
254 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
255 % Get Default Plist
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
256 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
257 function plout = getDefaultPlist(set)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
258 persistent pl;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
259 persistent lastset;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
260 if ~exist('pl', 'var') || isempty(pl) || ~strcmp(lastset, set)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
261 pl = buildplist(set);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
262 lastset = set;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
263 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
264 plout = pl;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
265 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
266
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
267
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
268 function pl = buildplist(set)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
269 pl = plist();
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
270
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
271 % Model
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
272 p = param({'model', ['A model estimation technique in the case of tsdata input:<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
273 '<li>PSD - using <tt>psd</tt> + <tt>bin_data</tt></li>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
274 '<li>LPSD - using <tt>lpsd</tt></li>']}, {1, {'PSD', 'LPSD'}, paramValue.SINGLE});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
275 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
276
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
277 % Range
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
278 p = param({'range', ['The frequency range to evaluate the fitting.<br>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
279 'An empty value or [-inf inf] will include the whole range.<br>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
280 'The remaining part of the model will be completed according<br>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
281 'to the option chosen in the ''complete'' parameter.<br>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
282 ]}, paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
283 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
284
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
285 % Complete
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
286 p = param({'complete_hf', ['Choose how to complete the frequency range up to fs/2.<ol>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
287 '<li>Assumes flat response</li>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
288 '<li>Assumes 4 poles low-pass type response</li>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
289 ]}, {1,{'flat', 'lowpass'}, paramValue.SINGLE});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
290 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
291
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
292 % fs
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
293 p = param({'fs', ['The sampling frequency to design the output filter on.<br>' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
294 'If it is not a positive number, it will be taken from the model' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
295 ]}, paramValue.EMPTY_DOUBLE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
296 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
297
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
298 % MaxIter
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
299 p = param({'MaxIter', 'Maximum number of iterations in fit routine.'}, paramValue.DOUBLE_VALUE(30));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
300 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
301
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
302 % PoleType
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
303 p = param({'PoleType', ['Choose the pole type for fitting:<ol>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
304 '<li>use real starting poles</li>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
305 '<li>generates complex conjugate poles of the<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
306 'type <tt>a.*exp(theta*pi*j)</tt>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
307 'with <tt>theta = linspace(0,pi,N/2+1)</tt></li>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
308 '<li>generates complex conjugate poles of the type<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
309 '<tt>a.*exp(theta*pi*j)</tt><br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
310 'with <tt>theta = linspace(0,pi,N/2+2)</tt></li></ol>']}, {1, {1, 2, 3}, paramValue.SINGLE});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
311 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
312
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
313 % MinOrder
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
314 p = param({'MinOrder', 'Minimum order to fit with.'}, paramValue.DOUBLE_VALUE(2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
315 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
316
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
317 % MaxOrder
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
318 p = param({'MaxOrder', 'Maximum order to fit with.'}, paramValue.DOUBLE_VALUE(25));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
319 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
320
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
321 % Weights
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
322 p = param({'Weights', ['Choose weighting method:<ol>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
323 '<li>equal weights for each point</li>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
324 '<li>weight with <tt>1/abs(model)</tt></li>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
325 '<li>weight with <tt>1/abs(model).^2</tt></li></ol>']}, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
326 {2, {'equal', '1/abs', '1/abs^2'}, paramValue.SINGLE});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
327 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
328
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
329 % Plot
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
330 p = param({'Plot', 'Plot results of each fitting step.'}, paramValue.FALSE_TRUE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
331 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
332
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
333 % Disp
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
334 p = param({'Disp', 'Display the progress of the fitting iteration.'}, paramValue.FALSE_TRUE);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
335 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
336
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
337 % MSEVARTOL
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
338 p = param({'MSEVARTOL', ['Mean Squared Error Variation - Check if the<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
339 'relative variation of the mean squared error is<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
340 'smaller than the value specified. This<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
341 'option is useful for finding the minimum of Chi-squared.']}, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
342 paramValue.DOUBLE_VALUE(1e-1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
343 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
344
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
345 % FITTOL
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
346 p = param({'FITTOL', ['Mean Squared Error Value - Check if the mean<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
347 'squared error value is lower than the value<br>'...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
348 'specified.']}, paramValue.DOUBLE_VALUE(1e-2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
349 pl.append(p);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
350
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
351 % Append sets of parameters according to the chosen spectral estimator
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
352 if ~utils.helper.ismember(lower(SETS), lower(set))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
353 error('### Unknown set [%s]', set);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
354 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
355
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
356 switch lower(set)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
357 case 'default'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
358 pl.remove('model');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
359 case 'psd'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
360 pl = combine(pl, ao.getInfo('psd').plists);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
361 pl.pset(...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
362 'model', 'PSD', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
363 'Navs', 16, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
364 'order', 1, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
365 'olap', 50 ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
366 );
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
367 pl = combine(pl, ao.getInfo('bin_data').plists);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
368 pl.pset(...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
369 'method', 'MEAN', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
370 'resolution', 50 ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
371 );
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
372 case 'lpsd'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
373 pl = combine(pl, ao.getInfo('lpsd').plists);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
374 pl.pset(...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
375 'model', 'LPSD' ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
376 );
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
377 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
378 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
379 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
380
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
381
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
382 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
383 % Estimate a model from the data or from user input
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
384 %--------------------------------------------------------------------------
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
385 function model = estimateModel(b, pl)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
386
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
387 import utils.const.*
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
388
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
389
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
390 % Estimate a model for the PSD
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
391 model_all = find(pl, 'model');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
392 if ischar(model_all)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
393 switch lower(model_all)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
394 case 'psd'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
395 % Select only the parameters associated to ao/psd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
396 pls = ao.getInfo('psd').plists;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
397 % Call ao/psd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
398 sp = psd(b, pl.subset(pls.getKeys()));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
399 % Select only the parameters associated to ao/bin_data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
400 pls = ao.getInfo('bin_data').plists;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
401 % Call ao/bin_data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
402 model_all = bin_data(sp, pl.subset(pls.getKeys()));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
403 case 'lpsd'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
404 % Select only the parameters associated to ao/lpsd
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
405 pls = ao.getInfo('lpsd').plists;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
406 model_all = lpsd(b, pl.subset(pls.getKeys()));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
407 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
408 error('### Unknown model [%s]', model_all);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
409 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
410 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
411
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
412 % Select only a limited frequency range
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
413 frange = find(pl, 'range');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
414 if isempty(frange)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
415 frange = [-inf inf];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
416 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
417 model = model_all.split(plist('frequencies', frange));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
418 f1 = frange(1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
419 f2 = frange(2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
420
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
421 if isfinite(f2)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
422 % Select a technique to complete the high frequency range
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
423 complete_up_opt = find(pl, 'complete_hf');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
424 switch complete_up_opt
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
425 case {'flat', 'allpass', 'all pass', 'all-pass'}
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
426 utils.helper.msg(msg.PROC1, 'Completing the frequency range from %s to %s with flat model', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
427 num2str(f2), num2str(b.fs/2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
428 % Build a flat model response
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
429 r = ones(size(model_all.x));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
430 case {'lowpass', 'low pass', 'low-pass'}
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
431 utils.helper.msg(msg.PROC1, 'Completing the frequency range from %s to %s with 4 poles low-pass model', ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
432 num2str(f2), num2str(b.fs/2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
433 % Build a 4 poles low pass resp
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
434 r = abs(resp(pzmodel(plist('gain', 1, 'poles', {10*f2,11*f2,12*f2,13*f2})), plist('f', model_all.x)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
435 r = r.y;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
436 otherwise
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
437 error('### Unknown option [%s] for high frequency completion', complete_up_opt);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
438 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
439 model_hf = r * model.y(end);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
440 model = join(model, ...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
441 ao(plist('type', 'fsdata', 'xvals', model_all.x, 'yvals', model_hf, 'fs', model_all.fs, 'yunits', model_all.yunits)));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
442 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
443
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
444 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
445
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
446
|